diff --git a/dask_glm/algorithms.py b/dask_glm/algorithms.py index b7bdbd0..e3d0bf6 100644 --- a/dask_glm/algorithms.py +++ b/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 diff --git a/dask_glm/estimators.py b/dask_glm/estimators.py index 9da2dcb..2ac91ad 100644 --- a/dask_glm/estimators.py +++ b/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) diff --git a/dask_glm/families.py b/dask_glm/families.py index 980b702..d944b5a 100644 --- a/dask_glm/families.py +++ b/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,65 +52,36 @@ 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. @@ -88,32 +89,18 @@ class Poisson(object): .. _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) diff --git a/dask_glm/regularizers.py b/dask_glm/regularizers.py index ccbb385..6719bc0 100644 --- a/dask_glm/regularizers.py +++ b/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.""" diff --git a/dask_glm/tests/test_admm.py b/dask_glm/tests/test_admm.py index 7b0373a..b5f21a3 100644 --- a/dask_glm/tests/test_admm.py +++ b/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)) diff --git a/dask_glm/tests/test_algos_families.py b/dask_glm/tests/test_algos_families.py index e04ffd1..1d3aa74 100644 --- a/dask_glm/tests/test_algos_families.py +++ b/dask_glm/tests/test_algos_families.py @@ -8,9 +8,33 @@ from dask_glm.algorithms import (newton, lbfgs, proximal_grad, gradient_descent, admm) -from dask_glm.families import Logistic, Normal, Poisson +from dask_glm.families import Family from dask_glm.regularizers import Regularizer from dask_glm.utils import sigmoid, make_y +from mock import patch + + +def test_family_pointwise_loss(): + beta = np.array([1, 2]) + X = np.array([[1, 2], [3, 4]]) + y = np.array([1, 2]) + with patch('dask_glm.families.Family.loglikelihood') as loglike: + Family().pointwise_loss(beta, X, y) + Xbeta, new_y = loglike.call_args[0] + np.testing.utils.assert_array_equal(Xbeta, X.dot(beta)) + np.testing.utils.assert_array_equal(new_y, y) + + +def test_family_pointwise_gradient(): + beta = np.array([1, 2]) + X = np.array([[1, 2], [3, 4]]) + y = np.array([1, 2]) + with patch('dask_glm.families.Family.gradient') as gradient: + Family().pointwise_gradient(beta, X, y) + Xbeta, new_x, new_y = gradient.call_args[0] + np.testing.utils.assert_array_equal(Xbeta, X.dot(beta)) + np.testing.utils.assert_array_equal(new_x, X) + np.testing.utils.assert_array_equal(new_y, y) def add_l1(f, lam): @@ -63,7 +87,7 @@ def test_methods(N, p, seed, opt): ]) @pytest.mark.parametrize('N', [1000]) @pytest.mark.parametrize('nchunks', [1, 10]) -@pytest.mark.parametrize('family', [Logistic, Normal, Poisson]) +@pytest.mark.parametrize('family', Family._all_children()) def test_basic_unreg_descent(func, kwargs, N, nchunks, family): beta = np.random.normal(size=2) M = len(beta) @@ -87,9 +111,9 @@ def test_basic_unreg_descent(func, kwargs, N, nchunks, family): ]) @pytest.mark.parametrize('N', [1000]) @pytest.mark.parametrize('nchunks', [1, 10]) -@pytest.mark.parametrize('family', [Logistic, Normal, Poisson]) +@pytest.mark.parametrize('family', Family._all_children()) @pytest.mark.parametrize('lam', [0.01, 1.2, 4.05]) -@pytest.mark.parametrize('reg', [r() for r in Regularizer.__subclasses__()]) +@pytest.mark.parametrize('reg', Regularizer._all_children()) def test_basic_reg_descent(func, kwargs, N, nchunks, family, lam, reg): beta = np.random.normal(size=2) M = len(beta) diff --git a/dask_glm/tests/test_regularizers.py b/dask_glm/tests/test_regularizers.py index 6c875af..5be3715 100644 --- a/dask_glm/tests/test_regularizers.py +++ b/dask_glm/tests/test_regularizers.py @@ -40,22 +40,6 @@ def foo(x): assert new_func(2) == 7 -def test_regularizer_get_passes_through_instance(): - x = FooRegularizer() - assert regs.Regularizer.get(x) == x - - -def test_regularizer_get_unnamed_raises(): - with pytest.raises(KeyError): - regs.Regularizer.get('foo') - - -def test_regularizer_gets_from_name(): - class Foo(regs.Regularizer): - name = 'foo' - assert isinstance(regs.Regularizer.get('foo'), Foo) - - @pytest.mark.parametrize('beta,expected', [ (np.array([0, 0, 0]), 0), (np.array([1, 2, 3]), 7) diff --git a/dask_glm/tests/test_utils.py b/dask_glm/tests/test_utils.py index 813b2fd..08bf360 100644 --- a/dask_glm/tests/test_utils.py +++ b/dask_glm/tests/test_utils.py @@ -78,3 +78,27 @@ def test_sparse(): assert utils.sum(x) == utils.sum(x.todense()) for func in [utils.sigmoid, utils.sum, utils.exp]: assert_eq(func(x), func(y)) + + +@pytest.fixture +def foo_class(): + class FooClass(utils.RegistryClass): + pass + + class BarClass(FooClass): + name = 'bar' + return FooClass + + +def test_registryclass_get_passes_through_instance(foo_class): + x = foo_class() + assert utils.RegistryClass.get(x) == x + + +def test_regularizer_get_unnamed_raises(foo_class): + with pytest.raises(KeyError): + foo_class().get('foobar') + + +def test_regularizer_gets_from_name(foo_class): + assert isinstance(foo_class().get('bar'), foo_class) diff --git a/dask_glm/utils.py b/dask_glm/utils.py index 592c74e..25686c2 100644 --- a/dask_glm/utils.py +++ b/dask_glm/utils.py @@ -194,3 +194,33 @@ def package_of(obj): return base, _sep, _stem = mod.__name__.partition('.') return sys.modules[base] + + +class RegistryClass(object): + """Convenience class for methods on subclasses.""" + + @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 {} object.'.format(cls.__name__)) + + @classmethod + def _all_children(cls): + """"Return a list of instances of all subclasses.""" + return [o() for o in cls.__subclasses__()]