Improved error handling and Bugfixes of MPs + AffineTransformation

ceml/optim/cvx.py

@@ -64,7 +64,7 @@ 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
+        return self.A is not None and self.b is not None
 
     def _apply_affine_preprocessing_to_var(self, var_x):
         if self.A is not None and self.b is not None:
@@ -90,6 +90,7 @@ class ConvexQuadraticProgram(ABC, SupportAffinePreprocessing):
     def __init__(self, **kwds):
         self.epsilon = 1e-2
         self.solver = cp.SCS
+        self.solver_verbosity = False
 
         super().__init__(**kwds)
 
@@ -112,7 +113,7 @@ def _build_constraints(self, var_x, y):
         raise NotImplementedError()
 
     def _solve(self, prob):
-        prob.solve(solver=self.solver, verbose=False)
+        prob.solve(solver=self.solver, verbose=self.solver_verbosity)
 
     def build_solve_opt(self, x_orig, y, features_whitelist=None, mad=None, optimizer_args=None):
         """Builds and solves the convex quadratic optimization problem.
@@ -152,6 +153,8 @@ def build_solve_opt(self, x_orig, y, features_whitelist=None, mad=None, optimize
                 self.epsilon = optimizer_args["epsilon"]
             if "solver" in optimizer_args:
                 self.solver = cvx_desc_to_solver(optimizer_args["solver"])
+            if "solver_verbosity" in optimizer_args:
+                self.solver_verbosity = optimizer_args["solver_verbosity"]
 
         dim = x_orig.shape[0]
 
@@ -218,6 +221,7 @@ class SDP(ABC):
     def __init__(self, **kwds):
         self.epsilon = 1e-2
         self.solver = cp.SCS
+        self.solver_verbosity = False
 
         super().__init__(**kwds)
 
@@ -242,7 +246,7 @@ def _build_constraints(self, var_X, var_x, y):
         raise NotImplementedError()
 
     def _solve(self, prob):
-        prob.solve(solver=self.solver, verbose=False)
+        prob.solve(solver=self.solver, verbose=self.solver_verbosity)
 
     def build_solve_opt(self, x_orig, y, features_whitelist=None, optimizer_args=None):
         """Builds and solves the SDP.
