Cleaning things up and adding more tests

Author: rvasav26

src/skmatter/decomposition/_pcovc.py

@@ -28,6 +28,7 @@
 # did a search of all classifiers that inherit from MultiOutputMixin - none of them implement
 # decision function, so I don't think we need to inherit
 
+
 class PCovC(LinearClassifierMixin, _BasePCov):
     r"""Principal Covariates Classification (PCovC).
 
@@ -178,16 +179,11 @@ class PCovC(LinearClassifierMixin, _BasePCov):
 
     pxz_ : ndarray of size :math:`({n_{features}, })`, :math:`({n_{features}, n_{classes}})`
         the projector, or weights, from the input space :math:`\mathbf{X}`
-        to the class confidence scores :math:`\mathbf{Z}`. In the multioutput case,
-        has shape , :math:`({n_{features}, n_{classes}*n_{outputs}})`, a flattened form
-        of a 3D tensor.
+        to the class confidence scores :math:`\mathbf{Z}`.
 
-    ptz_ : ndarray of size :math:`({n_{components}, })`, :math:`({n_{components}, n_{classes}})` \
-            or :math:`({n_{components}, n_{classes}*n_{outputs}})`
-        the projector, or weights, from the latent-space projection
-        :math:`\mathbf{T}` to the class confidence scores :math:`\mathbf{Z}`. 
-        In the multioutput case, has shape , :math:`({n_{components}, n_{classes}*n_{outputs}})`, 
-        a flattened form of a 3D tensor.
+    ptz_ : ndarray of size :math:`({n_{components}, })`, :math:`({n_{components}, n_{classes}})`
+        the projector, or weights, from from the latent-space projection
+        :math:`\mathbf{T}` to the class confidence scores :math:`\mathbf{Z}`.
 
     explained_variance_ : numpy.ndarray of shape (n_components,)
         The amount of variance explained by each of the selected components.
@@ -279,7 +275,7 @@ def fit(self, X, Y, W=None):
            `` W = np.hstack([est_.coef_.T for est_ in classifier.estimators_])``.
         """
         X, Y = validate_data(self, X, Y, multi_output=True, y_numeric=False)
-        
+
         check_classification_targets(Y)
         self.classes_ = np.unique(Y)
         self.n_outputs = 1 if Y.ndim == 1 else Y.shape[1]
@@ -305,33 +301,51 @@ def fit(self, X, Y, W=None):
                 "Classifier must be an instance of `"
                 f"{'`, `'.join(c.__name__ for c in compatible_classifiers)}`"
                 ", or `precomputed`"
-          )
+            )
 
-        # if self.n_outputs == 1:
-        #     classifier = LogisticRegression()
-        # else:
-        #     classifier = MultiOutputClassifier(estimator=LogisticRegression())
+        if self.n_outputs == 1 and isinstance(self.classifier, MultiOutputClassifier):
+            raise ValueError(
+                "Classifier cannot be an instance of `MultiOutputClassifier` when Y is 1D"
+            )
+
+        if (
+            self.n_outputs != 1
+            and self.classifier not in ["precomputed", None]
+            and not (
+                isinstance(self.classifier, MultiOutputClassifier)
+                or self.classifier == "precomputed"
+            )
+        ):
+            raise ValueError(
+                "Classifier must be an instance of `MultiOutputClassifier` when Y is 2D"
+            )
 
-        # if self.classifier == "precomputed":
-            
+        if self.n_outputs == 1:
+            if self.classifier != "precomputed":
+                classifier = self.classifier or LogisticRegression()
+                self.z_classifier_ = check_cl_fit(classifier, X, Y)
+                W = self.z_classifier_.coef_.T
 
-        if self.classifier != "precomputed":
-            if self.classifier is None:
-                classifier = LogisticRegression()
             else:
-                classifier = self.classifier
-
-            self.z_classifier_ = check_cl_fit(classifier, X, Y)
-            W = self.z_classifier_.coef_.T
+                # to be used later on as the classifier fit between T and Y
+                classifier = LogisticRegression()
+                if W is None:
+                    W = clone(classifier).fit(X, Y).coef_.T
 
         else:
