Affine preprocessing is supported by mathematical programs

ceml/optim/cvx.py

@@ -16,30 +16,50 @@ def solve(self):
         raise NotImplementedError()
 
 
-class ConvexQuadraticProgram(ABC):
-    """Base class for a convex quadratic program - for computing counterfactuals.
-
-    Attributes
-    ----------
-    epsilon : `float`
-        "Small" non-negative number for relaxing strict inequalities.
+class SupportAffinePreprocessing():
+    """Base class for a mathematical programs that support an affine preprocessing.
     """
     def __init__(self, **kwds):
-        self.epsilon = 1e-2
         self.A = None
         self.b = None
 
         super().__init__(**kwds)
-    
+
     def set_affine_preprocessing(self, A, b):
         self.A = A
         self.b = b
 
-    def _apply_affine_preprocessing(self, var_x):
+    def get_affine_preprocessing(self):
+        return {"A": self.A, "b": self.b}
+
+    def is_affine_preprocessing_set(self):
+        return self.A is None or self.b is None
+
+    def _apply_affine_preprocessing_to_var(self, var_x):
         if self.A is not None and self.b is not None:
             return self.A @ var_x + self.b
         else:
             return var_x
+
+    def _apply_affine_preprocessing_to_const(self, x):
+        if self.A is not None and self.b is not None:
+            return np.dot(self.A, x) + self.b
+        else:
+            return x
+
+
+class ConvexQuadraticProgram(ABC, SupportAffinePreprocessing):
+    """Base class for a convex quadratic program - for computing counterfactuals.
+
+    Attributes
+    ----------
+    epsilon : `float`
+        "Small" non-negative number for relaxing strict inequalities.
+    """
+    def __init__(self, **kwds):
+        self.epsilon = 1e-2
+
+        super().__init__(**kwds)
 
     @abstractmethod
     def _build_constraints(self, var_x, y):
@@ -244,7 +264,7 @@ def build_solve_opt(self, x_orig, y, features_whitelist=None):
         return x.value.reshape(dim)
 
 
-class DCQP():
+class DCQP(SupportAffinePreprocessing):
     """Class for a difference-of-convex-quadratic program (DCQP) - for computing counterfactuals.
 
     .. math:: \\underset{\\vec{x} \\in \\mathbb{R}^d}{\\min} \\vec{x}^\\top Q_0 \\vec{x} + \\vec{q}^\\top \\vec{x} + c - \\vec{x}^\\top Q_1 \\vec{x} \\quad \\text{s.t. } \\vec{x}^\\top A0_i \\vec{x} + \\vec{x}^\\top \\vec{b_i} + r_i - \\vec{x}^\\top A1_i \\vec{x} \\leq 0 \\; \\forall\\,i
