ABC for Families

dask_glm/algorithms.py

@@ -11,12 +11,12 @@
 
 
 from dask_glm.utils import dot, normalize
-from dask_glm.families import Logistic
+from dask_glm.families import Family
 from dask_glm.regularizers import Regularizer
 
 
 def compute_stepsize_dask(beta, step, Xbeta, Xstep, y, curr_val,
-                          family=Logistic, stepSize=1.0,
+                          family='logistic', stepSize=1.0,
                           armijoMult=0.1, backtrackMult=0.1):
     """Compute the optimal stepsize
 
@@ -38,8 +38,8 @@ def compute_stepsize_dask(beta, step, Xbeta, Xstep, y, curr_val,
     xBeta : array-like
     func : callable
     """
-
-    loglike = family.loglike
+    family = Family.get(family)
+    loglike = family.loglikelihood
     beta, step, Xbeta, Xstep, y, curr_val = persist(beta, step, Xbeta, Xstep, y, curr_val)
     obeta, oXbeta = beta, Xbeta
     (step,) = compute(step)
@@ -65,7 +65,7 @@ def compute_stepsize_dask(beta, step, Xbeta, Xstep, y, curr_val,
 
 
 @normalize
-def gradient_descent(X, y, max_iter=100, tol=1e-14, family=Logistic, **kwargs):
+def gradient_descent(X, y, max_iter=100, tol=1e-14, family='logistic', **kwargs):
     """
     Michael Grant's implementation of Gradient Descent.
 
@@ -85,8 +85,8 @@ def gradient_descent(X, y, max_iter=100, tol=1e-14, family=Logistic, **kwargs):
     -------
     beta : array-like, shape (n_features,)
     """
-
-    loglike, gradient = family.loglike, family.gradient
+    family = Family.get(family)
+    loglike, gradient = family.loglikelihood, family.gradient
     n, p = X.shape
     firstBacktrackMult = 0.1
     nextBacktrackMult = 0.5
@@ -138,7 +138,7 @@ def gradient_descent(X, y, max_iter=100, tol=1e-14, family=Logistic, **kwargs):
 
 
 @normalize
-def newton(X, y, max_iter=50, tol=1e-8, family=Logistic, **kwargs):
+def newton(X, y, max_iter=50, tol=1e-8, family='logistic', **kwargs):
     """Newtons Method for Logistic Regression.
 
     Parameters
@@ -157,6 +157,7 @@ def newton(X, y, max_iter=50, tol=1e-8, family=Logistic, **kwargs):
     -------
     beta : array-like, shape (n_features,)
     """
+    family = Family.get(family)
     gradient, hessian = family.gradient, family.hessian
     n, p = X.shape
     beta = np.zeros(p)  # always init to zeros?
@@ -194,7 +195,7 @@ def newton(X, y, max_iter=50, tol=1e-8, family=Logistic, **kwargs):
 
 @normalize
 def admm(X, y, regularizer='l1', lamduh=0.1, rho=1, over_relax=1,
-         max_iter=250, abstol=1e-4, reltol=1e-2, family=Logistic, **kwargs):
+         max_iter=250, abstol=1e-4, reltol=1e-2, family='logistic', **kwargs):
     """
     Alternating Direction Method of Multipliers
 
@@ -216,7 +217,7 @@ def admm(X, y, regularizer='l1', lamduh=0.1, rho=1, over_relax=1,
     -------
     beta : array-like, shape (n_features,)
     """
-
+    family = Family.get(family)
     pointwise_loss = family.pointwise_loss
     pointwise_gradient = family.pointwise_gradient
     regularizer = Regularizer.get(regularizer)
@@ -303,7 +304,7 @@ def local_update(X, y, beta, z, u, rho, f, fprime, solver=fmin_l_bfgs_b):
 
 @normalize
 def lbfgs(X, y, regularizer=None, lamduh=1.0, max_iter=100, tol=1e-4,
-          family=Logistic, verbose=False, **kwargs):
+          family='logistic', verbose=False, **kwargs):
     """L-BFGS solver using scipy.optimize implementation
 
     Parameters
@@ -322,7 +323,7 @@ def lbfgs(X, y, regularizer=None, lamduh=1.0, max_iter=100, tol=1e-4,
     -------
     beta : array-like, shape (n_features,)
     """
-
+    family = Family.get(family)
     pointwise_loss = family.pointwise_loss
     pointwise_gradient = family.pointwise_gradient
     if regularizer is not None:
@@ -349,7 +350,7 @@ def compute_loss_grad(beta, X, y):
 
 
 @normalize