-            # If precomputed, use default classifier to predict Y from T
-            classifier = LogisticRegression()
-            if W is None:
-                W = LogisticRegression().fit(X, Y).coef_.T
+            if self.classifier != "precomputed":
+                classifier = self.classifier or MultiOutputClassifier(
+                    estimator=LogisticRegression()
+                )
+                self.z_classifier_ = check_cl_fit(classifier, X, Y)
+                W = np.hstack([est_.coef_.T for est_ in self.z_classifier_.estimators_])
 
-        print(f"X: {X.shape}")
-        print(f"W: {W.shape}")
+            else:
+                # to be used later on as the classifier fit between T and Y
+                classifier = MultiOutputClassifier(estimator=LogisticRegression())
+                if W is None:
+                    _ = clone(classifier).fit(X, Y)
+                    W = np.hstack([_.coef_.T for _ in _.estimators_])
 
         Z = X @ W
 
@@ -344,7 +358,11 @@ def fit(self, X, Y, W=None):
         # classifier and steal weights to get pxz and ptz
         self.classifier_ = clone(classifier).fit(X @ self.pxt_, Y)
 
-        if isinstance(self.classifier_, MultiOutputClassifier):
+        if self.n_outputs == 1:
+            self.ptz_ = self.classifier_.coef_.T
+            # print(self.ptz_.shape)
+            self.pxz_ = self.pxt_ @ self.ptz_
+        else:
             self.ptz_ = np.hstack(
                 [est_.coef_.T for est_ in self.classifier_.estimators_]
             )
@@ -353,12 +371,7 @@ def fit(self, X, Y, W=None):
             self.pxz_ = self.pxt_ @ self.ptz_
             # print(f"pxz {self.pxz_.shape}")
 
