Switch covariance estimation to numpy.ma.cov

This may give slightly different results from `pd.DataFrame.cov` in
the presence of missing data, because the calculation of the mean is
done differently, see pandas-dev/pandas#16837
As both implementations are correct, we switch to the more performant
one here

Author: HippocampusGirl

src/gwas/src/gwas/pheno.py

@@ -348,24 +348,24 @@ def from_txt(
 
     def covariance_to_txt(
         self, path: UPath, compression_method: CompressionMethod, num_threads: int
-    ) -> None:
+    ) -> UPath:
         path = path.with_suffix(compression_method.suffix)
         if path.is_file():
             logger.debug("Skip writing covariance matrix because it already exists")
-            return
+            return path
 
         phenotype_array = self.phenotypes.to_numpy()
         covariate_array = self.covariates.to_numpy()
 
         array = np.hstack((phenotype_array, covariate_array))
         names = [*self.phenotype_names, *self.covariate_names]
 
-        data_frame = pd.DataFrame(array, index=self.samples, columns=names)
-
         logger.debug("Calculating covariance matrix")
-        covariance: npt.NDArray[np.float64] = np.asfortranarray(
-            data_frame.cov().to_numpy(dtype=np.float64)
-        )
+        masked_array = np.ma.array(array, mask=np.isnan(array))
+        numpy_covariance = np.ma.cov(
+            masked_array, rowvar=False, allow_masked=True
+        ).filled(np.nan)
+        covariance: npt.NDArray[np.float64] = np.asfortranarray(numpy_covariance)
 
         writer: FileArrayWriter[np.float64] = FileArray.create(
             path,
@@ -382,6 +382,8 @@ def covariance_to_txt(
         with writer:
             writer[:, :] = covariance
 
+        return writer.file_path
+
 
 @dataclass
 class VariableSummary:

src/gwas/tests/score/test_pheno.py

@@ -8,7 +8,10 @@
 from pytest import FixtureRequest
 from upath import UPath
 
-from gwas.compression.arr.base import Blosc2CompressionMethod, compression_methods
+from gwas.compression.arr.base import (
+    Blosc2CompressionMethod,
+    compression_methods,
+)
 from gwas.mem.wkspace import SharedWorkspace
 from gwas.pheno import VariableCollection
 from gwas.utils import cpu_count
@@ -209,7 +212,7 @@ def test_covariance(
     request.addfinalizer(variable_collection.free)
 
     covariance_path = tmp_path / "covariance.tsv"
-    variable_collection.covariance_to_txt(
+    covariance_path = variable_collection.covariance_to_txt(
         covariance_path, compression_method, num_threads=cpu_count()
     )