all unit tests are working

fixed remaining interfaces to batchglm

Author: davidsebfischer

diffxpy/testing/base.py

@@ -53,6 +53,16 @@ def design_loc(self) -> np.ndarray:
     def design_scale(self) -> np.ndarray:
         pass
 
+    @property
+    @abc.abstractmethod
+    def constraints_loc(self) -> np.ndarray:
+        pass
+
+    @property
+    @abc.abstractmethod
+    def constraints_scale(self) -> np.ndarray:
+        pass
+
     @property
     @abc.abstractmethod
     def num_observations(self) -> int:
@@ -73,12 +83,9 @@ def features(self) -> np.ndarray:
     def observations(self) -> np.ndarray:
         pass
 
+    @property
     @abc.abstractmethod
-    def probs(self) -> np.ndarray:
-        pass
-
-    @abc.abstractmethod
-    def log_probs(self) -> np.ndarray:
+    def log_likelihood(self, **kwargs) -> np.ndarray:
         pass
 
     @property
@@ -88,7 +95,7 @@ def loss(self, **kwargs) -> np.ndarray:
 
     @property
     @abc.abstractmethod
-    def gradient(self, **kwargs) -> np.ndarray:
+    def gradients(self, **kwargs) -> np.ndarray:
         pass
 
     @property
@@ -410,24 +417,21 @@ def X(self):
 
     @property
     def reduced_model_gradient(self):
-        return self.reduced_estim.gradient
+        return self.reduced_estim.gradients
 
     @property
     def full_model_gradient(self):
-        return self.full_estim.gradient
+        return self.full_estim.gradients
 
     def _test(self):
-        full = np.sum(self.full_estim.log_probs(), axis=0)
-        reduced = np.sum(self.reduced_estim.log_probs(), axis=0)
-
-        if np.any(full < reduced):
-            logger.warning("Test assumption failed: full model is (partially) less probable than reduced model!")
+        if np.any(self.full_estim.log_likelihood < self.reduced_estim.log_likelihood):
+            logger.warning("Test assumption failed: full model is (partially) less probable than reduced model")
 
         return stats.likelihood_ratio_test(
-            ll_full=full,
-            ll_reduced=reduced,
-            df_full=self.full_estim.design_loc.shape[-1] + self.full_estim.design_scale.shape[-1],
-            df_reduced=self.reduced_estim.design_loc.shape[-1] + self.reduced_estim.design_scale.shape[-1],
+            ll_full=self.full_estim.log_likelihood,
+            ll_reduced=self.reduced_estim.log_likelihood,
+            df_full=self.full_estim.constraints_loc.shape[1] + self.full_estim.constraints_scale.shape[1],
+            df_reduced=self.reduced_estim.constraints_loc.shape[1] + self.reduced_estim.constraints_scale.shape[1],
         )
 
     def _ave(self):
@@ -643,7 +647,7 @@ def X(self):
 
     @property
     def model_gradient(self):
-        return self.model_estim.gradient
+        return self.model_estim.gradients
 
     def log_fold_change(self, base=np.e, **kwargs):
         """
@@ -1244,7 +1248,7 @@ def log_probs(self):
 
     @property
     def model_gradient(self):
-        return self.model_estim.gradient
+        return self.model_estim.gradients
 
     def _ave(self):
         """
@@ -1449,7 +1453,7 @@ def log_probs(self):
 
     @property
     def model_gradient(self):
-        return self.model_estim.gradient
+        return self.model_estim.gradients
 
     def _ave(self):
         """
@@ -2611,6 +2615,9 @@ def lrt(
     """
     Perform log-likelihood ratio test for differential expression for each gene.
 
+    Note that lrt() does not support constraints in its current form. Please
+    use wald() for constraints.
+
     :param data: input data
     :param reduced_formula: formula
         Reduced model formula for location and scale parameter models.
