Merge pull request #86 from theislab/dev

Dev

Author: davidsebfischer

batchglm/data.py

@@ -21,7 +21,7 @@ def design_matrix(
         as_categorical: Union[bool, list] = True,
         dmat: Union[pd.DataFrame, None] = None,
         return_type: str = "patsy",
-) -> Union[patsy.design_info.DesignMatrix, pd.DataFrame]:
+) -> Tuple[Union[patsy.design_info.DesignMatrix, pd.DataFrame], List[str]]:
     """
     Create a design matrix from some sample description.
 
@@ -62,6 +62,7 @@ def design_matrix(
                     sample_description[col] = sample_description[col].astype("category")
 
         dmat = patsy.dmatrix(formula, sample_description)
+        coef_names = dmat.design_info.column_names
 
         if return_type == "dataframe":
             df = pd.DataFrame(dmat, columns=dmat.design_info.column_names)
@@ -70,12 +71,12 @@ def design_matrix(
 
             return df
         elif return_type == "patsy":
-            return dmat
+            return dmat, coef_names
         else:
             raise ValueError("return type %s not recognized" % return_type)
     else:
         if return_type == "dataframe":
-            return dmat
+            return dmat, dmat.columns
         elif return_type == "patsy":
             raise ValueError("return type 'patsy' not supported for input (dmat is not None)")
         else:
@@ -105,7 +106,7 @@ def preview_coef_names(
         sample_description: pd.DataFrame,
         formula: str,
         as_categorical: Union[bool, list] = True
-) -> np.ndarray:
+) -> List[str]:
     """
     Return coefficient names of model.
 