@@ -276,6 +280,8 @@ def build_solve_opt(self, x_orig, y, features_whitelist=None, optimizer_args=Non
                 self.epsilon = optimizer_args["epsilon"]
             if "solver" in optimizer_args:
                 self.solver = cvx_desc_to_solver(optimizer_args["solver"])
+            if "solver_verbosity" in optimizer_args:
+                self.solver_verbosity = optimizer_args["solver_verbosity"]
 
         dim = x_orig.shape[0]
 
@@ -418,6 +424,7 @@ def __init__(self, model, Q0, Q1, q, c, A0_i, A1_i, b_i, r_i, features_whitelist
         self.mu = 1.5
 
         self.solver = cp.SCS
+        self.solver_verbosity = False
 
         if optimizer_args is not None:
             if "epsilon" in optimizer_args:
@@ -430,14 +437,16 @@ def __init__(self, model, Q0, Q1, q, c, A0_i, A1_i, b_i, r_i, features_whitelist
                 self.mu = optimizer_args["mu"]
             if "solver" in optimizer_args:
                 self.solver = cvx_desc_to_solver(optimizer_args["solver"])
+            if "solver_verbosity" in optimizer_args:
+                self.solver_verbosity = optimizer_args["solver_verbosity"]
 
         if not(len(self.A0s) == len(self.A1s) and len(self.A0s) == len(self.bs) and len(self.rs) == len(self.bs)):
             raise ValueError("Inconsistent number of constraint parameters")
 
         super().__init__(**kwds)
 
     def _solve(self, prob):
-        prob.solve(solver=self.solver, verbose=False)
+        prob.solve(solver=self.solver, verbose=self.solver_verbosity)
 
     def solve_aux(self, xcf, tao, x_orig):
         try:
@@ -459,7 +468,7 @@ def solve_aux(self, xcf, tao, x_orig):
             for i in range(len(self.A0s)):
                 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]
+                c = self.rs[i] + np.dot(xcf, np.dot(xcf, self.A1s[i])) - 2. * var_x_prime.T @ np.dot(xcf, self.A1s[i]) - s[i]
 
                 constraints.append(A + q + c + self.epsilon <= 0)
 
@@ -516,10 +525,14 @@ def compute_counterfactual(self, x_orig, y_target, x0):
 
         # Solve a bunch of CCPs
         while cur_tao < self.tao_max:
-            xcf_ = self.solve_aux(xcf, cur_tao, x_orig)
+            cur_xcf = xcf
+            if cur_xcf.shape == x_orig.shape:   # Apply transformation is necessary - xcf is computed in the original space which can be different from the space the model works on!
+                cur_xcf = self._apply_affine_preprocessing_to_const(cur_xcf)
+
+            xcf_ = self.solve_aux(cur_xcf, cur_tao, x_orig)
             xcf = xcf_
 
-            if y_target == self.model.predict([xcf_])[0]:
+            if y_target == self.model.predict([self._apply_affine_preprocessing_to_const(xcf_)])[0]:
                 break
 
             # Increase penalty parameter

ceml/sklearn/lda.py

@@ -156,7 +156,7 @@ def solve(self, x_orig, y_target, regularization, features_whitelist, return_as_
         xcf = self.build_solve_opt(x_orig, y_target, features_whitelist, mad=mad, optimizer_args=optimizer_args)
         delta = x_orig - xcf
 
-        if self._model_predict([xcf]) != y_target:
+        if self._model_predict([self._apply_affine_preprocessing_to_const(xcf)]) != y_target:
             raise Exception("No counterfactual found - Consider changing parameters 'regularization', 'features_whitelist', 'optimizer' and try again")
 
         if return_as_dict is True:

ceml/sklearn/lvq.py

@@ -181,7 +181,7 @@ def _build_constraints(self, var_x, y):
     def solve(self, target_prototype, features_whitelist=None):
         self.target_prototype = target_prototype
 
-        return self.build_solve_opt(self.x_orig, self.y_target, features_whitelist, mad=None if self.regularization != "l1" else np.ones(self.mymodel.dim), optimizer_args=self.optimizer_args)
+        return self.build_solve_opt(self.x_orig, self.y_target, features_whitelist, mad=None if self.regularization != "l1" else np.ones(self.x_orig.shape[0]), optimizer_args=self.optimizer_args)
 
 
 class LvqCounterfactual(SklearnCounterfactual, MathematicalProgram, DCQP):
@@ -220,8 +220,10 @@ def _compute_counterfactual_via_convex_quadratic_programming(self, x_orig, y_tar
         xcf_dist = float("inf")
 
         dist = lambda x: np.linalg.norm(x - x_orig, 2)
+        dist_ = lambda x: np.linalg.norm(x - self._apply_affine_preprocessing_to_const(x_orig), 2)
         if regularization == "l1":
             dist = lambda x: np.sum(np.abs(x - x_orig))
+            dist_ = lambda x: np.sum(np.abs(x - self._apply_affine_preprocessing_to_const(x_orig)))
 
         # Search for suitable prototypes
         target_prototypes = []
@@ -238,7 +240,9 @@ def _compute_counterfactual_via_convex_quadratic_programming(self, x_orig, y_tar
         for i in range(len(target_prototypes)):
             try:
                 xcf_ = solver.solve(target_prototype=i, features_whitelist=features_whitelist)
-                ycf_ = self.mymodel.model.predict([xcf_])[0]
+                if xcf_ is None:
+                    raise Exception("Optimization algorithm failed - program might be infeasible.\nRerun with 'sovler_verbosity=True' and check the output.")
+                ycf_ = self.mymodel.model.predict([self._apply_affine_preprocessing_to_const(xcf_)])[0]
 
                 if ycf_ == y_target:
                     if dist(xcf_) < xcf_dist:
@@ -249,8 +253,11 @@ def _compute_counterfactual_via_convex_quadratic_programming(self, x_orig, y_tar
 
         if xcf is None:
             # It might happen that the solver (for a specific set of parameter values) does not find a counterfactual, although the feasible region is always non-empty
-            j = np.argmin([dist(self.mymodel.prototypes[proto]) for proto in target_prototypes]) # Select the nearest prototype!
-            xcf = self.mymodel.prototypes[j]
+            if self.is_affine_preprocessing_set() is False:
+                j = np.argmin([dist_(self.mymodel.prototypes[proto]) for proto in target_prototypes]) # Select the nearest prototype!
+                xcf = self.mymodel.prototypes[j]
+            else:
+                raise Exception("Did not find a counterfactual")  # Note: In case of an affine preprocessing, we can not simply select a prototype as a counterfactual!
 
         return xcf
 
@@ -290,7 +297,7 @@ def _build_solve_dcqp(self, x_orig, y_target, target_prototype_id, other_prototy
 
         self.build_program(self.model, x_orig, y_target, Q0, Q1, q, c, A0_i, A1_i, b_i, r_i, features_whitelist=features_whitelist, mad=None if regularization != "l1" else np.ones(self.mymodel.dim))
 
-        return DCQP.solve(self, x0=p_i, tao=1.2, tao_max=100, mu=1.5)
+        return DCQP.solve(self, x0=p_i)
 
     def _compute_counterfactual_via_dcqp(self, x_orig, y_target, features_whitelist, regularization, optimizer_args):
         xcf = None
@@ -313,7 +320,7 @@ def _compute_counterfactual_via_dcqp(self, x_orig, y_target, features_whitelist,
         for i in range(len(target_prototypes)):
             try:
                 xcf_ = self._build_solve_dcqp(x_orig, y_target, i, other_prototypes, features_whitelist, regularization)
-                ycf_ = self.model.predict([xcf_])[0]
+                ycf_ = self.model.predict([self._apply_affine_preprocessing_to_const(xcf_)])[0]
 
                 if ycf_ == y_target:
                     if dist(xcf_) < xcf_dist:
@@ -324,8 +331,11 @@ def _compute_counterfactual_via_dcqp(self, x_orig, y_target, features_whitelist,
 
         if xcf is None:
             # It might happen that the solver (for a specific set of parameter values) does not find a counterfactual, although the feasible region is always non-empty
-            j = np.argmin([dist(self.mymodel.prototypes[proto]) for proto in target_prototypes]) # Select the nearest prototype!
-            xcf = self.mymodel.prototypes[j]
+            if self.is_affine_preprocessing_set() is False:
+                j = np.argmin([dist(self.mymodel.prototypes[proto]) for proto in target_prototypes]) # Select the nearest prototype!
+                xcf = self.mymodel.prototypes[j]
+            else:
+                raise Exception("Did not find a counterfactual")  # Note: In case of an affine preprocessing, we can not simply select a prototype as a counterfactual!
 
         return xcf
 

ceml/sklearn/naivebayes.py

@@ -154,7 +154,7 @@ def _build_solve_dcqp(self, x_orig, y_target, regularization, features_whitelist
             b_i.append((self.mymodel.means[j, :] / self.mymodel.variances[j, :]) - (self.mymodel.means[i, :] / self.mymodel.variances[i, :]))
             r_i.append(np.log(self.mymodel.class_priors[j] / self.mymodel.class_priors[i]) + np.sum([np.log(1. / np.sqrt(2.*np.pi*self.mymodel.variances[j,k])) - ((self.mymodel.means[j,k]**2) / (2.*self.mymodel.variances[j,k])) for k in range(self.mymodel.dim)]) - np.sum([np.log(1. / np.sqrt(2.*np.pi*self.mymodel.variances[i,k])) - ((self.mymodel.means[i,k]**2) / (2.*self.mymodel.variances[i,k])) for k in range(self.mymodel.dim)]))
 
-        self.build_program(self.model, x_orig, y_target, Q0, Q1, q, c, A0_i, A1_i, b_i, r_i, features_whitelist=features_whitelist, mad=None if regularization != "l1" else np.ones(self.mymodel.dim), optimizer_args=optimizer_args)
+        self.build_program(self.model, x_orig, y_target, Q0, Q1, q, c, A0_i, A1_i, b_i, r_i, features_whitelist=features_whitelist, mad=None if regularization != "l1" else np.ones(x_orig.shape[0]), optimizer_args=optimizer_args)
 
         return DCQP.solve(self, x0=self.mymodel.means[i, :])
 

ceml/sklearn/qda.py

@@ -159,7 +159,7 @@ def _build_solve_dcqp(self, x_orig, y_target, regularization, features_whitelist
             b_i.append(np.dot(self.mymodel.sigma_inv[j], self.mymodel.means[j]) - np.dot(self.mymodel.sigma_inv[i], self.mymodel.means[i]))
             r_i.append(np.log(self.mymodel.class_priors[j] / self.mymodel.class_priors[i]) + .5 * np.log(np.linalg.det(self.mymodel.sigma_inv[j]) / np.linalg.det(self.mymodel.sigma_inv[i])) + .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])))
 
-        self.build_program(self.model, x_orig, y_target, Q0, Q1, q, c, A0_i, A1_i, b_i, r_i, features_whitelist=features_whitelist, mad=None if regularization != "l1" else np.ones(self.mymodel.dim), optimizer_args=optimizer_args)
+        self.build_program(self.model, x_orig, y_target, Q0, Q1, q, c, A0_i, A1_i, b_i, r_i, features_whitelist=features_whitelist, mad=None if regularization != "l1" else np.ones(x_orig.shape[0]), optimizer_args=optimizer_args)
 
         return DCQP.solve(self, x0=self.mymodel.means[i])