@@ -2734,6 +2741,8 @@ def lrt(
         data=X,
         design_loc=reduced_design_loc,
         design_scale=reduced_design_scale,
+        constraints_loc=None,
+        constraints_scale=None,
         init_a=init_a,
         init_b=init_b,
         as_numeric=as_numeric,
@@ -2750,6 +2759,8 @@ def lrt(
         data=X,
         design_loc=full_design_loc,
         design_scale=full_design_scale,
+        constraints_loc=None,
+        constraints_scale=None,
         gene_names=gene_names,
         init_a="init_model",
         init_b="init_model",
@@ -4035,7 +4046,7 @@ def continuous_1d(
         noise_model: str = 'nb',
         size_factors: np.ndarray = None,
         batch_size: int = None,
-        training_strategy: Union[str, List[Dict[str, object]], Callable] = "CONTINUOUS",
+        training_strategy: Union[str, List[Dict[str, object]], Callable] = "DEFAULT",
         quick_scale: bool = None,
         dtype="float64",
         **kwargs

diffxpy/testing/continuous.py

@@ -75,11 +75,10 @@ def test_null_distribution_wald_constrained(self, n_genes: int = 100):
         # Compare p-value distribution under null model against uniform distribution.
         pval_h0 = stats.kstest(test.pval, 'uniform').pvalue
 
-        print('KS-test pvalue for null model match of wald(): %f' % pval_h0)
-
+        logging.getLogger("diffxpy").info('KS-test pvalue for null model match of wald(): %f' % pval_h0)
         assert pval_h0 > 0.05, "KS-Test failed: pval_h0 is <= 0.05!"
 
-        return pval_h0
+        return True
 
     def test_null_distribution_wald_constrained_2layer(self, n_genes: int = 100):
         """
@@ -162,11 +161,10 @@ def test_null_distribution_wald_constrained_2layer(self, n_genes: int = 100):
         # Compare p-value distribution under null model against uniform distribution.
         pval_h0 = stats.kstest(test.pval, 'uniform').pvalue
 
-        print('KS-test pvalue for null model match of wald(): %f' % pval_h0)
-
+        logging.getLogger("diffxpy").info('KS-test pvalue for null model match of wald(): %f' % pval_h0)
         assert pval_h0 > 0.05, "KS-Test failed: pval_h0 is <= 0.05!"
 
-        return pval_h0
+        return True
 
     def test_null_distribution_wald_multi_constrained_2layer(self, n_genes: int = 50):
         """
@@ -181,8 +179,8 @@ def test_null_distribution_wald_multi_constrained_2layer(self, n_genes: int = 50
         :param n_genes: Number of genes to simulate (number of tests).
         """
         logging.getLogger("tensorflow").setLevel(logging.ERROR)
-        logging.getLogger("batchglm").setLevel(logging.INFO)
-        logging.getLogger("diffxpy").setLevel(logging.INFO)
+        logging.getLogger("batchglm").setLevel(logging.WARNING)
+        logging.getLogger("diffxpy").setLevel(logging.WARNING)
 
         n_cells = 3000
 
@@ -239,11 +237,10 @@ def test_null_distribution_wald_multi_constrained_2layer(self, n_genes: int = 50
         # Compare p-value distribution under null model against uniform distribution.
         pval_h0 = stats.kstest(test.pval, 'uniform').pvalue
 
-        print('KS-test pvalue for null model match of wald(): %f' % pval_h0)
-
+        logging.getLogger("diffxpy").info('KS-test pvalue for null model match of wald(): %f' % pval_h0)
         assert pval_h0 > 0.05, "KS-Test failed: pval_h0 is <= 0.05!"
 
-        return pval_h0
+        return True
 
 
 if __name__ == '__main__':

diffxpy/unit_test/test_constrained.py

@@ -3,11 +3,9 @@
 import numpy as np
 import pandas as pd
 import scipy.stats as stats
-import scipy.sparse
-import anndata
 import logging
 
-from batchglm.api.models.glm_nb import Simulator, Estimator, InputData
+from batchglm.api.models.glm_nb import Simulator
 import diffxpy.api as de
 
 
@@ -20,7 +18,9 @@ def test_forfatal_functions(self):
         :param n_cells: Number of cells to simulate (number of observations per test).
         :param n_genes: Number of genes to simulate (number of tests).
         """
-        logging.getLogger('diffxpy').addFilter('DEBUG')
+        logging.getLogger("tensorflow").setLevel(logging.ERROR)
+        logging.getLogger("batchglm").setLevel(logging.WARNING)
+        logging.getLogger("diffxpy").setLevel(logging.WARNING)
 
         num_observations = 10
         num_features = 2
@@ -45,6 +45,7 @@ def test_forfatal_functions(self):
             sample_description=random_sample_description,
             quick_scale=True,
             batch_size=None,
+            training_strategy="DEFAULT",
             dtype="float64"
         )
 
@@ -67,8 +68,9 @@ def test_forfatal_functions(self):
         temp = test.argmin(genes=ids, nonnumeric=True)
         temp = test.summary(nonnumeric=True)
 
+        return True
 
-    def test_null_distribution_wald(self, n_cells: int = 2000, n_genes: int = 500):
+    def test_null_distribution_wald(self, n_cells: int = 2000, n_genes: int = 100):
         """
         Test if de.test.continuous() generates a uniform p-value distribution in the wald test
         if it is given data simulated based on the null model. Returns the p-value
@@ -78,7 +80,9 @@ def test_null_distribution_wald(self, n_cells: int = 2000, n_genes: int = 500):
         :param n_cells: Number of cells to simulate (number of observations per test).
         :param n_genes: Number of genes to simulate (number of tests).
         """
-        logging.getLogger('diffxpy').addFilter('DEBUG')
+        logging.getLogger("tensorflow").setLevel(logging.INFO)
+        logging.getLogger("batchglm").setLevel(logging.INFO)
+        logging.getLogger("diffxpy").setLevel(logging.WARNING)
 
         sim = Simulator(num_observations=n_cells, num_features=n_genes)
         sim.generate_sample_description(num_batches=0, num_conditions=0)
@@ -99,20 +103,20 @@ def test_null_distribution_wald(self, n_cells: int = 2000, n_genes: int = 500):
             sample_description=random_sample_description,
             quick_scale=True,
             batch_size=None,
+            training_strategy="DEFAULT",
             dtype="float64"
         )
         summary = test.summary()
 
         # Compare p-value distribution under null model against uniform distribution.
         pval_h0 = stats.kstest(test.pval, 'uniform').pvalue
 
-        print('KS-test pvalue for null model match of wald(): %f' % pval_h0)
-
+        logging.getLogger("diffxpy").info('KS-test pvalue for null model match of wald(): %f' % pval_h0)
         assert pval_h0 > 0.05, "KS-Test failed: pval_h0 is <= 0.05!"
 
-        return pval_h0
+        return True
 
-    def test_null_distribution_lrt(self, n_cells: int = 2000, n_genes: int = 500):
+    def test_null_distribution_lrt(self, n_cells: int = 2000, n_genes: int = 100):
         """
         Test if de.test.continuous() generates a uniform p-value distribution in lrt
         if it is given data simulated based on the null model. Returns the p-value
@@ -122,7 +126,9 @@ def test_null_distribution_lrt(self, n_cells: int = 2000, n_genes: int = 500):
         :param n_cells: Number of cells to simulate (number of observations per test).
         :param n_genes: Number of genes to simulate (number of tests).
         """
-        logging.getLogger('diffxpy').addFilter('DEBUG')
+        logging.getLogger("tensorflow").setLevel(logging.INFO)
+        logging.getLogger("batchglm").setLevel(logging.INFO)
+        logging.getLogger("diffxpy").setLevel(logging.WARNING)
 
         sim = Simulator(num_observations=n_cells, num_features=n_genes)
         sim.generate_sample_description(num_batches=0, num_conditions=0)
@@ -143,18 +149,18 @@ def test_null_distribution_lrt(self, n_cells: int = 2000, n_genes: int = 500):
             sample_description=random_sample_description,
             quick_scale=False,
             batch_size=None,
+            training_strategy="DEFAULT",
             dtype="float64"
         )
         summary = test.summary()
 
         # Compare p-value distribution under null model against uniform distribution.
         pval_h0 = stats.kstest(test.pval, 'uniform').pvalue
 
-        print('KS-test pvalue for null model match of wald(): %f' % pval_h0)
-        
+        logging.getLogger("diffxpy").info('KS-test pvalue for null model match of wald(): %f' % pval_h0)
         assert pval_h0 > 0.05, "KS-Test failed: pval_h0 is <= 0.05!"
 
-        return pval_h0
+        return True
 
 if __name__ == '__main__':
     unittest.main()

diffxpy/unit_test/test_continuous.py

@@ -0,0 +1,16 @@
+import unittest
+import numpy as np
+import pandas as pd
+import scipy.stats as stats
+import logging
+
+from batchglm.api.models.glm_nb import Simulator, Estimator, InputData
+import diffxpy.api as de
+
+
+
+class TestCorrection(unittest.TestCase):
+    pass
+
+if __name__ == '__main__':
+    unittest.main()

diffxpy/unit_test/test_correction.py

@@ -0,0 +1,61 @@
+import unittest
+import logging
+
+import numpy as np
+import pandas as pd
+import scipy.stats as stats
+import scipy.sparse
+import anndata
+
+from batchglm.api.models.glm_nb import Simulator
+import diffxpy.api as de
+
+
+class TestDataTypes(unittest.TestCase):
+
+    def test_sparse_anndata(self, n_cells: int = 2000, n_genes: int = 100):
+        """
+        Test if de.wald() generates a uniform p-value distribution
+        if it is given data simulated based on the null model. Returns the p-value
+        of the two-side Kolmgorov-Smirnov test for equality of the observed
+        p-value distribution and a uniform distribution.
+
+        :param n_cells: Number of cells to simulate (number of observations per test).
+        :param n_genes: Number of genes to simulate (number of tests).
+        """
+        logging.getLogger("tensorflow").setLevel(logging.ERROR)
+        logging.getLogger("batchglm").setLevel(logging.WARNING)
+        logging.getLogger("diffxpy").setLevel(logging.WARNING)
+
+        sim = Simulator(num_observations=n_cells, num_features=n_genes)
+        sim.generate_sample_description(num_batches=0, num_conditions=0)
+        sim.generate()
+
+        random_sample_description = pd.DataFrame({
+            "condition": np.random.randint(2, size=sim.num_observations)
+        })
+
+        adata = anndata.AnnData(scipy.sparse.csr_matrix(sim.X.values))
+        # X = adata.X
+        test = de.test.wald(
+            data=adata,
+            factor_loc_totest="condition",
+            formula="~ 1 + condition",
+            sample_description=random_sample_description,
+            quick_scale=True,
+            training_strategy="DEFAULT",
+            dtype="float64"
+        )
+        summary = test.summary()
+
+        # Compare p-value distribution under null model against uniform distribution.
+        pval_h0 = stats.kstest(test.pval, 'uniform').pvalue
+
+        logging.getLogger("diffxpy").info('KS-test pvalue for null model match of wald(): %f' % pval_h0)
+        assert pval_h0 > 0.05, "KS-Test failed: pval_h0 is <= 0.05!"
+
+        return True
+
+
+if __name__ == '__main__':
+    unittest.main()