From 1f98c44eeb528f6008fc93cb1d0dce42b26bda42 Mon Sep 17 00:00:00 2001 From: Thomas Wiecki Date: Fri, 25 Nov 2016 11:03:37 +0100 Subject: [PATCH 1/4] ENH Do not compute gradient in find_MAP() if it's not required. Closes #639. --- pymc3/tests/test_starting.py | 19 +++++++++++++-- pymc3/tuning/starting.py | 46 ++++++++++++++++++++++++------------ 2 files changed, 48 insertions(+), 17 deletions(-) diff --git a/pymc3/tests/test_starting.py b/pymc3/tests/test_starting.py index b3719e11f1..a3e21dadcb 100644 --- a/pymc3/tests/test_starting.py +++ b/pymc3/tests/test_starting.py @@ -1,9 +1,10 @@ from .checks import close_to import numpy as np from pymc3.tuning import starting -from pymc3 import Model, Uniform, Normal, Beta, Binomial, find_MAP, Point +from pymc3 import Model, Uniform, Normal, Beta, Binomial, find_MAP, Point, Poisson from .models import simple_model, non_normal, exponential_beta - +from theano.compile.ops import as_op +import theano.tensor as tt def test_accuracy_normal(): _, model, (mu, _) = simple_model() @@ -53,6 +54,20 @@ def test_find_MAP_discrete(): assert map_est2['ss'] == 14 +def test_find_MAP_no_gradient(): + # as_op without gradient definition + @as_op(itypes=[tt.dscalar], otypes=[tt.dscalar]) + def arbitrary_det(value): + return value + + with Model() as model_deterministic: + a = Normal('a') + b = arbitrary_det(a) + c = Normal('obs', mu=b.astype('float64'), + observed=np.array([1, 3, 5])) + find_MAP() + + def test_find_MAP(): tol = 2.0**-11 # 16 bit machine epsilon, a low bar data = np.random.randn(100) diff --git a/pymc3/tuning/starting.py b/pymc3/tuning/starting.py index 823c6a1987..2a1a9a83ca 100644 --- a/pymc3/tuning/starting.py +++ b/pymc3/tuning/starting.py @@ -17,8 +17,8 @@ __all__ = ['find_MAP'] -def find_MAP(start=None, vars=None, fmin=None, return_raw=False, - model=None, *args, **kwargs): +def find_MAP(start=None, vars=None, fmin=None, + return_raw=False,model=None, *args, **kwargs): """ Sets state to the local maximum a posteriori point given a model. Current default of fmin_Hessian does not deal well with optimizing close @@ -55,8 +55,15 @@ def find_MAP(start=None, vars=None, fmin=None, return_raw=False, disc_vars = list(typefilter(vars, discrete_types)) - if disc_vars: - pm._log.warning("Warning: vars contains discrete variables. MAP " + + try: + model.fastdlogp(vars) + gradient_avail = True + except AttributeError: + gradient_avail = False + + if disc_vars or not gradient_avail : + pm._log.warning("Warning: gradient not available." + + "(E.g. vars contains discrete variables). MAP " + "estimates may not be accurate for the default " + "parameters. Defaulting to non-gradient minimization " + "fmin_powell.") @@ -74,19 +81,21 @@ def find_MAP(start=None, vars=None, fmin=None, return_raw=False, bij = DictToArrayBijection(ArrayOrdering(vars), start) logp = bij.mapf(model.fastlogp) - dlogp = bij.mapf(model.fastdlogp(vars)) - def logp_o(point): return nan_to_high(-logp(point)) - def grad_logp_o(point): - return nan_to_num(-dlogp(point)) - # Check to see if minimization function actually uses the gradient if 'fprime' in getargspec(fmin).args: + dlogp = bij.mapf(model.fastdlogp(vars)) + def grad_logp_o(point): + return nan_to_num(-dlogp(point)) + r = fmin(logp_o, bij.map( start), fprime=grad_logp_o, *args, **kwargs) + compute_gradient = True else: + compute_gradient = False + # Check to see if minimization function uses a starting value if 'x0' in getargspec(fmin).args: r = fmin(logp_o, bij.map(start), *args, **kwargs) @@ -100,17 +109,24 @@ def grad_logp_o(point): mx = bij.rmap(mx0) - if (not allfinite(mx0) or - not allfinite(model.logp(mx)) or - not allfinite(model.dlogp()(mx))): + allfinite_mx0 = allfinite(mx0) + allfinite_logp = (model.logp(mx)) + if compute_gradient: + allfinite_dlogp = allfinite(model.dlogp()(mx)) + else: + allfinite_dlogp = True + + if (not allfinite_mx0 or + not allfinite_logp or + not allfinite_dlogp): messages = [] for var in vars: - vals = { "value": mx[var.name], - "logp": var.logp(mx), - "dlogp": var.dlogp()(mx)} + "logp": var.logp(mx)} + if compute_gradient: + vals["dlogp"] = var.dlogp()(mx) def message(name, values): if np.size(values) < 10: From 5d52473ee56741001313c9477c16c3eed53fceeb Mon Sep 17 00:00:00 2001 From: Thomas Wiecki Date: Fri, 25 Nov 2016 11:05:03 +0100 Subject: [PATCH 2/4] MAINT Remove unused Poisson import. --- pymc3/tests/test_starting.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pymc3/tests/test_starting.py b/pymc3/tests/test_starting.py index a3e21dadcb..db312126f5 100644 --- a/pymc3/tests/test_starting.py +++ b/pymc3/tests/test_starting.py @@ -1,7 +1,7 @@ from .checks import close_to import numpy as np from pymc3.tuning import starting -from pymc3 import Model, Uniform, Normal, Beta, Binomial, find_MAP, Point, Poisson +from pymc3 import Model, Uniform, Normal, Beta, Binomial, find_MAP, Point from .models import simple_model, non_normal, exponential_beta from theano.compile.ops import as_op import theano.tensor as tt From 5018a3f207021839d19a055713e82b8807f40343 Mon Sep 17 00:00:00 2001 From: Thomas Wiecki Date: Fri, 25 Nov 2016 11:45:49 +0100 Subject: [PATCH 3/4] BUG Add missing allfinite() call. Refactor model to models.py --- pymc3/tests/models.py | 15 +++++++++++++++ pymc3/tests/test_starting.py | 14 ++------------ pymc3/tuning/starting.py | 2 +- 3 files changed, 18 insertions(+), 13 deletions(-) diff --git a/pymc3/tests/models.py b/pymc3/tests/models.py index 2e3a9b2db3..d8a195fd46 100644 --- a/pymc3/tests/models.py +++ b/pymc3/tests/models.py @@ -3,6 +3,7 @@ import pymc3 as pm from itertools import product import theano.tensor as tt +from theano.compile.ops import as_op def simple_model(): @@ -34,6 +35,20 @@ def multidimensional_model(): return model.test_point, model, (mu, tau ** -1) +def simple_arbitrary_det(): + @as_op(itypes=[tt.dscalar], otypes=[tt.dscalar]) + def arbitrary_det(value): + return value + + with Model() as model: + a = Normal('a') + b = arbitrary_det(a) + c = Normal('obs', mu=b.astype('float64'), + observed=np.array([1, 3, 5])) + + return model.test_point, model + + def simple_init(): start, model, moments = simple_model() step = Metropolis(model.vars, np.diag([1.]), model=model) diff --git a/pymc3/tests/test_starting.py b/pymc3/tests/test_starting.py index db312126f5..c8c2995613 100644 --- a/pymc3/tests/test_starting.py +++ b/pymc3/tests/test_starting.py @@ -3,8 +3,6 @@ from pymc3.tuning import starting from pymc3 import Model, Uniform, Normal, Beta, Binomial, find_MAP, Point from .models import simple_model, non_normal, exponential_beta -from theano.compile.ops import as_op -import theano.tensor as tt def test_accuracy_normal(): _, model, (mu, _) = simple_model() @@ -55,16 +53,8 @@ def test_find_MAP_discrete(): def test_find_MAP_no_gradient(): - # as_op without gradient definition - @as_op(itypes=[tt.dscalar], otypes=[tt.dscalar]) - def arbitrary_det(value): - return value - - with Model() as model_deterministic: - a = Normal('a') - b = arbitrary_det(a) - c = Normal('obs', mu=b.astype('float64'), - observed=np.array([1, 3, 5])) + _, model = simple_arbitrary_det() + with model: find_MAP() diff --git a/pymc3/tuning/starting.py b/pymc3/tuning/starting.py index 2a1a9a83ca..7313f96d51 100644 --- a/pymc3/tuning/starting.py +++ b/pymc3/tuning/starting.py @@ -110,7 +110,7 @@ def grad_logp_o(point): mx = bij.rmap(mx0) allfinite_mx0 = allfinite(mx0) - allfinite_logp = (model.logp(mx)) + allfinite_logp = allfinite(model.logp(mx)) if compute_gradient: allfinite_dlogp = allfinite(model.dlogp()(mx)) else: From 107a483b582770696805dfc3801f50218eb5a5dd Mon Sep 17 00:00:00 2001 From: Thomas Wiecki Date: Fri, 25 Nov 2016 13:43:54 +0100 Subject: [PATCH 4/4] TST Missed import. --- pymc3/tests/test_starting.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pymc3/tests/test_starting.py b/pymc3/tests/test_starting.py index c8c2995613..b97d2114ac 100644 --- a/pymc3/tests/test_starting.py +++ b/pymc3/tests/test_starting.py @@ -2,7 +2,7 @@ import numpy as np from pymc3.tuning import starting from pymc3 import Model, Uniform, Normal, Beta, Binomial, find_MAP, Point -from .models import simple_model, non_normal, exponential_beta +from .models import simple_model, non_normal, exponential_beta, simple_arbitrary_det def test_accuracy_normal(): _, model, (mu, _) = simple_model()