-def proximal_grad(X, y, regularizer='l1', lamduh=0.1, family=Logistic,
+def proximal_grad(X, y, regularizer='l1', lamduh=0.1, family='logistic',
                   max_iter=100, tol=1e-8, **kwargs):
     """
 
@@ -371,7 +372,7 @@ def proximal_grad(X, y, regularizer='l1', lamduh=0.1, family=Logistic,
     -------
     beta : array-like, shape (n_features,)
     """
-
+    family = Family.get(family)
     n, p = X.shape
     firstBacktrackMult = 0.1
     nextBacktrackMult = 0.5
@@ -387,7 +388,7 @@ def proximal_grad(X, y, regularizer='l1', lamduh=0.1, family=Logistic,
         # Compute the gradient
         if k % recalcRate == 0:
             Xbeta = X.dot(beta)
-            func = family.loglike(Xbeta, y)
+            func = family.loglikelihood(Xbeta, y)
 
         gradient = family.gradient(Xbeta, X, y)
 
@@ -405,7 +406,7 @@ def proximal_grad(X, y, regularizer='l1', lamduh=0.1, family=Logistic,
 
             Xbeta, beta = persist(Xbeta, beta)
 
-            func = family.loglike(Xbeta, y)
+            func = family.loglikelihood(Xbeta, y)
             func = persist(func)[0]
             func = compute(func)[0]
             df = lf - func

dask_glm/estimators.py

@@ -112,7 +112,7 @@ class LogisticRegression(_GLM):
 
     @property
     def family(self):
-        return families.Logistic
+        return families.Logistic()
 
     def predict(self, X):
         return self.predict_proba(X) > .5  # TODO: verify, multiclass broken
@@ -167,7 +167,7 @@ class LinearRegression(_GLM):
     """
     @property
     def family(self):
-        return families.Normal
+        return families.Normal()
 
     def predict(self, X):
         X_ = self._maybe_add_intercept(X)
@@ -219,7 +219,7 @@ class PoissonRegression(_GLM):
     """
     @property
     def family(self):
-        return families.Poisson
+        return families.Poisson()
 
     def predict(self, X):
         X_ = self._maybe_add_intercept(X)

dask_glm/families.py

@@ -1,18 +1,48 @@
 from __future__ import absolute_import, division, print_function
 
-from dask_glm.utils import dot, exp, log1p, sigmoid
+from dask_glm.utils import dot, exp, log1p, sigmoid, RegistryClass
 
 
-class Logistic(object):
+class Family(RegistryClass):
+    """Base class methods for distribution.
+
+    This class represents the required methods to add a new distribution to work
+    with the algorithms.
+    """
+
+    def loglikelihood(self, Xbeta, y):
+        raise NotImplementedError
+
+    def pointwise_loss(self, beta, X, y):
+        """Loss, evaluated point-wise."""
+        beta, y = beta.ravel(), y.ravel()
+        Xbeta = X.dot(beta)
+        return self.loglikelihood(Xbeta, y)
+
+    def pointwise_gradient(self, beta, X, y):
+        """Loss, evaluated point-wise."""
+        beta, y = beta.ravel(), y.ravel()
+        Xbeta = X.dot(beta)
+        return self.gradient(Xbeta, X, y)
+
+    def gradient(self, Xbeta, x, y):
+        raise NotImplementedError
+
+    def hessian(self, Xbeta, x):
+        raise NotImplementedError
+
+
+class Logistic(Family):
     """Implements methods for `Logistic regression`_,
     useful for classifying binary outcomes.
 
     .. _Logistic regression: https://en.wikipedia.org/wiki/Logistic_regression
     """
-    @staticmethod
-    def loglike(Xbeta, y):
+    name = 'logistic'
+
+    def loglikelihood(self, Xbeta, y):
         """
-        Evaluate the logistic loglikeliehood
+        Evaluate the logistic loglikelihood
 
         Parameters
         ----------
@@ -22,98 +52,55 @@ def loglike(Xbeta, y):
         enXbeta = exp(-Xbeta)
         return (Xbeta + log1p(enXbeta)).sum() - dot(y, Xbeta)
 
-    @staticmethod
-    def pointwise_loss(beta, X, y):
-        """Logistic Loss, evaluated point-wise."""
-        beta, y = beta.ravel(), y.ravel()
-        Xbeta = X.dot(beta)
-        return Logistic.loglike(Xbeta, y)
-
-    @staticmethod
-    def pointwise_gradient(beta, X, y):
-        """Logistic gradient, evaluated point-wise."""
-        beta, y = beta.ravel(), y.ravel()
-        Xbeta = X.dot(beta)
-        return Logistic.gradient(Xbeta, X, y)
-
-    @staticmethod
-    def gradient(Xbeta, X, y):
+    def gradient(self, Xbeta, X, y):
         """Logistic gradient"""
         p = sigmoid(Xbeta)
         return dot(X.T, p - y)
 
-    @staticmethod
-    def hessian(Xbeta, X):
+    def hessian(self, Xbeta, X):
         """Logistic hessian"""
         p = sigmoid(Xbeta)
         return dot(p * (1 - p) * X.T, X)
 
 
-class Normal(object):
+class Normal(Family):
     """Implements methods for `Linear regression`_,
     useful for modeling continuous outcomes.
 
     .. _Linear regression: https://en.wikipedia.org/wiki/Linear_regression
     """
-    @staticmethod
-    def loglike(Xbeta, y):
-        return ((y - Xbeta) ** 2).sum()
+    name = 'normal'
 
-    @staticmethod
-    def pointwise_loss(beta, X, y):
-        beta, y = beta.ravel(), y.ravel()
-        Xbeta = X.dot(beta)
-        return Normal.loglike(Xbeta, y)
-
-    @staticmethod
-    def pointwise_gradient(beta, X, y):
-        beta, y = beta.ravel(), y.ravel()
-        Xbeta = X.dot(beta)
-        return Normal.gradient(Xbeta, X, y)
+    def loglikelihood(self, Xbeta, y):
+        return ((y - Xbeta) ** 2).sum()
 
-    @staticmethod
-    def gradient(Xbeta, X, y):
+    def gradient(self, Xbeta, X, y):
         return 2 * dot(X.T, Xbeta) - 2 * dot(X.T, y)
 
-    @staticmethod
-    def hessian(Xbeta, X):
+    def hessian(self, Xbeta, X):
         return 2 * dot(X.T, X)
 
 
-class Poisson(object):
+class Poisson(Family):
     """
     This implements `Poisson regression`_, useful for
     modelling count data.
 
 
     .. _Poisson regression: https://en.wikipedia.org/wiki/Poisson_regression
     """
+    name = 'poisson'
 
-    @staticmethod
-    def loglike(Xbeta, y):
+    def loglikelihood(self, Xbeta, y):
         eXbeta = exp(Xbeta)
         yXbeta = y * Xbeta
         return (eXbeta - yXbeta).sum()
 
-    @staticmethod
-    def pointwise_loss(beta, X, y):
-        beta, y = beta.ravel(), y.ravel()
-        Xbeta = X.dot(beta)
-        return Poisson.loglike(Xbeta, y)
-
-    @staticmethod
-    def pointwise_gradient(beta, X, y):
-        beta, y = beta.ravel(), y.ravel()
-        Xbeta = X.dot(beta)
-        return Poisson.gradient(Xbeta, X, y)
-
-    @staticmethod
-    def gradient(Xbeta, X, y):
+    def gradient(self, Xbeta, X, y):
         eXbeta = exp(Xbeta)
         return dot(X.T, eXbeta - y)
 
-    @staticmethod
-    def hessian(Xbeta, X):
+    def hessian(self, Xbeta, X):
         eXbeta = exp(Xbeta)
         x_diag_eXbeta = eXbeta[:, None] * X
         return dot(X.T, x_diag_eXbeta)

dask_glm/regularizers.py

@@ -1,9 +1,10 @@
 from __future__ import absolute_import, division, print_function
 
 import numpy as np
+from dask_glm.utils import RegistryClass
 
 
-class Regularizer(object):
+class Regularizer(RegistryClass):
     """Abstract base class for regularization object.
 
     Defines the set of methods required to create a new regularization object. This includes
@@ -121,27 +122,6 @@ def wrapped(beta, *args):
             return hess(beta, *args) + lam * self.hessian(beta)
         return wrapped
 
-    @classmethod
-    def get(cls, obj):
-        """Get the concrete instance for the name ``obj``.
-
-        Parameters
-        ----------
-        obj : Regularizer or str
-            Valid instances of ``Regularizer`` are passed through.
-            Strings are looked up according to ``obj.name`` and a
-            new instance is created
-
-        Returns
-        -------
-        obj : Regularizer
-        """
-        if isinstance(obj, cls):
-            return obj
-        elif isinstance(obj, str):
-            return {o.name: o for o in cls.__subclasses__()}[obj]()
-        raise TypeError('Not a valid regularizer object.')
-
 
 class L2(Regularizer):
     """L2 regularization."""

dask_glm/tests/test_admm.py

@@ -15,7 +15,7 @@
                          [np.array([-1.5, 3]),
                           np.array([35, 2, 0, -3.2]),
                           np.array([-1e-2, 1e-4, 1.0, 2e-3, -1.2])])
-@pytest.mark.parametrize('family', [Logistic, Normal])
+@pytest.mark.parametrize('family', [Logistic(), Normal()])
 def test_local_update(N, beta, family):
     M = beta.shape[0]
     X = np.random.random((N, M))