Add support to SCS solver

clogistic/__init__.py

@@ -1,5 +1,6 @@
+from ._version import __version__
 from .clogistic import LogisticRegression
 
-__version__ = "0.1.0"
 
-__all__ = ['LogisticRegression']
+__all__ = ['__version__',
+           'LogisticRegression']

clogistic/_version.py

@@ -0,0 +1,3 @@
+"""Version information."""
+
+__version__ = "0.2.0"

clogistic/clogistic.py

@@ -65,9 +65,9 @@ def _check_parameters(penalty, tol, C, fit_intercept, class_weight, solver,
             raise ValueError('Invalid value for class_weight. Allowed string '
                              'value is "balanced".')
 
-    if solver not in ("ecos", "lbfgs"):
+    if solver not in ("ecos", "lbfgs", "scs"):
         raise ValueError('Invalid value for solver. Allowed string '
-                         'values are "ecos" and "lbfgs".')
+                         'values are "ecos", "lbfgs" and "scs".')
 
     if not isinstance(max_iter, numbers.Number) or max_iter < 0:
         raise ValueError("max_iter must be positive; got {}.".format(max_iter))
@@ -85,10 +85,12 @@ def _check_solver(solver, penalty, bounds, constraints, warm_start):
         if penalty in ("l1", "elasticnet") and bounds is not None:
             raise ValueError('Solver "lbfgs" does not support bound '
                              'constraints with penalty "l1" and '
-                             '"elasticnet"; choose "ecos" solver.')
+                             '"elasticnet"; choose either "ecos" or "scs" '
+                             'solver.')
 
         if constraints is not None:
-            raise ValueError('Only "ecos" solver supports linear constraints.')
+            raise ValueError('"lbfgs" solver does not  supports linear '
+                             'constraints.')
 
         if penalty in ("l1", "elasticnet") and warm_start:
             raise ValueError('Solver "lbfgs" does not support warm start with '
@@ -269,9 +271,9 @@ def _fit_lbfgs(penalty, tol, C, fit_intercept, max_iter, l1_ratio,
     return coef_, intercept_
 
 
-def _fit_ecos(penalty, tol, C, fit_intercept, max_iter, l1_ratio,
-              warm_start_coef, verbose, X, y, sample_weight, bounds=None,
-              constraints=None):
+def _fit_cvxpy(solver, penalty, tol, C, fit_intercept, max_iter, l1_ratio,
+               warm_start_coef, verbose, X, y, sample_weight, bounds=None,
+               constraints=None):
 
     m, n = X.shape
 
@@ -343,8 +345,14 @@ def _fit_ecos(penalty, tol, C, fit_intercept, max_iter, l1_ratio,
 
     obj = cp.Maximize(log_likelihood - penalty)
     problem = cp.Problem(obj, cons)
-    problem.solve(max_iters=max_iter, verbose=verbose, abstol=tol,
-                  warm_start=warm_start)
+
+    if solver == "ecos":
+        solve_options = {'solver': cp.ECOS, 'abstol': tol}
+    elif solver == "scs":
+        solve_options = {'solver': cp.SCS, 'eps': tol}
+
+    problem.solve(max_iters=max_iter, verbose=verbose, warm_start=warm_start,
+                  **solve_options)
 
     if fit_intercept:
         intercept_ = beta.value[-1]
@@ -362,12 +370,12 @@ class LogisticRegression(BaseEstimator, LinearClassifierMixin,
     Constrained Logistic Regression (aka logit, MaxEnt) classifier.
 
     This class implements regularized logistic regression supported bound
-    and linear constraints using the 'ecos' and 'lbfgs' solvers.
+    and linear constraints using the 'ecos', 'scs' and 'lbfgs' solvers.
 
     All solvers support only L1, L2 and Elastic-Net regularization or no
     regularization. The 'lbfgs' solver supports bound constraints for L2
-    regularization. The 'ecos' solver supports bound constraints and linear
-    constraints for all regularizations.
+    regularization. The 'ecos' and 'scs' solver support bound constraints and
+    linear constraints for all regularizations.
 
     Parameters
     ----------
@@ -399,14 +407,15 @@ class LogisticRegression(BaseEstimator, LinearClassifierMixin,
         Note that these weights will be multiplied with sample_weight (passed
         through the fit method) if sample_weight is specified.
 
-    solver : {'ecos', 'lbfgs'}, default='lbfgs'
+    solver : {'ecos', 'lbfgs', 'scs'}, default='ecos'
         Algorithm/solver to use in the optimization problem.
 
         - Unconstrained 'lbfgs' handles all regularizations.
         - Bound constrainted 'lbfgs' handles L2 or no penalty.
-        - For other cases, use 'ecos'.
+        - For other cases, use 'ecos' or 'scs'.
 
-        Note that 'ecos' uses the general-purpose solver ECOS via CVXPY.
+        Note that 'ecos' and 'scs' are general-purpose solvers called via
+        CVXPY.
 
     max_iter : int, default=100
         Maximum number of iterations taken for the solvers to converge.
@@ -520,11 +529,12 @@ def fit(self, X, y, sample_weight=None, bounds=None, constraints=None):
                                         self.intercept_[:, np.newaxis],
                                         axis=1)
 
-        if self.solver == "ecos":
-            coef_, intercept_ = _fit_ecos(
-                self.penalty, self.tol, self.C, self.fit_intercept,
-                self.max_iter, self.l1_ratio, warm_start_coef, self.verbose,
-                X, y, sample_weight, bounds, constraints)
+        if self.solver in ("ecos", "scs"):
+            coef_, intercept_ = _fit_cvxpy(
+                self.solver, self.penalty, self.tol, self.C,
+                self.fit_intercept, self.max_iter, self.l1_ratio,
+                warm_start_coef, self.verbose, X, y, sample_weight, bounds,
+                constraints)
         else:
             coef_, intercept_ = _fit_lbfgs(
                 self.penalty, self.tol, self.C, self.fit_intercept,

setup.py

@@ -24,9 +24,15 @@ def read(fname):
 ]
 
 
+# Read version file
+version_info = {}
+with open("clogistic/_version.py") as f:
+    exec(f.read(), version_info)
+
+
 setup(
     name='clogistic',
-    version='0.1.0',
+    version=version_info['__version__'],
     description="Constrained Logistic Regression",
     long_description=long_description,
     author="Guillermo Navas-Palencia",

tests/test_logistic.py

@@ -151,31 +151,25 @@ def test_predict_breast_cancer():
 
     # Test that all solvers with all regularizations score (>0.93) for the
     # training data
-    clf_none_lbfgs = LogisticRegression(solver="lbfgs", penalty="none")
-    clf_l1_lbfgs = LogisticRegression(solver="lbfgs", penalty="l1")
-    clf_l2_lbfgs = LogisticRegression(solver="lbfgs", penalty="l2")
-    clf_en_lbfgs = LogisticRegression(solver="lbfgs", penalty="elasticnet",
-                                      l1_ratio=0.5)
-    clf_none_ecos = LogisticRegression(solver="ecos", penalty="none")
-    clf_l1_ecos = LogisticRegression(solver="ecos", penalty="l1")
-    clf_l2_ecos = LogisticRegression(solver="ecos", penalty="l2")
-    clf_en_ecos = LogisticRegression(solver="ecos", penalty="elasticnet",
-                                     l1_ratio=0.5)
-
-    for clf in (clf_none_lbfgs, clf_l1_lbfgs, clf_l2_lbfgs, clf_en_lbfgs,
-                clf_none_ecos, clf_l1_ecos, clf_l2_ecos, clf_en_ecos):
+    for solver in ("lbfgs", "ecos", "scs"):
+        for penalty in ("none", "l1", "l2", "elasticnet"):
+            if penalty == "elasticnet":
+                clf = LogisticRegression(solver=solver, penalty=penalty,
+                                         l1_ratio=0.5)
+            else:
+                clf = LogisticRegression(solver=solver, penalty=penalty)
 
-        clf.fit(X, y)
-        assert np.all(np.unique(y) == clf.classes_)
+            clf.fit(X, y)
+            assert np.all(np.unique(y) == clf.classes_)
 
-        pred = clf.predict(X)
-        assert np.mean(pred == y) > 0.93
+            pred = clf.predict(X)
+            assert np.mean(pred == y) > 0.93
 
-        probabilities = clf.predict_proba(X)
-        assert probabilities.sum(axis=1) == approx(np.ones(X.shape[0]))
+            probabilities = clf.predict_proba(X)
+            assert probabilities.sum(axis=1) == approx(np.ones(X.shape[0]))
 
-        pred = clf.classes_[np.argmax(clf.predict_log_proba(X), axis=1)]
-        assert np.mean(pred == y) > .9
+            pred = clf.classes_[np.argmax(clf.predict_log_proba(X), axis=1)]
+            assert np.mean(pred == y) > .9
 
 
 def test_predict_breast_cancer_no_intercept():
@@ -184,33 +178,26 @@ def test_predict_breast_cancer_no_intercept():
 
     # Test that all solvers with all regularizations score (>0.93) for the
     # training data without intercept
-    clf_l1_lbfgs = LogisticRegression(solver="lbfgs", penalty="l1",
-                                      fit_intercept=False)
-    clf_l2_lbfgs = LogisticRegression(solver="lbfgs", penalty="l2",
-                                      fit_intercept=False)
-    clf_en_lbfgs = LogisticRegression(solver="lbfgs", penalty="elasticnet",
-                                      fit_intercept=False, l1_ratio=0.5)
-    clf_l1_ecos = LogisticRegression(solver="ecos", penalty="l1",
-                                     fit_intercept=False)
-    clf_l2_ecos = LogisticRegression(solver="ecos", penalty="l2",
-                                     fit_intercept=False)
-    clf_en_ecos = LogisticRegression(solver="ecos", penalty="elasticnet",
-                                     fit_intercept=False, l1_ratio=0.5)
-
-    for clf in (clf_l1_lbfgs, clf_l2_lbfgs, clf_en_lbfgs, clf_l1_ecos,
-                clf_l2_ecos, clf_en_ecos):
+    for solver in ("lbfgs", "ecos"):
+        for penalty in ("none", "l1", "l2", "elasticnet"):
+            if penalty == "elasticnet":
+                clf = LogisticRegression(solver=solver, penalty=penalty,
+                                         l1_ratio=0.5, fit_intercept=False)
+            else:
+                clf = LogisticRegression(solver=solver, penalty=penalty,
+                                         fit_intercept=False)
 
-        clf.fit(X, y)
-        assert np.all(np.unique(y) == clf.classes_)
+            clf.fit(X, y)
+            assert np.all(np.unique(y) == clf.classes_)
 
-        pred = clf.predict(X)
-        assert np.mean(pred == y) > 0.93
+            pred = clf.predict(X)
+            assert np.mean(pred == y) > 0.93
 
-        probabilities = clf.predict_proba(X)
-        assert probabilities.sum(axis=1) == approx(np.ones(X.shape[0]))
+            probabilities = clf.predict_proba(X)
+            assert probabilities.sum(axis=1) == approx(np.ones(X.shape[0]))
 
-        pred = clf.classes_[np.argmax(clf.predict_log_proba(X), axis=1)]
-        assert np.mean(pred == y) > .9
+            pred = clf.classes_[np.argmax(clf.predict_log_proba(X), axis=1)]
+            assert np.mean(pred == y) > .9
 
 
 def test_predict_breast_cancer_bounds_constraints():
@@ -229,25 +216,25 @@ def test_predict_breast_cancer_bounds_constraints():
 
     # Test that all solvers with all regularizations score (>0.93) for the
     # training data
-    clf_none_ecos = LogisticRegression(solver="ecos", penalty="none")
-    clf_l1_ecos = LogisticRegression(solver="ecos", penalty="l1")
-    clf_l2_ecos = LogisticRegression(solver="ecos", penalty="l2")
-    clf_en_ecos = LogisticRegression(solver="ecos", penalty="elasticnet",
-                                     l1_ratio=0.5)
+    for solver in ("ecos", "scs"):
+        for penalty in ("none", "l1", "l2", "elasticnet"):
+            if penalty == "elasticnet":
+                clf = LogisticRegression(solver=solver, penalty=penalty,
+                                         l1_ratio=0.5)
+            else:
+                clf = LogisticRegression(solver=solver, penalty=penalty)
 
-    for clf in (clf_none_ecos, clf_l1_ecos, clf_l2_ecos, clf_en_ecos):
+            clf.fit(X, y, bounds=bounds, constraints=constraints)
+            assert np.all(np.unique(y) == clf.classes_)
 
-        clf.fit(X, y, bounds=bounds, constraints=constraints)
-        assert np.all(np.unique(y) == clf.classes_)
-
-        pred = clf.predict(X)
-        assert np.mean(pred == y) > 0.93
+            pred = clf.predict(X)
+            assert np.mean(pred == y) > 0.93
 
-        probabilities = clf.predict_proba(X)
-        assert probabilities.sum(axis=1) == approx(np.ones(X.shape[0]))
+            probabilities = clf.predict_proba(X)
+            assert probabilities.sum(axis=1) == approx(np.ones(X.shape[0]))
 
-        pred = clf.classes_[np.argmax(clf.predict_log_proba(X), axis=1)]
-        assert np.mean(pred == y) > .9
+            pred = clf.classes_[np.argmax(clf.predict_log_proba(X), axis=1)]
+            assert np.mean(pred == y) > .9
 
 
 def test_warm_start():