@@ -326,10 +346,11 @@ def solve(self, x0, tao=1.2, tao_max=100, mu=1.5):
 
             The default is 1.5
         """
+        self.pccp.set_affine_preprocessing(**self.get_affine_preprocessing())
         return self.pccp.compute_counterfactual(self.x_orig, self.y_target, x0, tao=1.2, tao_max=100, mu=1.5)
 
 
-class PenaltyConvexConcaveProcedure():
+class PenaltyConvexConcaveProcedure(SupportAffinePreprocessing):
     """Implementation of the penalty convex-concave procedure for approximately solving a DCQP.
     """
     def __init__(self, model, Q0, Q1, q, c, A0_i, A1_i, b_i, r_i, features_whitelist=None, mad=None, **kwds):      
@@ -362,9 +383,11 @@ def solve_aux(self, xcf, tao, x_orig):
             self.dim = x_orig.shape[0]
 
             # Variables
-            x = cp.Variable(self.dim)
+            var_x = cp.Variable(self.dim)
             s = cp.Variable(len(self.A0s))
 
+            var_x_prime = self._apply_affine_preprocessing_to_var(var_x)
+
             # Constants
             s_z = np.zeros(len(self.A0s))
             s_c = np.ones(len(self.A0s))
@@ -373,9 +396,9 @@ def solve_aux(self, xcf, tao, x_orig):
             # Build constraints
             constraints = []
             for i in range(len(self.A0s)):
-                A = cp.quad_form(x, self.A0s[i])
-                q = x.T @ self.bs[i]
-                c = self.rs[i] + np.dot(xcf, np.dot(xcf, self.A1s[i])) - 2. * x.T @ np.dot(xcf, self.A1s[i]) - s[i]
+                A = cp.quad_form(var_x_prime, self.A0s[i])
+                q = var_x_prime.T @ self.bs[i]
+                c = self.rs[i] + np.dot(xcf, np.dot(xcf, self.A1s[i])) - 2. * var_x.T @ np.dot(xcf, self.A1s[i]) - s[i]
 
                 constraints.append(A + q + c + self.epsilon <= 0)
 
@@ -395,25 +418,25 @@ def solve_aux(self, xcf, tao, x_orig):
                     A = np.array(A)
                     a = np.array(a)
 
-                    constraints += [A @ x == a]
+                    constraints += [A @ var_x == a]
 
             # Build the final program
             f = None
             if self.mad is not None:    # TODO: Right now, mad != 1 is not supported.
-                f = cp.Minimize(cp.norm(x - x_orig, 1) + s.T @ (tao*s_c))
+                f = cp.Minimize(cp.norm(var_x - x_orig, 1) + s.T @ (tao*s_c))
             else:
-                f = cp.Minimize(cp.quad_form(x, self.Q0) + self.q.T @ x + self.c + np.dot(xcf, np.dot(xcf, self.Q1)) - 2. * x.T @ np.dot(xcf, self.Q1) + s.T @ (tao*s_c))
+                f = cp.Minimize(cp.quad_form(var_x, self.Q0) + self.q.T @ var_x + self.c + np.dot(xcf, np.dot(xcf, self.Q1)) - 2. * var_x.T @ np.dot(xcf, self.Q1) + s.T @ (tao*s_c))
             constraints += [s >= s_z]
 
             prob = cp.Problem(f, constraints)
 
             # Solve it!
             self._solve(prob)
 
-            if x.value is None:
+            if var_x.value is None:
                 raise Exception("No solution found!")
             else:
-                return x.value
+                return var_x.value
         except Exception as ex:
             logging.debug(str(ex))
 

ceml/sklearn/lda.py

@@ -130,7 +130,7 @@ def _build_constraints(self, var_x, y):
         constraints = []
 
         # If set, a apply an affine preprocessing to x
-        var_x_ = self._apply_affine_preprocessing(var_x)
+        var_x_ = self._apply_affine_preprocessing_to_var(var_x)
 
         # Build constraints
         i = y

ceml/sklearn/linearregression.py

@@ -114,7 +114,7 @@ def _build_constraints(self, var_x, y):
         constraints = []
 
         # If set, a apply an affine preprocessing to x
-        var_x_ = self._apply_affine_preprocessing(var_x)
+        var_x_ = self._apply_affine_preprocessing_to_var(var_x)
 
         # Build box constraints
         constraints.append(self.mymodel.w @ var_x_ + self.mymodel.b - y <= self.epsilon)

ceml/sklearn/lvq.py

@@ -163,6 +163,8 @@ def __init__(self, mymodel, x_orig, y_target, indices_other_prototypes, features
         super().__init__(**kwds)
 
     def _build_constraints(self, var_x, y):
+        var_x_prime = self._apply_affine_preprocessing_to_var(var_x)
+
         Omega = self.mymodel._get_omega()
         p_i = self.mymodel.prototypes[self.target_prototype]
 
@@ -171,7 +173,7 @@ def _build_constraints(self, var_x, y):
             p_j = self.other_prototypes[k]
             qj = np.dot(Omega, p_j - p_i)
             bj = -.5 * (np.dot(p_i, np.dot(Omega, p_i)) - np.dot(p_j, np.dot(Omega, p_j)))
-            results.append(qj.T @ var_x + bj + self.epsilon <= 0)
+            results.append(qj.T @ var_x_prime + bj + self.epsilon <= 0)
 
         return results
 
@@ -231,6 +233,7 @@ def _compute_counterfactual_via_convex_quadratic_programming(self, x_orig, y_tar
 
         # Compute a counterfactual for each prototype
         solver = CQPHelper(mymodel=self.mymodel, x_orig=x_orig, y_target=y_target, indices_other_prototypes=other_prototypes, features_whitelist=features_whitelist, regularization=regularization)
+        solver.set_affine_preprocessing(**self.get_affine_preprocessing())
         for i in range(len(target_prototypes)):
             try:
                 xcf_ = solver.solve(target_prototype=i, features_whitelist=features_whitelist)
@@ -251,6 +254,8 @@ def _compute_counterfactual_via_convex_quadratic_programming(self, x_orig, y_tar
         return xcf
 
     def _build_solve_dcqp(self, x_orig, y_target, target_prototype_id, other_prototypes, features_whitelist, regularization):
+        x_orig_prime = self._apply_affine_preprocessing_to_const(x_orig)
+
         p_i = self.mymodel.prototypes[target_prototype_id]
         o_i = self.mymodel.dist_mats[target_prototype_id] if not self.mymodel.classwise else self.mymodel.omegas[y_target]
         ri = .5 * np.dot(p_i, np.dot(o_i, p_i))
@@ -263,7 +268,7 @@ def _build_solve_dcqp(self, x_orig, y_target, target_prototype_id, other_prototy
         if regularization == "l2":
             Q0 = np.eye(self.mymodel.dim)
             Q1 = np.zeros((self.mymodel.dim, self.mymodel.dim))
-            q = -x_orig
+            q = -x_orig_prime
             c = 0.0
 
         A0_i = []

ceml/sklearn/naivebayes.py

@@ -135,9 +135,11 @@ def _build_constraints(self, var_X, var_x, y):
         return [cp.trace(A @ var_X) + b @ var_x + c + self.epsilon <= 0]
 
     def _build_solve_dcqp(self, x_orig, y_target, regularization, features_whitelist):
+        x_orig_prime = self._apply_affine_preprocessing_to_const(x_orig)
+
         Q0 = np.eye(self.mymodel.dim)   # TODO: Can be ignored if regularization != l2
         Q1 = np.zeros((self.mymodel.dim, self.mymodel.dim))
-        q = -2. * x_orig
+        q = -2. * x_orig_prime
         c = 0.0
 
         A0_i = []
@@ -158,7 +160,7 @@ def _build_solve_dcqp(self, x_orig, y_target, regularization, features_whitelist
 
     def solve(self, x_orig, y_target, regularization, features_whitelist, return_as_dict):
         xcf = None
-        if self.mymodel.is_binary:
+        if self.mymodel.is_binary and not self.is_affine_preprocessing_set():
             xcf = self.build_solve_opt(x_orig, y_target)
         else:
             xcf = self._build_solve_dcqp(x_orig, y_target, regularization, features_whitelist)

ceml/sklearn/pipeline.py

@@ -4,17 +4,18 @@
 import sklearn_lvq
 
 from .softmaxregression import SoftmaxRegression, SoftmaxCounterfactual
-from .naivebayes import GaussianNB
+from .naivebayes import GaussianNB, GaussianNbCounterfactual
 from .linearregression import LinearRegression, LinearRegressionCounterfactual
 from .knn import KNN
-from .lvq import LVQ
+from .lvq import LVQ, LvqCounterfactual
 from .lda import Lda, LdaCounterfactual
-from .qda import Qda
+from .qda import Qda, QdaCounterfactual
 from ..model import ModelWithLoss
 from ..backend.jax.preprocessing import StandardScaler, PCA, PolynomialFeatures, Normalizer, MinMaxScaler, AffinePreprocessing, concatenate_affine_mappings
 from ..backend.jax.costfunctions import CostFunctionDifferentiable, RegularizedCost
 from ..costfunctions import RegularizedCost as RegularizedCostNonDifferentiable
 from ..optim.cvx import ConvexQuadraticProgram
+from ..optim.cvx import SupportAffinePreprocessing
 from .utils import build_regularization_loss
 from .counterfactual import SklearnCounterfactual
 
@@ -142,14 +143,14 @@ def wrap_model(self, model, return_sklearn_counterfactual=False):
             if return_sklearn_counterfactual is False:
                 return GaussianNB(model)
             else:
-                raise NotImplementedError()
+                return GaussianNbCounterfactual(model)
         elif isinstance(model, sklearn.discriminant_analysis.LinearDiscriminantAnalysis):
             return Lda(model) if return_sklearn_counterfactual is False else LdaCounterfactual(model)
         elif isinstance(model, sklearn.discriminant_analysis.QuadraticDiscriminantAnalysis):
             if return_sklearn_counterfactual is False:
                 return Qda(model)
             else:
-                raise NotImplementedError()
+                return QdaCounterfactual(model)
         elif isinstance(model, sklearn.tree.DecisionTreeClassifier) or isinstance(model, sklearn.tree.DecisionTreeRegressor):
             raise NotImplementedError()
         elif isinstance(model, sklearn.ensemble.RandomForestClassifier) or isinstance(model, sklearn.ensemble.RandomForestRegressor):
@@ -163,7 +164,7 @@ def wrap_model(self, model, return_sklearn_counterfactual=False):
             if return_sklearn_counterfactual is False:
                 return LVQ(model)
             else:
-                raise NotImplementedError()
+                return LvqCounterfactual(model)
         else:
             raise ValueError(f"{type(model)} is not supported")
 
@@ -323,7 +324,7 @@ def compute_counterfactual(self, x, y_target, features_whitelist=None, regulariz
         """
         if optimizer == "auto":
             # Check if we can use a mathematical program
