Merge pull request #79 from theislab/dev

Dev

Author: davidsebfischer

batchglm/api/data.py

@@ -1,9 +1,5 @@
-from batchglm.data import design_matrix
-from batchglm.data import design_matrix_from_xarray
-from batchglm.data import design_matrix_from_anndata
-from batchglm.data import sample_description_from_xarray
-from batchglm.data import sample_description_from_anndata
-from batchglm.data import load_mtx_to_adata
-from batchglm.data import load_mtx_to_xarray
-from batchglm.data import load_recursive_mtx
-from batchglm.data import xarray_from_data
+from batchglm.data import design_matrix, design_matrix_from_xarray, design_matrix_from_anndata
+from batchglm.data import sample_description_from_xarray, sample_description_from_anndata
+from batchglm.data import load_mtx_to_adata, load_mtx_to_xarray, load_recursive_mtx, xarray_from_data
+from batchglm.data import constraint_matrix_from_dict, constraint_matrix_from_string, string_constraints_from_dict
+from batchglm.data import view_coef_names, preview_coef_names

batchglm/data.py

@@ -145,6 +145,7 @@ def _plot_coef_vs_ref(
         :param ncols: Number of columns in plot grid if multiple genes are plotted.
         :param row_gap: Vertical gap between panel rows relative to panel height.
         :param col_gap: Horizontal gap between panel columns relative to panel width.
+        :param title: Plot title.
         :param return_axs: Whether to return axis objects.
         :return: Matplotlib axis objects.
         """
@@ -196,8 +197,8 @@ def _plot_coef_vs_ref(
                 legend=False
             )
             sns.lineplot(
-                x=np.array([np.min([np.min(x), np.min(y), np.max([np.max(x), np.max(y)])])]),
-                y=np.array([np.min([np.min(x), np.min(y), np.max([np.max(x), np.max(y)])])]),
+                x=np.array([np.min([np.min(x), np.min(y)]), np.max([np.max(x), np.max(y)])]),
+                y=np.array([np.min([np.min(x), np.min(y)]), np.max([np.max(x), np.max(y)])]),
                 ax=ax
             )
 

batchglm/models/base/estimator.py

@@ -22,20 +22,20 @@
 TF_CONFIG_PROTO.gpu_options.allow_growth = True
 TF_CONFIG_PROTO.graph_options.optimizer_options.global_jit_level = tf.OptimizerOptions.ON_1
 
-TF_CONFIG_PROTO.inter_op_parallelism_threads = 0 if TF_NUM_THREADS == 0 else 1
+TF_CONFIG_PROTO.inter_op_parallelism_threads = TF_NUM_THREADS
 TF_CONFIG_PROTO.intra_op_parallelism_threads = TF_NUM_THREADS
 
 if TF_NUM_THREADS == 0:
     TF_NUM_THREADS = multiprocessing.cpu_count()
 
 # Trust region hyper parameters:
-TRUST_REGION_RADIUS_INIT = 4.
+TRUST_REGION_RADIUS_INIT = 100.
 TRUST_REGION_ETA0 = 0.
 TRUST_REGION_ETA1 = 0.25
 TRUST_REGION_ETA2 = 0.25  # Allow expansion if not shrinking.
-TRUST_REGION_T1 = 0.1  # Fast collapse to avoid trailing.
-TRUST_REGION_T2 = 2.  # Very conservative expansion to run updates once valid region is reached.
-TRUST_REGION_UPPER_BOUND = 1e4  # Low upper limit so that collapse to valid region does not cause feature to trail.
+TRUST_REGION_T1 = 0.01  # Fast collapse to avoid trailing.
+TRUST_REGION_T2 = 10.
+TRUST_REGION_UPPER_BOUND = 1e5
 
 # Convergence hyper-parameters:
 LLTOL_BY_FEATURE = 1e-10

batchglm/pkg_constants.py

@@ -236,7 +236,7 @@ def train(
         """
         if train_loc is None:
             # check if mu was initialized with MLE
-            train_mu = self._train_loc
+            train_loc = self._train_loc
         if train_scale is None:
             # check if r was initialized with MLE
             train_scale = self._train_scale

batchglm/train/tf/base_glm_all/estimator.py

@@ -11,8 +11,6 @@
 from .model import ProcessModel
 from .training_strategies import TrainingStrategies
 
-logger = logging.getLogger("batchglm")
-
 
 class Estimator(EstimatorAll, AbstractEstimator, ProcessModel):
     """
@@ -185,7 +183,7 @@ def init_par(
                 init_a_str = init_a.lower()
                 # Chose option if auto was chosen
                 if init_a.lower() == "auto":
-                    init_a = "closed_form"
+                    init_a = "standard"
 
                 if init_a.lower() == "closed_form":
                     groupwise_means, init_a, rmsd_a = closedform_nb_glm_logmu(
@@ -203,8 +201,8 @@ def init_par(
                         if np.any(input_data.size_factors != 1):
                             self._train_loc = True
 
-                    logger.debug("Using closed-form MLE initialization for mean")
-                    logger.debug("Should train mu: %s", self._train_loc)
+                    logging.getLogger("batchglm").debug("Using closed-form MLE initialization for mean")
+                    logging.getLogger("batchglm").debug("Should train mu: %s", self._train_loc)
                 elif init_a.lower() == "standard":
                     if isinstance(input_data.X, SparseXArrayDataArray):
                         overall_means = input_data.X.mean(dim="observations")
@@ -216,14 +214,14 @@ def init_par(
                     init_a[0, :] = np.log(overall_means)
                     self._train_loc = True
 
-                    logger.debug("Using standard initialization for mean")
-                    logger.debug("Should train mu: %s", self._train_loc)
+                    logging.getLogger("batchglm").debug("Using standard initialization for mean")
+                    logging.getLogger("batchglm").debug("Should train mu: %s", self._train_loc)
                 elif init_a.lower() == "all_zero":
                     init_a = np.zeros([input_data.num_loc_params, input_data.num_features])
                     self._train_loc = True
 
-                    logger.debug("Using all_zero initialization for mean")
-                    logger.debug("Should train mu: %s", self._train_loc)
+                    logging.getLogger("batchglm").debug("Using all_zero initialization for mean")
+                    logging.getLogger("batchglm").debug("Should train mu: %s", self._train_loc)
                 else:
                     raise ValueError("init_a string %s not recognized" % init_a)
 
@@ -243,8 +241,8 @@ def init_par(
                     init_b = np.zeros([input_data.num_scale_params, input_data.X.shape[1]])
                     init_b[0, :] = init_b_intercept
 
-                    logger.debug("Using standard-form MME initialization for dispersion")
-                    logger.debug("Should train r: %s", self._train_scale)
+                    logging.getLogger("batchglm").debug("Using standard-form MME initialization for dispersion")
+                    logging.getLogger("batchglm").debug("Should train r: %s", self._train_scale)
                 elif init_b.lower() == "closed_form":
                     dmats_unequal = False
                     if input_data.design_loc.shape[1] == input_data.design_scale.shape[1]:
@@ -269,13 +267,13 @@ def init_par(
                         link_fn=lambda r: np.log(self.np_clip_param(r, "r"))
                     )
 
-                    logger.debug("Using closed-form MME initialization for dispersion")
-                    logger.debug("Should train r: %s", self._train_scale)
+                    logging.getLogger("batchglm").debug("Using closed-form MME initialization for dispersion")
+                    logging.getLogger("batchglm").debug("Should train r: %s", self._train_scale)
                 elif init_b.lower() == "all_zero":
                     init_b = np.zeros([input_data.num_scale_params, input_data.X.shape[1]])
 
-                    logger.debug("Using standard initialization for dispersion")
-                    logger.debug("Should train r: %s", self._train_scale)
+                    logging.getLogger("batchglm").debug("Using standard initialization for dispersion")
+                    logging.getLogger("batchglm").debug("Should train r: %s", self._train_scale)
                 else:
                     raise ValueError("init_b string %s not recognized" % init_b)
         else:
@@ -291,7 +289,7 @@ def init_par(
                     init_loc[my_idx] = init_model.a_var[init_idx]
 
                 init_a = init_loc
-                logger.debug("Using initialization based on input model for mean")
+                logging.getLogger("batchglm").debug("Using initialization based on input model for mean")
 
             # Scale model:
             if isinstance(init_b, str) and (init_b.lower() == "auto" or init_b.lower() == "init_model"):
@@ -305,7 +303,7 @@ def init_par(
                     init_scale[my_idx] = init_model.b_var[init_idx]
 
                 init_b = init_scale
-                logger.debug("Using initialization based on input model for dispersion")
+                logging.getLogger("batchglm").debug("Using initialization based on input model for dispersion")
 
         return init_a, init_b
 

batchglm/train/tf/glm_nb/estimator.py

@@ -1,6 +1,5 @@
 import batchglm.data as data_utils
 
-from batchglm.models.base.input import SparseXArrayDataSet, SparseXArrayDataArray
 from batchglm.models.glm_nb import AbstractEstimator, EstimatorStoreXArray, InputData, Model
 from batchglm.models.base_glm.utils import closedform_glm_mean, closedform_glm_scale
 from batchglm.models.glm_nb.utils import closedform_nb_glm_logmu, closedform_nb_glm_logphi
@@ -17,4 +16,5 @@
 
 import batchglm.utils.random as rand_utils
 from batchglm.utils.linalg import groupwise_solve_lm
+from batchglm.xarray_sparse.base import SparseXArrayDataSet, SparseXArrayDataArray
 from batchglm import pkg_constants

batchglm/train/tf/glm_nb/external.py

@@ -246,15 +246,15 @@ def init_par(
                     # Calculated variance via E(x)^2 or directly depending on whether `mu` was specified.
                     if isinstance(input_data.X, SparseXArrayDataArray):
                         variance = input_data.X.var(input_data.X.dims[0])
-                        variance = np.expand_dims(variance, axis=0)
                     else:
-                        expect_xsq = input_data.X.mean(input_data.X.dims[0])
+                        expect_xsq = np.square(input_data.X).mean(input_data.X.dims[0])
                         mean_model = np.matmul(
                             np.matmul(input_data.design_loc.values, input_data.constraints_loc.values),
                             init_a
                         )
-                        expect_x_sq = np.square(mean_model).mean(input_data.X.dims[0])
-                        variance = expect_xsq - expect_x_sq
+                        expect_x_sq = np.mean(np.square(mean_model), axis=0)  # for xr compatibility input_data.X.dims[0])
+                        variance = (expect_xsq - expect_x_sq).values
+                    variance = np.expand_dims(variance, axis=0)
                     init_b = np.log(np.sqrt(variance))
 
                     self._train_scale = False

batchglm/train/tf/glm_norm/estimator.py

@@ -52,7 +52,8 @@ def _weight_hessian_bb(
         else:
             X_minus_loc = X - loc
 
-        const = - tf.multiply(scalar_two,
+        const = - tf.multiply(
+            scalar_two,
             tf.square(
                 tf.divide(
                     X_minus_loc,

batchglm/train/tf/glm_norm/hessians.py

@@ -1,6 +1,4 @@
 from .test_acc_glm import Test_Accuracy_GLM, _Test_Accuracy_GLM_Estim
-from .test_acc_analytic_glm import Test_AccuracyAnalytic_GLM, _Test_AccuracyAnalytic_GLM_Estim
-from .test_acc_constrained_vglm import Test_AccuracyConstrained_VGLM, _Test_AccuracyConstrained_VGLM_Estim
 from .test_acc_sizefactors_glm import Test_AccuracySizeFactors_GLM, _Test_AccuracySizeFactors_GLM_Estim
 from .test_graph_glm import Test_Graph_GLM, _Test_Graph_GLM_Estim
 from .test_data_types_glm import Test_DataTypes_GLM, _Test_DataTypes_GLM_Estim