@@ -125,21 +126,22 @@ def preview_coef_names(
         Set to false, if columns should not be changed.
     :return: A list of coefficient names.
     """
-    return view_coef_names(dmat=design_matrix(
+    _, coef_names = design_matrix(
         sample_description=sample_description,
         formula=formula,
         as_categorical=as_categorical,
         return_type="patsy"
-    ))
+    )
+    return coef_names
 
 
 def constraint_system_from_star(
-        dmat: Union[None, np.ndarray] = None,
+        dmat: Union[None, patsy.design_info.DesignMatrix, pd.DataFrame] = None,
         sample_description: Union[None, pd.DataFrame] = None,
         formula: Union[None, str] = None,
         as_categorical: Union[bool, list] = True,
         constraints: Union[None, List[str], Tuple[str], dict, np.ndarray] = None,
-        return_type: str = "patsy",
+        return_type: str = "patsy"
 ) -> Tuple:
     """
     Wrap different constraint matrix building formats with building of design matrix.
@@ -202,7 +204,7 @@ def constraint_system_from_star(
         raise ValueError("supply either sample_description or dmat")
 
     if dmat is None and not isinstance(constraints, dict):
-       dmat = design_matrix(
+       dmat, coef_names = design_matrix(
             sample_description=sample_description,
             formula=formula,
             as_categorical=as_categorical,
@@ -213,39 +215,69 @@ def constraint_system_from_star(
         raise ValueError("dmat was supplied even though constraints were given as dict")
 
     if isinstance(constraints, dict):
-        dmat, cmat = constraint_matrix_from_dict(
+        dmat, coef_names, cmat, term_names = constraint_matrix_from_dict(
             sample_description=sample_description,
             formula=formula,
             as_categorical=as_categorical,
             constraints=constraints,
-            return_type="dataframe"
+            return_type="patsy"
         )
     elif isinstance(constraints, tuple) or isinstance(constraints, list):
-        cmat = constraint_matrix_from_string(
+        cmat, coef_names = constraint_matrix_from_string(
             dmat=dmat,
+            coef_names=dmat.design_info.column_names,
             constraints=constraints
         )
+        term_names = None  # not supported yet.
     elif isinstance(constraints, np.ndarray):
         cmat = constraints
+        term_names = None
+        if isinstance(dmat, pd.DataFrame):
+            coef_names = dmat.columns
+            dmat = dmat.values
     elif constraints is None:
         cmat = None
+        term_names = None
+        if isinstance(dmat, pd.DataFrame):
+            coef_names = dmat.columns
+            dmat = dmat.values
     else:
         raise ValueError("constraint format %s not recognized" % type(constraints))
 
-    return dmat, cmat
+    # Test full design matrix for being full rank before returning:
+    if cmat is None:
+        if np.linalg.matrix_rank(dmat) != dmat.shape[1]:
+            raise ValueError(
+                "constrained design matrix is not full rank: %i %i" %
+                (np.linalg.matrix_rank(dmat), dmat.shape[1])
+            )
+    else:
+        if np.linalg.matrix_rank(np.matmul(dmat, cmat)) != cmat.shape[1]:
+            raise ValueError(
+                "constrained design matrix is not full rank: %i %i" %
+                (np.linalg.matrix_rank(np.matmul(dmat, cmat)), cmat.shape[1])
+            )
+
+    return dmat, coef_names, cmat, term_names
 
 
 def constraint_matrix_from_dict(
         sample_description: pd.DataFrame,
         formula: str,
         as_categorical: Union[bool, list] = True,
         constraints: dict = {},
-        return_type: str = "dataframe"
+        return_type: str = "patsy"
 ) -> Tuple:
     """
     Create a design matrix from some sample description and a constraint matrix
     based on factor encoding of constrained parameter sets.
 
+    Note that we build a dataframe instead of a pasty.DesignMatrix here if constraints are used.
+    This is done because we were not able to build a patsy.DesignMatrix of the constrained form
+    required in this context. In those cases in which the return type cannot be patsy, we encourage the
+    use of the returned term_names to perform term-wise slicing which is not supported by other
+    design matrix return types.
+
     :param sample_description: pandas.DataFrame of length "num_observations" containing explanatory variables as columns
     :param formula: model formula as string, describing the relations of the explanatory variables.
 
@@ -270,7 +302,9 @@ def constraint_matrix_from_dict(
 
         Can only group by non-constrained effects right now, use constraint_matrix_from_string
         for other cases.
-    :return: a model design matrix
+    :return:
+        - model design matrix
+        - term_names to allow slicing by factor if return type cannot be patsy.DesignMatrix
     """
     assert len(constraints) > 0, "supply constraints"
     sample_description: pd.DataFrame = sample_description.copy()
@@ -287,10 +321,11 @@ def constraint_matrix_from_dict(
     # absorption of the first level of each factor for each constrained factor onto the
     # core matrix.
     formula_unconstrained = formula.split("+")
-    formula_unconstrained = [x for x in formula_unconstrained if x not in constraints.keys()]
+    formula_unconstrained = [x for x in formula_unconstrained if x.strip(" ") not in constraints.keys()]
     formula_unconstrained = "+".join(formula_unconstrained)
     dmat = patsy.dmatrix(formula_unconstrained, sample_description)
     coef_names = dmat.design_info.column_names
+    term_names = dmat.design_info.term_names
 
     constraints_ls = string_constraints_from_dict(
         sample_description=sample_description,
@@ -301,6 +336,7 @@ def constraint_matrix_from_dict(
         dmat_constrained_temp = patsy.highlevel.dmatrix("0+" + x, sample_description)
         dmat = np.hstack([dmat, dmat_constrained_temp])
         coef_names.extend(dmat_constrained_temp.design_info.column_names)
+        term_names.extend(dmat_constrained_temp.design_info.term_names)
 
     # Build constraint matrix.
     constraints_ar = constraint_matrix_from_string(
@@ -312,8 +348,7 @@ def constraint_matrix_from_dict(
     # Format return type
     if return_type == "dataframe":
         dmat = pd.DataFrame(dmat, columns=coef_names)
-
-    return dmat, constraints_ar
+    return dmat, coef_names, constraints_ar, term_names
 
 
 def string_constraints_from_dict(
@@ -388,6 +423,10 @@ def constraint_matrix_from_string(
 
     di = patsy.DesignInfo(coef_names)
     constraint_ls = [di.linear_constraint(x).coefs[0] for x in constraints]
+    # Check that constraints are sensible:
+    for constraint_i in constraint_ls:
+        if np.sum(constraint_i != 0) == 1:
+            raise ValueError("a zero-equality constraint only involved one parameter: remove this parameter")
     idx_constr = np.asarray([np.where(x == 1)[0][0] for x in constraint_ls])
     idx_depending = [np.where(x == 1)[0][1:] for x in constraint_ls]
     idx_unconstr = np.asarray(list(
@@ -407,8 +446,10 @@ def constraint_matrix_from_string(
             constraint_mat[i, idx_unconstr_i] = 1
 
     # Test unconstrained subset design matrix for being full rank before returning constraints:
-    dmat_var =dmat[:, idx_unconstr]
-    if np.linalg.matrix_rank(dmat_var) != np.linalg.matrix_rank(dmat_var.T):
-        logging.getLogger("batchglm").error("constrained design matrix is not full rank")
+    if np.linalg.matrix_rank(dmat[:, idx_unconstr]) != np.linalg.matrix_rank(dmat[:, idx_unconstr].T):
+        raise ValueError(
+            "unconstrained sub-design matrix is not full rank" %
+            np.linalg.matrix_rank(dmat[:, idx_unconstr]), np.linalg.matrix_rank(dmat[:, idx_unconstr].T)
+        )
 
     return constraint_mat

batchglm/models/base_glm/input.py

@@ -88,7 +88,7 @@ def __init__(
             design_matrix=design_loc,
             param_names=design_loc_names
         )
-        design_scale, design_scale_names =  parse_design(
+        design_scale, design_scale_names = parse_design(
             design_matrix=design_scale,
             param_names=design_scale_names
         )

batchglm/models/base_glm/simulator.py

@@ -3,7 +3,7 @@
 import numpy as np
 import pandas
 import patsy
-from typing import Union
+from typing import Union, Tuple
 
 from .model import _ModelGLM
 from .external import _SimulatorBase
@@ -14,7 +14,7 @@ def generate_sample_description(
         num_conditions: int = 2,
         num_batches: int = 4,
         shuffle_assignments=False
-) -> np.ndarray:
+) -> Tuple[patsy.DesignMatrix, pandas.DataFrame]:
     """ Build a sample description.
 
     :param num_observations: Number of observations to simulate.
@@ -45,7 +45,7 @@ def generate_sample_description(
             observations=np.random.permutation(sample_description.observations.values)
         )
 
-    return np.asarray(patsy.dmatrix("~1+condition+batch", sample_description)), sample_description
+    return patsy.dmatrix("~1+condition+batch", sample_description), sample_description
 
 
 class _SimulatorGLM(_SimulatorBase, metaclass=abc.ABCMeta):
@@ -73,6 +73,7 @@ def __init__(
         self.sample_description = None
         self.sim_a_var = None
         self.sim_b_var = None
+        self._size_factors = None
 
     def generate_sample_description(
             self,
@@ -139,6 +140,10 @@ def _generate_params(
             rand_fn_scale((self.sim_design_scale.shape[1], self.nfeatures))
         ], axis=0)
 
+    @property
+    def size_factors(self):
+        return self._size_factors
+
     @property
     def a_var(self):
         return self.sim_a_var

batchglm/models/base_glm/utils.py

@@ -37,8 +37,8 @@ def parse_design(
         assert False
 
     if param_names is not None:
-        if params is not None:
-            assert len(param_names) == len(params)
+        if params is None:
+            assert len(param_names) == dmat.shape[1]
             params = param_names
 
     return dmat, params
@@ -62,8 +62,11 @@ def parse_constraints(
         constraints = np.identity(n=dmat.shape[1])
         constraint_params = dmat_par_names
     else:
-        # Cannot use given parameter names if constraint matrix is not identity: Make up new ones.
-        constraint_params = ["var_"+str(x) for x in range(constraints.shape[1])]
+        # Cannot use all parameter names if constraint matrix is not identity: Make up new ones.
+        # Use variable names that can be mapped (unconstrained).
+        constraint_params = ["var_"+str(i) if np.sum(constraints[:, i] != 0) > 1
+                             else dmat_par_names[np.where(constraints[:, i] != 0)[0][0]]
+                             for i in range(constraints.shape[1])]
         assert constraints.shape[0] == dmat.shape[1], "constraint dimension mismatch"
 
     if constraint_par_names is not None:

batchglm/models/glm_nb/simulator.py

@@ -15,6 +15,10 @@ def __init__(
             num_observations=1000,
             num_features=100
     ):
+        Model.__init__(
+            self=self,
+            input_data=None
+        )
         _SimulatorGLM.__init__(
             self=self,
             model=None,
@@ -53,3 +57,4 @@ def generate_data(self):
             design_loc_names=None,
             design_scale_names=None
         )
+

batchglm/unit_test/test_acc_constrained_vglm_all.py

@@ -24,8 +24,11 @@ def simulate(self):
         constraints[1, 1] = 1
         constraints[2, 2] = 1
         constraints[3, 2] = -1
+        new_coef_names = ['Intercept', 'condition[T.1]', 'batch[1]', 'batch[2]']
         self.sim1.input_data.design_loc = dmat
         self.sim1.input_data.design_scale = dmat
+        self.sim1.input_data._design_loc_names = new_coef_names
+        self.sim1.input_data._design_scale_names = new_coef_names
         self.sim1.input_data.constraints_loc = constraints
         self.sim1.input_data.constraints_scale = constraints
 

batchglm/unit_test/test_acc_glm_all.py

@@ -50,6 +50,8 @@ def __init__(
                 data=scipy.sparse.csr_matrix(simulator.input_data.x),
                 design_loc=simulator.input_data.design_loc,
                 design_scale=simulator.input_data.design_scale,
+                design_loc_names=simulator.input_data.design_loc_names,
+                design_scale_names=simulator.input_data.design_scale_names,
                 constraints_loc=simulator.input_data.constraints_loc,
                 constraints_scale=simulator.input_data.constraints_scale,
                 size_factors=simulator.input_data.size_factors
@@ -59,6 +61,8 @@ def __init__(
                 data=simulator.input_data.x,
                 design_loc=simulator.input_data.design_loc,
                 design_scale=simulator.input_data.design_scale,
+                design_loc_names=simulator.input_data.design_loc_names,
+                design_scale_names=simulator.input_data.design_scale_names,
                 constraints_loc=simulator.input_data.constraints_loc,
                 constraints_scale=simulator.input_data.constraints_scale,
                 size_factors=simulator.input_data.size_factors