diff --git a/pymc3/data.py b/pymc3/data.py index 944ce55458..d40d58ca0a 100644 --- a/pymc3/data.py +++ b/pymc3/data.py @@ -1,7 +1,10 @@ +import itertools import pkgutil import io -__all__ = ['get_data_file'] +import theano.tensor as tt + +__all__ = ['get_data_file', 'DataGenerator'] def get_data_file(pkg, path): @@ -19,3 +22,30 @@ def get_data_file(pkg, path): """ return io.BytesIO(pkgutil.get_data(pkg, path)) + + +class DataGenerator(object): + """ + Helper class that helps to infer data type of generator with looking + at the first item, preserving the order of the resulting generator + """ + def __init__(self, generator): + self.test_value = next(generator) + self.gen = itertools.chain([self.test_value], generator) + self.tensortype = tt.TensorType( + self.test_value.dtype, + ((False, ) * self.test_value.ndim)) + + def __next__(self): + return next(self.gen) + + next = __next__ + + def __iter__(self): + return self + + def __eq__(self, other): + return id(self) == id(other) + + def __hash__(self): + return hash(id(self)) diff --git a/pymc3/distributions/distribution.py b/pymc3/distributions/distribution.py index 2ad2243f74..cae9740740 100644 --- a/pymc3/distributions/distribution.py +++ b/pymc3/distributions/distribution.py @@ -30,8 +30,9 @@ def __new__(cls, name, *args, **kwargs): if isinstance(name, string_types): data = kwargs.pop('observed', None) + total_size = kwargs.pop('total_size', None) dist = cls.dist(*args, **kwargs) - return model.Var(name, dist, data) + return model.Var(name, dist, data, total_size) else: raise TypeError("Name needs to be a string but got: %s" % name) diff --git a/pymc3/examples/arbitrary_stochastic.py b/pymc3/examples/arbitrary_stochastic.py index 9739c56dac..a342a453e5 100644 --- a/pymc3/examples/arbitrary_stochastic.py +++ b/pymc3/examples/arbitrary_stochastic.py @@ -20,7 +20,7 @@ def run(n_samples=3000): start = model.test_point h = pm.find_hessian(start, model=model) step = pm.Metropolis(model.vars, h, blocked=True, model=model) - trace = pm.sample(n_samples, step, start, model=model) + trace = pm.sample(n_samples, step=step, start=start, model=model) return trace if __name__ == "__main__": diff --git a/pymc3/examples/factor_potential.py b/pymc3/examples/factor_potential.py index 17b1d6c743..a189d729e2 100644 --- a/pymc3/examples/factor_potential.py +++ b/pymc3/examples/factor_potential.py @@ -13,7 +13,7 @@ def run(n=3000): if n == "short": n = 50 with model: - trace = sample(n, step, start) + trace = sample(n, step=step, start=start) if __name__ == '__main__': run() diff --git a/pymc3/model.py b/pymc3/model.py index 6017484d06..0a0d5a04b4 100644 --- a/pymc3/model.py +++ b/pymc3/model.py @@ -1,5 +1,6 @@ import threading import six +import types import numpy as np import theano @@ -8,7 +9,7 @@ import pymc3 as pm from .memoize import memoize -from .theanof import gradient, hessian, inputvars +from .theanof import gradient, hessian, inputvars, generator from .vartypes import typefilter, discrete_types, continuous_types from .blocking import DictToArrayBijection, ArrayOrdering @@ -458,7 +459,7 @@ def cont_vars(self): """All the continuous variables in the model""" return list(typefilter(self.vars, continuous_types)) - def Var(self, name, dist, data=None): + def Var(self, name, dist, data=None, total_size=None): """Create and add (un)observed random variable to the model with an appropriate prior distribution. @@ -469,6 +470,8 @@ def Var(self, name, dist, data=None): data : array_like (optional) If data is provided, the variable is observed. If None, the variable is unobserved. + total_size : scalar + upscales logp of variable with :math:`coef = total_size/var.shape[0]` Returns ------- @@ -477,11 +480,13 @@ def Var(self, name, dist, data=None): name = self.name_for(name) if data is None: if getattr(dist, "transform", None) is None: - var = FreeRV(name=name, distribution=dist, model=self) + var = FreeRV(name=name, distribution=dist, model=self, + total_size=total_size) self.free_RVs.append(var) else: var = TransformedRV(name=name, distribution=dist, model=self, - transform=dist.transform) + transform=dist.transform, + total_size=total_size) pm._log.debug('Applied {transform}-transform to {name}' ' and added transformed {orig_name} to model.'.format( transform=dist.transform.name, @@ -491,7 +496,7 @@ def Var(self, name, dist, data=None): return var elif isinstance(data, dict): var = MultiObservedRV(name=name, data=data, distribution=dist, - model=self) + model=self, total_size=total_size) self.observed_RVs.append(var) if var.missing_values: self.free_RVs += var.missing_values @@ -500,7 +505,8 @@ def Var(self, name, dist, data=None): self.named_vars[v.name] = v else: var = ObservedRV(name=name, data=data, - distribution=dist, model=self) + distribution=dist, model=self, + total_size=total_size) self.observed_RVs.append(var) if var.missing_values: self.free_RVs.append(var.missing_values) @@ -717,7 +723,7 @@ class FreeRV(Factor, TensorVariable): """Unobserved random variable that a model is specified in terms of.""" def __init__(self, type=None, owner=None, index=None, name=None, - distribution=None, model=None): + distribution=None, model=None, total_size=None): """ Parameters ---------- @@ -725,7 +731,10 @@ def __init__(self, type=None, owner=None, index=None, name=None, owner : theano owner (optional) name : str distribution : Distribution - model : Model""" + model : Model + total_size : scalar Tensor (optional) + needed for upscaling logp + """ if type is None: type = distribution.type super(FreeRV, self).__init__(type, owner, index, name) @@ -736,7 +745,14 @@ def __init__(self, type=None, owner=None, index=None, name=None, self.distribution = distribution self.tag.test_value = np.ones( distribution.shape, distribution.dtype) * distribution.default() - self.logp_elemwiset = distribution.logp(self) + logp_elemwiset = distribution.logp(self) + if total_size is None: + coef = tt.constant(1) + else: + assert logp_elemwiset.ndim >= 1, ('Variable with scaled density ' + 'needs to be at least 1 dimensional') + coef = tt.as_tensor(total_size) / logp_elemwiset.shape[0] + self.logp_elemwiset = logp_elemwiset * coef self.model = model incorporate_methods(source=distribution, destination=self, @@ -759,6 +775,9 @@ def pandas_to_array(data): return data elif isinstance(data, theano.gof.graph.Variable): return data + elif (hasattr(data, '__next__') or + isinstance(data, types.GeneratorType)): + return generator(data) else: return np.asarray(data) @@ -792,7 +811,7 @@ class ObservedRV(Factor, TensorVariable): """ def __init__(self, type=None, owner=None, index=None, name=None, data=None, - distribution=None, model=None): + distribution=None, model=None, total_size=None): """ Parameters ---------- @@ -801,6 +820,8 @@ def __init__(self, type=None, owner=None, index=None, name=None, data=None, name : str distribution : Distribution model : Model + total_size : scalar Tensor (optional) + needed for upscaling logp """ from .distributions import TensorType if type is None: @@ -814,7 +835,14 @@ def __init__(self, type=None, owner=None, index=None, name=None, data=None, data = as_tensor(data, name, model, distribution) self.missing_values = data.missing_values - self.logp_elemwiset = distribution.logp(data) + logp_elemwiset = distribution.logp(data) + if total_size is None: + coef = tt.constant(1) + else: + assert logp_elemwiset.ndim >= 1, ('Variable with scaled density ' + 'needs to be at least 1 dimensional') + coef = tt.as_tensor(total_size) / logp_elemwiset.shape[0] + self.logp_elemwiset = logp_elemwiset * coef self.model = model self.distribution = distribution @@ -835,7 +863,7 @@ class MultiObservedRV(Factor): Potentially partially observed. """ - def __init__(self, name, data, distribution, model): + def __init__(self, name, data, distribution, model, total_size=None): """ Parameters ---------- @@ -844,6 +872,8 @@ def __init__(self, name, data, distribution, model): name : str distribution : Distribution model : Model + total_size : scalar Tensor (optional) + needed for upscaling logp """ self.name = name self.data = {name: as_tensor(data, name, model, distribution) @@ -851,7 +881,14 @@ def __init__(self, name, data, distribution, model): self.missing_values = [datum.missing_values for datum in self.data.values() if datum.missing_values is not None] - self.logp_elemwiset = distribution.logp(**self.data) + logp_elemwiset = distribution.logp(**self.data) + if total_size is None: + coef = tt.constant(1) + else: + assert logp_elemwiset.ndim >= 1, ('Variable with scaled density ' + 'needs to be at least 1 dimensional') + coef = tt.as_tensor(total_size) / logp_elemwiset.shape[0] + self.logp_elemwiset = logp_elemwiset * coef self.model = model self.distribution = distribution @@ -896,17 +933,20 @@ def Potential(name, var, model=None): class TransformedRV(TensorVariable): def __init__(self, type=None, owner=None, index=None, name=None, - distribution=None, model=None, transform=None): + distribution=None, model=None, transform=None, + total_size=None): """ Parameters ---------- type : theano type (optional) owner : theano owner (optional) - name : str distribution : Distribution - model : Model""" + model : Model + total_size : scalar Tensor (optional) + needed for upscaling logp + """ if type is None: type = distribution.type super(TransformedRV, self).__init__(type, owner, index, name) @@ -916,7 +956,7 @@ def __init__(self, type=None, owner=None, index=None, name=None, transformed_name = "{}_{}_".format(name, transform.name) self.transformed = model.Var( - transformed_name, transform.apply(distribution)) + transformed_name, transform.apply(distribution), total_size=total_size) normalRV = transform.backward(self.transformed) diff --git a/pymc3/sampling.py b/pymc3/sampling.py index 6c80b49472..95f5660e89 100644 --- a/pymc3/sampling.py +++ b/pymc3/sampling.py @@ -98,16 +98,20 @@ def sample(draws, step=None, init='advi', n_init=200000, start=None, A step function or collection of functions. If no step methods are specified, or are partially specified, they will be assigned automatically (defaults to None). - init : str {'advi', 'advi_map', 'map', 'nuts', None} + init : str {'ADVI', 'ADVI_MAP', 'MAP', 'NUTS', None} Initialization method to use. - * advi : Run ADVI to estimate posterior mean and diagonal covariance matrix. - * advi_map: Initialize ADVI with MAP and use MAP as starting point. - * map : Use the MAP as starting point. - * nuts : Run NUTS and estimate posterior mean and covariance matrix. + * ADVI : Run ADVI to estimate starting points and diagonal covariance + matrix. If njobs > 1 it will sample starting points from the estimated + posterior, otherwise it will use the estimated posterior mean. + * ADVI_MAP: Initialize ADVI with MAP and use MAP as starting point. + * MAP : Use the MAP as starting point. + * NUTS : Run NUTS to estimate starting points and covariance matrix. If + njobs > 1 it will sample starting points from the estimated posterior, + otherwise it will use the estimated posterior mean. * None : Do not initialize. n_init : int Number of iterations of initializer - If 'advi', number of iterations, if 'nuts', number of draws. + If 'ADVI', number of iterations, if 'nuts', number of draws. start : dict Starting point in parameter space (or partial point) Defaults to trace.point(-1)) if there is a trace provided and @@ -142,11 +146,14 @@ def sample(draws, step=None, init='advi', n_init=200000, start=None, MultiTrace object with access to sampling values """ model = modelcontext(model) + + if init is not None: + init = init.lower() if step is None and init is not None and pm.model.all_continuous(model.vars): # By default, use NUTS sampler pm._log.info('Auto-assigning NUTS sampler...') - start_, step = init_nuts(init=init, n_init=n_init, model=model) + start_, step = init_nuts(init=init, njobs=njobs, n_init=n_init, model=model, random_seed=random_seed) if start is None: start = start_ else: @@ -393,7 +400,8 @@ def sample_ppc(trace, samples=None, model=None, vars=None, size=None, random_see return {k: np.asarray(v) for k, v in ppc.items()} -def init_nuts(init='advi', n_init=500000, model=None, **kwargs): +def init_nuts(init='ADVI', njobs=1, n_init=500000, model=None, + random_seed=-1, **kwargs): """Initialize and sample from posterior of a continuous model. This is a convenience function. NUTS convergence and sampling speed is extremely @@ -403,15 +411,17 @@ def init_nuts(init='advi', n_init=500000, model=None, **kwargs): Parameters ---------- - init : str {'advi', 'advi_map', 'map', 'nuts'} + init : str {'ADVI', 'ADVI_MAP', 'MAP', 'NUTS'} Initialization method to use. - * advi : Run ADVI to estimate posterior mean and diagonal covariance matrix. - * advi_map: Initialize ADVI with MAP and use MAP as starting point. - * map : Use the MAP as starting point. - * nuts : Run NUTS and estimate posterior mean and covariance matrix. + * ADVI : Run ADVI to estimate posterior mean and diagonal covariance matrix. + * ADVI_MAP: Initialize ADVI with MAP and use MAP as starting point. + * MAP : Use the MAP as starting point. + * NUTS : Run NUTS and estimate posterior mean and covariance matrix. + njobs : int + Number of parallel jobs to start. n_init : int Number of iterations of initializer - If 'advi', number of iterations, if 'metropolis', number of draws. + If 'ADVI', number of iterations, if 'metropolis', number of draws. model : Model (optional if in `with` context) **kwargs : keyword arguments Extra keyword arguments are forwarded to pymc3.NUTS. @@ -430,23 +440,34 @@ def init_nuts(init='advi', n_init=500000, model=None, **kwargs): pm._log.info('Initializing NUTS using {}...'.format(init)) + random_seed = int(np.atleast_1d(random_seed)[0]) + + if init is not None: + init = init.lower() + if init == 'advi': - v_params = pm.variational.advi(n=n_init) - start = pm.variational.sample_vp(v_params, 1, progressbar=False, hide_transformed=False)[0] + v_params = pm.variational.advi(n=n_init, random_seed=random_seed) + start = pm.variational.sample_vp(v_params, njobs, progressbar=False, + hide_transformed=False, + random_seed=random_seed) + if njobs == 1: + start = start[0] cov = np.power(model.dict_to_array(v_params.stds), 2) elif init == 'advi_map': start = pm.find_MAP() - v_params = pm.variational.advi(n=n_init, start=start) + v_params = pm.variational.advi(n=n_init, start=start, + random_seed=random_seed) cov = np.power(model.dict_to_array(v_params.stds), 2) elif init == 'map': start = pm.find_MAP() cov = pm.find_hessian(point=start) - elif init == 'nuts': - init_trace = pm.sample(step=pm.NUTS(), draws=n_init) - cov = pm.trace_cov(init_trace[n_init//2:]) - - start = {varname: np.mean(init_trace[varname]) for varname in init_trace.varnames} + init_trace = pm.sample(step=pm.NUTS(), draws=n_init, + random_seed=random_seed)[n_init // 2:] + cov = np.atleast_1d(pm.trace_cov(init_trace)) + start = np.random.choice(init_trace, njobs) + if njobs == 1: + start = start[0] else: raise NotImplemented('Initializer {} is not supported.'.format(init)) diff --git a/pymc3/stats.py b/pymc3/stats.py index 5bea76ea8c..7a95622cc3 100644 --- a/pymc3/stats.py +++ b/pymc3/stats.py @@ -105,7 +105,7 @@ def log_post_trace(trace, model): return np.vstack([obs.logp_elemwise(pt) for obs in model.observed_RVs] for pt in trace) -def waic(trace, model=None, n_eff=False): +def waic(trace, model=None, n_eff=False, pointwise=False): """ Calculate the widely available information criterion, its standard error and the effective number of parameters of the samples in trace from model. @@ -121,6 +121,9 @@ def waic(trace, model=None, n_eff=False): n_eff: bool if True the effective number parameters will be returned. Default False + pointwise: bool + if True the pointwise predictive accuracy will be returned. + Default False Returns @@ -152,6 +155,8 @@ def waic(trace, model=None, n_eff=False): if n_eff: return waic, waic_se, p_waic + elif pointwise: + return waic, waic_se, waic_i, p_waic else: return waic, waic_se diff --git a/pymc3/tests/test_diagnostics.py b/pymc3/tests/test_diagnostics.py index d5e6987385..d85738012d 100644 --- a/pymc3/tests/test_diagnostics.py +++ b/pymc3/tests/test_diagnostics.py @@ -20,9 +20,9 @@ def get_ptrace(self, n_samples): # Run sampler step1 = Slice([model.early_mean_log_, model.late_mean_log_]) step2 = Metropolis([model.switchpoint]) - start = {'early_mean': 2., 'late_mean': 3., 'switchpoint': 90} - ptrace = sample(n_samples, [step1, step2], start, njobs=2, progressbar=False, - random_seed=[1, 4]) + start = {'early_mean': 7., 'late_mean': 1., 'switchpoint': 100} + ptrace = sample(n_samples, step=[step1, step2], start=start, njobs=2, + progressbar=False, random_seed=[20090425, 19700903]) return ptrace def test_good(self): @@ -53,7 +53,7 @@ def test_right_shape_python_float(self, shape=None, test_shape=None): # start sampling at the MAP start = find_MAP() step = NUTS(scaling=start) - ptrace = sample(n_samples, step, start, + ptrace = sample(n_samples, step=step, start=start, njobs=n_jobs, random_seed=42) rhat = gelman_rubin(ptrace)['x'] @@ -92,7 +92,7 @@ def get_switchpoint(self, n_samples): # Run sampler step1 = Slice([model.early_mean_log_, model.late_mean_log_]) step2 = Metropolis([model.switchpoint]) - trace = sample(n_samples, [step1, step2], progressbar=False, random_seed=1) + trace = sample(n_samples, step=[step1, step2], progressbar=False, random_seed=1) return trace['switchpoint'] def test_geweke_negative(self): @@ -153,7 +153,7 @@ def test_effective_n(self): # start sampling at the MAP start = find_MAP() step = NUTS(scaling=start) - ptrace = sample(n_samples, step, start, + ptrace = sample(n_samples, step=step, start=start, njobs=n_jobs, random_seed=42) n_effective = effective_n(ptrace)['x'] @@ -174,7 +174,7 @@ def test_effective_n_right_shape_python_float(self, # start sampling at the MAP start = find_MAP() step = NUTS(scaling=start) - ptrace = sample(n_samples, step, start, + ptrace = sample(n_samples, step=step, start=start, njobs=n_jobs, random_seed=42) n_effective = effective_n(ptrace)['x'] diff --git a/pymc3/tests/test_examples.py b/pymc3/tests/test_examples.py index 2ce8f1a793..d071826cce 100644 --- a/pymc3/tests/test_examples.py +++ b/pymc3/tests/test_examples.py @@ -54,7 +54,7 @@ def test_run(self): h = H(start) step = pm.HamiltonianMC(model.vars, h) - pm.sample(50, step, start) + pm.sample(50, step=step, start=start) class TestARM12_6(SeededTest): @@ -88,7 +88,7 @@ def too_slow(self): start = pm.find_MAP(start=start, vars=[model['groupmean'], model['sd_interval_'], model['floor_m']]) step = pm.NUTS(model.vars, scaling=start) - pm.sample(50, step, start) + pm.sample(50, step=step, start=start) class TestARM12_6Uranium(SeededTest): @@ -129,7 +129,7 @@ def too_slow(self): h = np.diag(H(start)) step = pm.HamiltonianMC(model.vars, h) - pm.sample(50, step, start) + pm.sample(50, step=step, start=start) def build_disaster_model(masked=False): @@ -263,7 +263,7 @@ def test_run(self): start = {'psi': 0.5, 'z': (self.y > 0).astype(int), 'theta': 5} step_one = pm.Metropolis([model.theta_interval_, model.psi_logodds_]) step_two = pm.BinaryMetropolis([model.z]) - pm.sample(50, [step_one, step_two], start) + pm.sample(50, step=[step_one, step_two], start=start) class TestRSV(SeededTest): diff --git a/pymc3/tests/test_glm.py b/pymc3/tests/test_glm.py index 259ffdb48b..b4cc107ed7 100644 --- a/pymc3/tests/test_glm.py +++ b/pymc3/tests/test_glm.py @@ -34,7 +34,7 @@ def test_linear_component(self): Normal('y_obs', mu=y_est, sd=sigma, observed=self.y_linear) start = find_MAP(vars=[sigma]) step = Slice(model.vars) - trace = sample(500, step, start, progressbar=False, random_seed=self.random_seed) + trace = sample(500, step=step, start=start, progressbar=False, random_seed=self.random_seed) self.assertAlmostEqual(np.mean(trace['Intercept']), self.intercept, 1) self.assertAlmostEqual(np.mean(trace['x']), self.slope, 1) diff --git a/pymc3/tests/test_model.py b/pymc3/tests/test_model.py index 4eb96e6789..42d56e667e 100644 --- a/pymc3/tests/test_model.py +++ b/pymc3/tests/test_model.py @@ -1,5 +1,6 @@ import unittest -import theano.tensor as tt +from theano import theano, tensor as tt +import numpy as np import pymc3 as pm from pymc3.distributions import HalfCauchy, Normal from pymc3 import Potential, Deterministic @@ -118,3 +119,41 @@ def test_model_root(self): self.assertTrue(model is model.root) with pm.Model() as sub: self.assertTrue(model is sub.root) + + def test_density_scaling(self): + with pm.Model() as model1: + Normal('n', observed=[[1]], total_size=1) + p1 = theano.function([], model1.logpt) + + with pm.Model() as model2: + Normal('n', observed=[[1]], total_size=2) + p2 = theano.function([], model2.logpt) + self.assertEqual(p1() * 2, p2()) + + def test_density_scaling_with_genarator(self): + # We have different size generators + def gen1(): + i = 0 + while True: + yield np.ones((10, 100)) * i + i += 1 + + def gen2(): + i = 0 + while True: + yield np.ones((20, 100)) * i + i += 1 + + # We have same size models + with pm.Model() as model1: + Normal('n', observed=gen1(), total_size=100) + p1 = theano.function([], model1.logpt) + + with pm.Model() as model2: + Normal('n', observed=gen2(), total_size=100) + p2 = theano.function([], model2.logpt) + + # We want densities to be equal + for _ in range(10): + np.testing.assert_almost_equal(p1(), p2()) + # Done diff --git a/pymc3/tests/test_models_linear.py b/pymc3/tests/test_models_linear.py index fac66fc0de..b0a00953e7 100644 --- a/pymc3/tests/test_models_linear.py +++ b/pymc3/tests/test_models_linear.py @@ -44,7 +44,7 @@ def test_linear_component(self): Normal('y_obs', mu=lm.y_est, sd=sigma, observed=self.y_linear) # yields y_obs start = find_MAP(vars=[sigma]) step = Slice(model.vars) - trace = sample(500, step, start, progressbar=False, random_seed=self.random_seed) + trace = sample(500, step=step, start=start, progressbar=False, random_seed=self.random_seed) self.assertAlmostEqual(np.mean(trace['lm_Intercept']), self.intercept, 1) self.assertAlmostEqual(np.mean(trace['lm_x0']), self.slope, 1) @@ -58,7 +58,7 @@ def test_linear_component_from_formula(self): Normal('y_obs', mu=lm.y_est, sd=sigma, observed=self.y_linear) start = find_MAP(vars=[sigma]) step = Slice(model.vars) - trace = sample(500, step, start, progressbar=False, random_seed=self.random_seed) + trace = sample(500, step=step, start=start, progressbar=False, random_seed=self.random_seed) self.assertAlmostEqual(np.mean(trace['Intercept']), self.intercept, 1) self.assertAlmostEqual(np.mean(trace['x']), self.slope, 1) @@ -79,7 +79,7 @@ def test_glm(self): ) start = find_MAP() step = Slice(model.vars) - trace = sample(500, step, start, progressbar=False, random_seed=self.random_seed) + trace = sample(500, step=step, start=start, progressbar=False, random_seed=self.random_seed) self.assertAlmostEqual(np.mean(trace['glm_Intercept']), self.intercept, 1) self.assertAlmostEqual(np.mean(trace['glm_x0']), self.slope, 1) self.assertAlmostEqual(np.mean(trace['glm_sd']), self.sd, 1) @@ -91,7 +91,7 @@ def test_glm_from_formula(self): Glm.from_formula('y ~ x', self.data_linear, name=NAME) start = find_MAP() step = Slice(model.vars) - trace = sample(500, step, start, progressbar=False, random_seed=self.random_seed) + trace = sample(500, step=step, start=start, progressbar=False, random_seed=self.random_seed) self.assertAlmostEqual(np.mean(trace['%s_Intercept' % NAME]), self.intercept, 1) self.assertAlmostEqual(np.mean(trace['%s_x' % NAME]), self.slope, 1) diff --git a/pymc3/tests/test_plots.py b/pymc3/tests/test_plots.py index 05dc8d5a46..a0d23464da 100644 --- a/pymc3/tests/test_plots.py +++ b/pymc3/tests/test_plots.py @@ -20,7 +20,7 @@ def test_plots(): start = model.test_point h = find_hessian(start) step = Metropolis(model.vars, h) - trace = sample(3000, step, start) + trace = sample(3000, step=step, start=start) traceplot(trace) forestplot(trace) @@ -34,7 +34,7 @@ def test_plots_multidimensional(): with model as model: h = np.diag(find_hessian(start)) step = Metropolis(model.vars, h) - trace = sample(3000, step, start) + trace = sample(3000, step=step, start=start) traceplot(trace) plot_posterior(trace) @@ -47,7 +47,7 @@ def test_multichain_plots(): step1 = Slice([model.early_mean_log_, model.late_mean_log_]) step2 = Metropolis([model.switchpoint]) start = {'early_mean': 2., 'late_mean': 3., 'switchpoint': 50} - ptrace = sample(1000, [step1, step2], start, njobs=2) + ptrace = sample(1000, step=[step1, step2], start=start, njobs=2) forestplot(ptrace, varnames=['early_mean', 'late_mean']) autocorrplot(ptrace, varnames=['switchpoint']) diff --git a/pymc3/tests/test_sampling.py b/pymc3/tests/test_sampling.py index 5c81be566c..7d85d2284d 100644 --- a/pymc3/tests/test_sampling.py +++ b/pymc3/tests/test_sampling.py @@ -59,7 +59,7 @@ def test_sample(self): with self.model: for njobs in test_njobs: for steps in [1, 10, 300]: - pm.sample(steps, self.step, {}, None, njobs=njobs, random_seed=self.random_seed) + pm.sample(steps, step=self.step, njobs=njobs, random_seed=self.random_seed) def test_sample_init(self): with self.model: diff --git a/pymc3/tests/test_theanof.py b/pymc3/tests/test_theanof.py new file mode 100644 index 0000000000..ebc7308781 --- /dev/null +++ b/pymc3/tests/test_theanof.py @@ -0,0 +1,50 @@ +import itertools +import unittest +import numpy as np +import theano +from ..theanof import DataGenerator, GeneratorOp, generator + + +def integers(): + i = 0 + while True: + yield np.float32(i) + i += 1 + + +def integers_ndim(ndim): + i = 0 + while True: + yield np.ones((2,) * ndim) * i + i += 1 + + +class TestGenerator(unittest.TestCase): + def test_basic(self): + generator = DataGenerator(integers()) + gop = GeneratorOp(generator)() + self.assertEqual(gop.tag.test_value, np.float32(0)) + f = theano.function([], gop) + self.assertEqual(f(), np.float32(0)) + self.assertEqual(f(), np.float32(1)) + for i in range(2, 100): + f() + self.assertEqual(f(), np.float32(100)) + + def test_ndim(self): + for ndim in range(10): + res = list(itertools.islice(integers_ndim(ndim), 0, 2)) + generator = DataGenerator(integers_ndim(ndim)) + gop = GeneratorOp(generator)() + f = theano.function([], gop) + self.assertEqual(ndim, res[0].ndim) + np.testing.assert_equal(f(), res[0]) + np.testing.assert_equal(f(), res[1]) + + def test_cloning_available(self): + gop = generator(integers()) + res = gop ** 2 + shared = theano.shared(np.float32(10)) + res1 = theano.clone(res, {gop: shared}) + f = theano.function([], res1) + self.assertEqual(f(), np.float32(100)) diff --git a/pymc3/theanof.py b/pymc3/theanof.py index 56bad3fb0e..6870119f7e 100644 --- a/pymc3/theanof.py +++ b/pymc3/theanof.py @@ -3,13 +3,17 @@ from .vartypes import typefilter, continuous_types from theano import theano, scalar, tensor as tt from theano.gof.graph import inputs +from theano.tensor import TensorVariable +from theano.gof import Op +from theano.configparser import change_flags from .memoize import memoize from .blocking import ArrayOrdering +from .data import DataGenerator __all__ = ['gradient', 'hessian', 'hessian_diag', 'inputvars', 'cont_inputs', 'floatX', 'jacobian', 'CallableTensor', 'join_nonshared_inputs', - 'make_shared_replacements'] + 'make_shared_replacements', 'generator'] def inputvars(a): @@ -242,3 +246,46 @@ def __call__(self, input): scalar_identity = IdentityOp(scalar.upgrade_to_float, name='scalar_identity') identity = tt.Elemwise(scalar_identity, name='identity') + + +class GeneratorOp(Op): + """ + Generaror Op is designed for storing python generators + inside theano graph. The main limitation is generator itself. + + There are some important cases when generator becomes exhausted + - not endless generator is passed + - exception is raised while `generator.__next__` is performed + Note: it is dangerous in simple python generators, but ok in + custom class based generators with explicit state + - data type on each iteration should be the same + """ + __props__ = ('generator',) + + def __init__(self, gen): + if not isinstance(gen, DataGenerator): + gen = DataGenerator(gen) + super(GeneratorOp, self).__init__() + self.generator = gen + self.itypes = [] + self.otypes = [self.generator.tensortype] + + def perform(self, node, inputs, output_storage, params=None): + try: + output_storage[0][0] = next(self.generator) + except StopIteration: + output_storage[0][0] = np.nan + + def do_constant_folding(self, node): + return False + + @change_flags(compute_test_value='off') + def __call__(self, *args, **kwargs): + rval = super(GeneratorOp, self).__call__(*args, **kwargs) + rval.tag.test_value = self.generator.test_value + return rval + + +def generator(gen): + """shortcut for `GeneratorOp`""" + return GeneratorOp(gen)()