-        else:
-            self.ptz_ = self.classifier_.coef_.T
-            # print(self.ptz_.shape)
-            self.pxz_ = self.pxt_ @ self.ptz_
-
-        print(self.ptz_.shape)
+        # print(self.ptz_.shape)
         if len(Y.shape) == 1 and type_of_target(Y) == "binary":
             self.pxz_ = self.pxz_.reshape(
                 X.shape[1],
@@ -460,7 +473,7 @@ def decision_function(self, X=None, T=None):
                 n_outputs such arrays if n_outputs > 1
             Confidence scores. For binary classification, has shape `(n_samples,)`,
             for multiclass classification, has shape `(n_samples, n_classes)`. If n_outputs > 1,
-            the list returned can contain arrays with differing shapes depending on the
+            the list can contain arrays with differing shapes depending on the
             number of classes in each output of Y.
         """
         check_is_fitted(self, attributes=["pxz_", "ptz_"])
@@ -471,25 +484,24 @@ def decision_function(self, X=None, T=None):
         if X is not None:
             X = validate_data(self, X, reset=False)
 
-            # this is similar to how MultiOutputClassifier handles predict_proba() if n_outputs > 1
-            if isinstance(self.classifier_, MultiOutputClassifier):
+            if self.n_outputs == 1:
+                # Or self.classifier_.decision_function(X @ self.pxt_)
+                return X @ self.pxz_ + self.classifier_.intercept_
+            else:
                 return [
                     est_.decision_function(X @ self.pxt_)
                     for est_ in self.classifier_.estimators_
                 ]
-
-            # Or self.classifier_.decision_function(X @ self.pxt_)
-            return X @ self.pxz_ + self.classifier_.intercept_
         else:
             T = check_array(T)
 
-            if isinstance(self.classifier_, MultiOutputClassifier):
+            if self.n_outputs == 1:
+                return T @ self.ptz_ + self.classifier_.intercept_
+            else:
                 return [
                     est_.decision_function(T) for est_ in self.classifier_.estimators_
                 ]
 
-            return T @ self.ptz_ + self.classifier_.intercept_
-
     def predict(self, X=None, T=None):
         """Predicts the property labels using classification on T."""
         check_is_fitted(self, attributes=["pxz_", "ptz_"])

tests/test_pcovc.py

@@ -3,7 +3,8 @@
 
 import numpy as np
 from sklearn import exceptions
-from sklearn.datasets import load_breast_cancer as get_dataset
+from sklearn.calibration import LinearSVC
+from sklearn.datasets import load_iris as get_dataset
 from sklearn.decomposition import PCA
 from sklearn.linear_model import LogisticRegression, RidgeClassifier
 from sklearn.svm import LinearSVC
@@ -98,7 +99,7 @@ def test_simple_prediction(self):
                 Yp = pcovc.predict(self.X)
 
                 self.assertLessEqual(
-                    np.linalg.norm(Yp - Yhat) ** 2.0 / np.linalg.norm(Yp) ** 2.0,
+                    np.linalg.norm(Yp - Yhat) ** 2.0 / np.linalg.norm(Yhat) ** 2.0,
                     self.error_tol,
                 )
 
@@ -580,15 +581,56 @@ def test_incompatible_coef_shape(self):
 
 class PCovCMultiOutputTest(PCovCBaseTest):
 
-    def test_projector_shapes(self):
-        pass
+    def test_prefit_multioutput(self):
+        """Check that PCovC works if a prefit classifier is passed when `n_ouputs > 1`."""
+        classifier = MultiOutputClassifier(estimator=LogisticRegression())
+        Y_double = np.column_stack((self.Y, self.Y))
 
-    def test_decision_function(self):
+        classifier.fit(self.X, Y_double)
+        pcovc = self.model(mixing=0.25, classifier=classifier)
+        pcovc.fit(self.X, Y_double)
+
+        W_classifier = np.hstack([est_.coef_.T for est_ in classifier.estimators_])
+        Z_classifier = self.X @ W_classifier
+
+        W_pcovc = np.hstack([est_.coef_.T for est_ in pcovc.z_classifier_.estimators_])
+        Z_pcovc = self.X @ W_pcovc
+
+        self.assertTrue(np.allclose(Z_classifier, Z_pcovc))
+        self.assertTrue(np.allclose(W_classifier, W_pcovc))
+
+    def test_precomputed_multioutput(self):
+        """Check that PCovC works if classifier=`precomputed` and `n_ouputs > 1`."""
+        classifier = MultiOutputClassifier(estimator=LogisticRegression())
+        Y_double = np.column_stack((self.Y, self.Y))
+
+        classifier.fit(self.X, Y_double)
+        W = np.hstack([est_.coef_.T for est_ in classifier.estimators_])
+        pcovc1 = self.model(mixing=0.5, classifier="precomputed", n_components=1)
+        pcovc1.fit(self.X, Y_double, W)
+        t1 = pcovc1.transform(self.X)
+
+        pcovc2 = self.model(mixing=0.5, classifier=classifier, n_components=1)
+        pcovc2.fit(self.X, Y_double)
+        t2 = pcovc2.transform(self.X)
+
+        self.assertTrue(np.linalg.norm(t1 - t2) < self.error_tol)
+
+        # Now check for match when W is not passed:
+        pcovc3 = self.model(mixing=0.5, classifier="precomputed", n_components=1)
+        pcovc3.fit(self.X, Y_double)
+        t3 = pcovc3.transform(self.X)
+
+        self.assertTrue(np.linalg.norm(t3 - t2) < self.error_tol)
+        self.assertTrue(np.linalg.norm(t3 - t1) < self.error_tol)
+
+    def test_Z_shape_multioutput(self):
+        """Check that PCovC returns the evidence Z in the desired form when `n_ouputs > 1`."""
         pcovc = PCovC(
             classifier=MultiOutputClassifier(LogisticRegression()), n_components=2
         )
 
-        Y_double = np.column_stack((self.Y, self.Y[::-1]))
+        Y_double = np.column_stack((self.Y, self.Y))
         pcovc.fit(self.X, Y_double)
 
         Z = pcovc.decision_function(self.X)
@@ -602,6 +644,20 @@ def test_decision_function(self):
                 self.assertEqual(self.X.shape[0], z_slice.shape[0])
                 self.assertEqual(est.coef_.shape[0], z_slice.shape[1])
 
+    def test_decision_function_multioutput(self):
+        """Check that PCovC's decision_function works in edge cases when `n_ouputs > 1`."""
+        pcovc = self.model(classifier=MultiOutputClassifier(estimator=LinearSVC()))
+        pcovc.fit(self.X, np.column_stack((self.Y, self.Y)))
+        with self.assertRaises(ValueError) as cm:
+            _ = pcovc.decision_function()
+        self.assertEqual(
+            str(cm.exception),
+            "Either X or T must be supplied.",
+        )
+
+        T = pcovc.transform(self.X)
+        _ = pcovc.decision_function(T=T)
+
 
 if __name__ == "__main__":
     unittest.main(verbosity=2)