-            if isinstance(self.last_model_sklearn_counterfactual, ConvexQuadraticProgram): # TODO: Support SDPs and DCPs
+            if isinstance(self.last_model_sklearn_counterfactual, SupportAffinePreprocessing):
                 if all([isinstance(m, AffinePreprocessing) for m in self.mymodel.models[:-1]]):
                     optimizer = "mp"
             else:   # Use Downhill-Simplex method otherwise
@@ -334,8 +335,8 @@ def compute_counterfactual(self, x, y_target, features_whitelist=None, regulariz
             preprocessings = self.mymodel.models[:-1]
 
             # Check types
-            if not isinstance(model, ConvexQuadraticProgram):
-                raise TypeError(f"The last model in the pipeline must be an instance of 'ceml.optim.ConvexQuadraticProgram' but not of {type(model)}")
+            if not isinstance(model, SupportAffinePreprocessing):
+                raise TypeError(f"The last model in the pipeline must be an instance of 'ceml.optim.cvx.SupportAffinePreprocessing' but not of {type(model)}")
             if not all([isinstance(m, AffinePreprocessing) for m in preprocessings]):
                 raise TypeError("All models (except the last one) in the pipeline must be an instance of an affine mapping('ceml.backend.jax.AffinePreprocessing')")
 

ceml/sklearn/qda.py

@@ -133,16 +133,18 @@ def _build_constraints(self, var_X, var_x, y):
         i = int(y)
         j = 0 if y == 1 else 1
 
-        A = .5 * ( self.mymodel.sigma_inv[i] - self.mymodel.sigma_inv[j])
+        A = .5 * (self.mymodel.sigma_inv[i] - self.mymodel.sigma_inv[j])
         b = np.dot(self.mymodel.sigma_inv[j], self.mymodel.means[j]) - np.dot(self.mymodel.sigma_inv[i], self.mymodel.means[i])
         c = np.log(self.mymodel.class_priors[j] / self.mymodel.class_priors[i]) + 0.5 * np.log(np.linalg.det(self.mymodel.sigma_inv[j]) / np.linalg.det(self.mymodel.sigma_inv[i])) + 0.5 * (self.mymodel.means[i].T.dot(self.mymodel.sigma_inv[i]).dot(self.mymodel.means[i]) - self.mymodel.means[j].T.dot(self.mymodel.sigma_inv[j]).dot(self.mymodel.means[j]))
 
         return [cp.trace(A @ var_X) + var_x.T @ b + c + self.epsilon <= 0]
 
     def _build_solve_dcqp(self, x_orig, y_target, regularization, features_whitelist):
+        x_orig_prime = self._apply_affine_preprocessing_to_const(x_orig)
+
         Q0 = np.eye(self.mymodel.dim)   # TODO: Can be ignored if regularization != l2
         Q1 = np.zeros((self.mymodel.dim, self.mymodel.dim))
-        q = -2. * x_orig
+        q = -2. * x_orig_prime
         c = 0.0
 
         A0_i = []
@@ -163,7 +165,7 @@ def _build_solve_dcqp(self, x_orig, y_target, regularization, features_whitelist
 
     def solve(self, x_orig, y_target, regularization, features_whitelist, return_as_dict):
         xcf = None
-        if self.mymodel.is_binary:
+        if self.mymodel.is_binary and not self.is_affine_preprocessing_set():
             xcf = self.build_solve_opt(x_orig, y_target, features_whitelist)
         else:
             xcf = self._build_solve_dcqp(x_orig, y_target, regularization, features_whitelist)

ceml/sklearn/softmaxregression.py

@@ -128,7 +128,7 @@ def _build_constraints(self, var_x, y):
         constraints = []
 
         # If set, a apply an affine preprocessing to x
-        var_x_ = self._apply_affine_preprocessing(var_x)
+        var_x_ = self._apply_affine_preprocessing_to_var(var_x)
 
         # Build constraints
         if self.mymodel.is_multiclass is True:  # Multiclass classifier