From 20aa34690805f1906a6e9a524d1382bc5d035fb2 Mon Sep 17 00:00:00 2001 From: Chris Fonnesbeck Date: Thu, 17 Nov 2016 13:43:51 -0500 Subject: [PATCH 01/17] Added mode argument to several step methods and advi to allow mode setting (e.g. nanguardmode) for Theano functions --- pymc3/step_methods/metropolis.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pymc3/step_methods/metropolis.py b/pymc3/step_methods/metropolis.py index c2ab7260cab..4fcbdc7866b 100644 --- a/pymc3/step_methods/metropolis.py +++ b/pymc3/step_methods/metropolis.py @@ -433,6 +433,6 @@ def delta_logp(logp, vars, shared): logp1 = pm.CallableTensor(logp0)(inarray1) - f = theano.function([inarray1, inarray0], logp1 - logp0) + f = theano.function([inarray1, inarray0], logp1 - logp0, mode=self.mode) f.trust_input = True return f From 12c556a13b913dbaf1478ad41f1012ae9fcae84f Mon Sep 17 00:00:00 2001 From: Maxim Kochurov Date: Wed, 11 Jan 2017 19:21:20 +0300 Subject: [PATCH 02/17] Created Generator Op with simple Test --- pymc3/data.py | 34 +++++++++++++++++++++++++++++++++- pymc3/tests/test_theanof.py | 18 ++++++++++++++++++ pymc3/theanof.py | 32 ++++++++++++++++++++++++++++++++ 3 files changed, 83 insertions(+), 1 deletion(-) create mode 100644 pymc3/tests/test_theanof.py diff --git a/pymc3/data.py b/pymc3/data.py index 944ce55458e..bd0a6e57350 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,32 @@ 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, iterable): + if hasattr(iterable, '__len__'): + generator = itertools.cycle(iterable) + else: + generator = iter(iterable) + 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) + + 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/tests/test_theanof.py b/pymc3/tests/test_theanof.py new file mode 100644 index 00000000000..c7fa50ab3e6 --- /dev/null +++ b/pymc3/tests/test_theanof.py @@ -0,0 +1,18 @@ +import unittest +import numpy as np +import theano +from ..theanof import DataGenerator, GeneratorOp + + +class TestGenerator(unittest.TestCase): + def test_basic(self): + def integers(): + i = 0 + while True: + yield np.float32(i) + i += 1 + generator = DataGenerator(integers()) + gop = GeneratorOp(generator)() + f = theano.function([], gop) + self.assertEqual(f(), np.float32(0)) + self.assertEqual(f(), np.float32(1)) diff --git a/pymc3/theanof.py b/pymc3/theanof.py index 56bad3fb0e9..8e8d4c6af81 100644 --- a/pymc3/theanof.py +++ b/pymc3/theanof.py @@ -3,8 +3,12 @@ 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, Apply +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', @@ -242,3 +246,31 @@ def __call__(self, input): scalar_identity = IdentityOp(scalar.upgrade_to_float, name='scalar_identity') identity = tt.Elemwise(scalar_identity, name='identity') + + +class GeneratorOp(Op): + __props__ = ('generator',) + + def __init__(self, generator): + if not isinstance(generator, DataGenerator): + raise ValueError('pymc3 requires special DataGenerator for consistency, ' + 'wrap your generator with it') + super(GeneratorOp, self).__init__() + self.generator = generator + self.itypes = [] + self.otypes = [self.generator.tensortype] + + def perform(self, node, inputs, output_storage, params=None): + output_storage[0][0] = self.generator.__next__() + + 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 + + + From b842b3353f17dd630b8f8803cc5c436292e8c8d5 Mon Sep 17 00:00:00 2001 From: Maxim Kochurov Date: Wed, 11 Jan 2017 19:37:25 +0300 Subject: [PATCH 03/17] added ndim test --- pymc3/tests/test_theanof.py | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/pymc3/tests/test_theanof.py b/pymc3/tests/test_theanof.py index c7fa50ab3e6..31a66a943c7 100644 --- a/pymc3/tests/test_theanof.py +++ b/pymc3/tests/test_theanof.py @@ -1,3 +1,4 @@ +import itertools import unittest import numpy as np import theano @@ -16,3 +17,18 @@ def integers(): f = theano.function([], gop) self.assertEqual(f(), np.float32(0)) self.assertEqual(f(), np.float32(1)) + + def test_ndim(self): + for ndim in range(10): + def integers(): + i = 0 + while True: + yield np.ones((2,) * ndim) * i + i += 1 + res = list(itertools.islice(integers(), 0, 2)) + generator = DataGenerator(integers()) + 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]) From 52b4f3a34096a19b4a02b138c696ba2ecb56acc6 Mon Sep 17 00:00:00 2001 From: Maxim Kochurov Date: Wed, 11 Jan 2017 19:39:13 +0300 Subject: [PATCH 04/17] updated test --- pymc3/tests/test_theanof.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/pymc3/tests/test_theanof.py b/pymc3/tests/test_theanof.py index 31a66a943c7..758b5a31734 100644 --- a/pymc3/tests/test_theanof.py +++ b/pymc3/tests/test_theanof.py @@ -17,6 +17,9 @@ def integers(): 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): From dced530b93673d406bbd04967060e4c1d69c9acf Mon Sep 17 00:00:00 2001 From: Maxim Kochurov Date: Wed, 11 Jan 2017 19:43:11 +0300 Subject: [PATCH 05/17] updated test, added test value check --- pymc3/tests/test_theanof.py | 1 + 1 file changed, 1 insertion(+) diff --git a/pymc3/tests/test_theanof.py b/pymc3/tests/test_theanof.py index 758b5a31734..86c32416a70 100644 --- a/pymc3/tests/test_theanof.py +++ b/pymc3/tests/test_theanof.py @@ -14,6 +14,7 @@ def integers(): i += 1 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)) From 76940e79292c10ba27e4fe57367c37cdaeea000c Mon Sep 17 00:00:00 2001 From: Maxim Kochurov Date: Wed, 11 Jan 2017 19:49:26 +0300 Subject: [PATCH 06/17] added test for replacing generator with shared variable --- pymc3/tests/test_theanof.py | 37 +++++++++++++++++++++++++------------ 1 file changed, 25 insertions(+), 12 deletions(-) diff --git a/pymc3/tests/test_theanof.py b/pymc3/tests/test_theanof.py index 86c32416a70..7907b98da6a 100644 --- a/pymc3/tests/test_theanof.py +++ b/pymc3/tests/test_theanof.py @@ -5,13 +5,22 @@ from ..theanof import DataGenerator, GeneratorOp +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): - def integers(): - i = 0 - while True: - yield np.float32(i) - i += 1 generator = DataGenerator(integers()) gop = GeneratorOp(generator)() self.assertEqual(gop.tag.test_value, np.float32(0)) @@ -24,15 +33,19 @@ def integers(): def test_ndim(self): for ndim in range(10): - def integers(): - i = 0 - while True: - yield np.ones((2,) * ndim) * i - i += 1 - res = list(itertools.islice(integers(), 0, 2)) - generator = DataGenerator(integers()) + 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): + generator = DataGenerator(integers()) + gop = GeneratorOp(generator)() + 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)) From a4a42075032c89f51f7eb40601b3a90a6f7ec717 Mon Sep 17 00:00:00 2001 From: Maxim Kochurov Date: Wed, 11 Jan 2017 20:01:17 +0300 Subject: [PATCH 07/17] added shortcut for generator op --- pymc3/theanof.py | 21 ++++++++++++--------- 1 file changed, 12 insertions(+), 9 deletions(-) diff --git a/pymc3/theanof.py b/pymc3/theanof.py index 8e8d4c6af81..666d66220cc 100644 --- a/pymc3/theanof.py +++ b/pymc3/theanof.py @@ -4,7 +4,7 @@ from theano import theano, scalar, tensor as tt from theano.gof.graph import inputs from theano.tensor import TensorVariable -from theano.gof import Op, Apply +from theano.gof import Op from theano.configparser import change_flags from .memoize import memoize from .blocking import ArrayOrdering @@ -13,7 +13,7 @@ __all__ = ['gradient', 'hessian', 'hessian_diag', 'inputvars', 'cont_inputs', 'floatX', 'jacobian', 'CallableTensor', 'join_nonshared_inputs', - 'make_shared_replacements'] + 'make_shared_replacements', 'generator'] def inputvars(a): @@ -251,17 +251,19 @@ def __call__(self, input): class GeneratorOp(Op): __props__ = ('generator',) - def __init__(self, generator): - if not isinstance(generator, DataGenerator): - raise ValueError('pymc3 requires special DataGenerator for consistency, ' - 'wrap your generator with it') + def __init__(self, gen): + if not isinstance(gen, DataGenerator): + gen = DataGenerator(gen) super(GeneratorOp, self).__init__() - self.generator = generator + self.generator = gen self.itypes = [] self.otypes = [self.generator.tensortype] def perform(self, node, inputs, output_storage, params=None): - output_storage[0][0] = self.generator.__next__() + try: + output_storage[0][0] = self.generator.__next__() + except StopIteration: + output_storage[0][0] = np.nan def do_constant_folding(self, node): return False @@ -273,4 +275,5 @@ def __call__(self, *args, **kwargs): return rval - +def generator(gen): + return GeneratorOp(gen)() From 13d5fb3f88594608baeed6ce6bb16f5de9b7425b Mon Sep 17 00:00:00 2001 From: Maxim Kochurov Date: Wed, 11 Jan 2017 20:07:49 +0300 Subject: [PATCH 08/17] refactored test --- pymc3/tests/test_theanof.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/pymc3/tests/test_theanof.py b/pymc3/tests/test_theanof.py index 7907b98da6a..ebc7308781d 100644 --- a/pymc3/tests/test_theanof.py +++ b/pymc3/tests/test_theanof.py @@ -2,7 +2,7 @@ import unittest import numpy as np import theano -from ..theanof import DataGenerator, GeneratorOp +from ..theanof import DataGenerator, GeneratorOp, generator def integers(): @@ -42,8 +42,7 @@ def test_ndim(self): np.testing.assert_equal(f(), res[1]) def test_cloning_available(self): - generator = DataGenerator(integers()) - gop = GeneratorOp(generator)() + gop = generator(integers()) res = gop ** 2 shared = theano.shared(np.float32(10)) res1 = theano.clone(res, {gop: shared}) From 7104be56b550999424dbabe7fcab7a624156a2b4 Mon Sep 17 00:00:00 2001 From: Maxim Kochurov Date: Wed, 11 Jan 2017 20:27:55 +0300 Subject: [PATCH 09/17] added population kwarg (no tests yet) --- pymc3/data.py | 6 +----- pymc3/distributions/distribution.py | 3 ++- pymc3/model.py | 29 +++++++++++++++++++++-------- 3 files changed, 24 insertions(+), 14 deletions(-) diff --git a/pymc3/data.py b/pymc3/data.py index bd0a6e57350..fcac7377197 100644 --- a/pymc3/data.py +++ b/pymc3/data.py @@ -29,11 +29,7 @@ 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, iterable): - if hasattr(iterable, '__len__'): - generator = itertools.cycle(iterable) - else: - generator = iter(iterable) + def __init__(self, generator): self.test_value = next(generator) self.gen = itertools.chain([self.test_value], generator) self.tensortype = tt.TensorType( diff --git a/pymc3/distributions/distribution.py b/pymc3/distributions/distribution.py index 2ad2243f745..d6ab06d6838 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) + population = kwargs.pop('population', None) dist = cls.dist(*args, **kwargs) - return model.Var(name, dist, data) + return model.Var(name, dist, data, population) else: raise TypeError("Name needs to be a string but got: %s" % name) diff --git a/pymc3/model.py b/pymc3/model.py index 6017484d06e..7c34f6047e2 100644 --- a/pymc3/model.py +++ b/pymc3/model.py @@ -8,7 +8,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 +458,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, population=None): """Create and add (un)observed random variable to the model with an appropriate prior distribution. @@ -491,7 +491,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, population=population) self.observed_RVs.append(var) if var.missing_values: self.free_RVs += var.missing_values @@ -500,7 +500,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, + population=population) self.observed_RVs.append(var) if var.missing_values: self.free_RVs.append(var.missing_values) @@ -780,6 +781,8 @@ def as_tensor(data, name, model, distribution): constant[data.mask.nonzero()], missing_values) dataTensor.missing_values = missing_values return dataTensor + elif hasattr(data, '__next__'): + return generator(data) else: data = tt.as_tensor_variable(data, name=name) data.missing_values = None @@ -792,7 +795,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, population=None): """ Parameters ---------- @@ -814,7 +817,12 @@ 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 population is None: + coef = tt.as_tensor(1) + else: + coef = tt.as_tensor(population) / logp_elemwiset.shape[0] + self.logp_elemwiset = logp_elemwiset * coef self.model = model self.distribution = distribution @@ -835,7 +843,7 @@ class MultiObservedRV(Factor): Potentially partially observed. """ - def __init__(self, name, data, distribution, model): + def __init__(self, name, data, distribution, model, population=None): """ Parameters ---------- @@ -851,7 +859,12 @@ 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 population is None: + coef = tt.as_tensor(1) + else: + coef = tt.as_tensor(population) / logp_elemwiset.shape[0] + self.logp_elemwiset = logp_elemwiset * coef self.model = model self.distribution = distribution From a530e7d2970f0e91b53cd70115651296b7a84830 Mon Sep 17 00:00:00 2001 From: Maxim Kochurov Date: Wed, 11 Jan 2017 21:20:23 +0300 Subject: [PATCH 10/17] added population kwarg for free var(autoencoder case) --- pymc3/model.py | 20 ++++++++++++++------ 1 file changed, 14 insertions(+), 6 deletions(-) diff --git a/pymc3/model.py b/pymc3/model.py index 7c34f6047e2..c437c55e750 100644 --- a/pymc3/model.py +++ b/pymc3/model.py @@ -477,11 +477,13 @@ def Var(self, name, dist, data=None, population=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, + population=population) self.free_RVs.append(var) else: var = TransformedRV(name=name, distribution=dist, model=self, - transform=dist.transform) + transform=dist.transform, + population=population) pm._log.debug('Applied {transform}-transform to {name}' ' and added transformed {orig_name} to model.'.format( transform=dist.transform.name, @@ -718,7 +720,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, population=None): """ Parameters ---------- @@ -737,7 +739,12 @@ 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 population is None: + coef = tt.as_tensor(1) + else: + coef = tt.as_tensor(population) / logp_elemwiset.shape[0] + self.logp_elemwiset = logp_elemwiset * coef self.model = model incorporate_methods(source=distribution, destination=self, @@ -909,7 +916,8 @@ 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, + population=None): """ Parameters ---------- @@ -929,7 +937,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), population=population) normalRV = transform.backward(self.transformed) From babd88c7927b38ab5c6018b817706becffc341eb Mon Sep 17 00:00:00 2001 From: Maxim Kochurov Date: Thu, 12 Jan 2017 13:16:28 +0300 Subject: [PATCH 11/17] Revert "Added mode argument to several step methods and advi to allow mode setting (e.g. nanguardmode) for Theano functions" This reverts commit 6bdb53a069c6d569fadcbc36a93489c0ddb334a3. --- pymc3/step_methods/metropolis.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pymc3/step_methods/metropolis.py b/pymc3/step_methods/metropolis.py index 4fcbdc7866b..c2ab7260cab 100644 --- a/pymc3/step_methods/metropolis.py +++ b/pymc3/step_methods/metropolis.py @@ -433,6 +433,6 @@ def delta_logp(logp, vars, shared): logp1 = pm.CallableTensor(logp0)(inarray1) - f = theano.function([inarray1, inarray0], logp1 - logp0, mode=self.mode) + f = theano.function([inarray1, inarray0], logp1 - logp0) f.trust_input = True return f From 50d2abfeda3009f36ed2ba62b32d91e8cdb4f36a Mon Sep 17 00:00:00 2001 From: Maxim Kochurov Date: Sat, 14 Jan 2017 13:25:48 +0300 Subject: [PATCH 12/17] add docstring to generator Op --- pymc3/theanof.py | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/pymc3/theanof.py b/pymc3/theanof.py index 666d66220cc..eeddacf5ce2 100644 --- a/pymc3/theanof.py +++ b/pymc3/theanof.py @@ -249,6 +249,17 @@ def __call__(self, input): 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): @@ -276,4 +287,5 @@ def __call__(self, *args, **kwargs): def generator(gen): + """shortcut for `GeneratorOp`""" return GeneratorOp(gen)() From 2d6f23874464d20ba379abdc6140e15b018968ad Mon Sep 17 00:00:00 2001 From: Maxim Kochurov Date: Sat, 14 Jan 2017 13:27:54 +0300 Subject: [PATCH 13/17] rename population -> total_size --- pymc3/distributions/distribution.py | 4 ++-- pymc3/model.py | 32 ++++++++++++++--------------- 2 files changed, 18 insertions(+), 18 deletions(-) diff --git a/pymc3/distributions/distribution.py b/pymc3/distributions/distribution.py index d6ab06d6838..cae97407407 100644 --- a/pymc3/distributions/distribution.py +++ b/pymc3/distributions/distribution.py @@ -30,9 +30,9 @@ def __new__(cls, name, *args, **kwargs): if isinstance(name, string_types): data = kwargs.pop('observed', None) - population = kwargs.pop('population', None) + total_size = kwargs.pop('total_size', None) dist = cls.dist(*args, **kwargs) - return model.Var(name, dist, data, population) + 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/model.py b/pymc3/model.py index c437c55e750..46ce029b61e 100644 --- a/pymc3/model.py +++ b/pymc3/model.py @@ -458,7 +458,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, population=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. @@ -478,12 +478,12 @@ def Var(self, name, dist, data=None, population=None): if data is None: if getattr(dist, "transform", None) is None: var = FreeRV(name=name, distribution=dist, model=self, - population=population) + total_size=total_size) self.free_RVs.append(var) else: var = TransformedRV(name=name, distribution=dist, model=self, transform=dist.transform, - population=population) + total_size=total_size) pm._log.debug('Applied {transform}-transform to {name}' ' and added transformed {orig_name} to model.'.format( transform=dist.transform.name, @@ -493,7 +493,7 @@ def Var(self, name, dist, data=None, population=None): return var elif isinstance(data, dict): var = MultiObservedRV(name=name, data=data, distribution=dist, - model=self, population=population) + model=self, total_size=total_size) self.observed_RVs.append(var) if var.missing_values: self.free_RVs += var.missing_values @@ -503,7 +503,7 @@ def Var(self, name, dist, data=None, population=None): else: var = ObservedRV(name=name, data=data, distribution=dist, model=self, - population=population) + total_size=total_size) self.observed_RVs.append(var) if var.missing_values: self.free_RVs.append(var.missing_values) @@ -720,7 +720,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, population=None): + distribution=None, model=None, total_size=None): """ Parameters ---------- @@ -740,10 +740,10 @@ def __init__(self, type=None, owner=None, index=None, name=None, self.tag.test_value = np.ones( distribution.shape, distribution.dtype) * distribution.default() logp_elemwiset = distribution.logp(self) - if population is None: + if total_size is None: coef = tt.as_tensor(1) else: - coef = tt.as_tensor(population) / logp_elemwiset.shape[0] + coef = tt.as_tensor(total_size) / logp_elemwiset.shape[0] self.logp_elemwiset = logp_elemwiset * coef self.model = model @@ -802,7 +802,7 @@ class ObservedRV(Factor, TensorVariable): """ def __init__(self, type=None, owner=None, index=None, name=None, data=None, - distribution=None, model=None, population=None): + distribution=None, model=None, total_size=None): """ Parameters ---------- @@ -825,10 +825,10 @@ def __init__(self, type=None, owner=None, index=None, name=None, data=None, self.missing_values = data.missing_values logp_elemwiset = distribution.logp(data) - if population is None: + if total_size is None: coef = tt.as_tensor(1) else: - coef = tt.as_tensor(population) / logp_elemwiset.shape[0] + coef = tt.as_tensor(total_size) / logp_elemwiset.shape[0] self.logp_elemwiset = logp_elemwiset * coef self.model = model self.distribution = distribution @@ -850,7 +850,7 @@ class MultiObservedRV(Factor): Potentially partially observed. """ - def __init__(self, name, data, distribution, model, population=None): + def __init__(self, name, data, distribution, model, total_size=None): """ Parameters ---------- @@ -867,10 +867,10 @@ def __init__(self, name, data, distribution, model, population=None): self.missing_values = [datum.missing_values for datum in self.data.values() if datum.missing_values is not None] logp_elemwiset = distribution.logp(**self.data) - if population is None: + if total_size is None: coef = tt.as_tensor(1) else: - coef = tt.as_tensor(population) / logp_elemwiset.shape[0] + coef = tt.as_tensor(total_size) / logp_elemwiset.shape[0] self.logp_elemwiset = logp_elemwiset * coef self.model = model self.distribution = distribution @@ -917,7 +917,7 @@ class TransformedRV(TensorVariable): def __init__(self, type=None, owner=None, index=None, name=None, distribution=None, model=None, transform=None, - population=None): + total_size=None): """ Parameters ---------- @@ -937,7 +937,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), population=population) + transformed_name, transform.apply(distribution), total_size=total_size) normalRV = transform.backward(self.transformed) From 9fdf17e91ca22955669521b25f8d43be4b959a88 Mon Sep 17 00:00:00 2001 From: Maxim Kochurov Date: Sat, 14 Jan 2017 13:40:53 +0300 Subject: [PATCH 14/17] update docstrings in model --- pymc3/model.py | 17 ++++++++++++++--- 1 file changed, 14 insertions(+), 3 deletions(-) diff --git a/pymc3/model.py b/pymc3/model.py index 46ce029b61e..bdfaaf0c15a 100644 --- a/pymc3/model.py +++ b/pymc3/model.py @@ -469,6 +469,8 @@ def Var(self, name, dist, data=None, total_size=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 ------- @@ -728,7 +730,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) @@ -811,6 +816,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: @@ -859,6 +866,8 @@ def __init__(self, name, data, distribution, model, total_size=None): 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) @@ -924,10 +933,12 @@ def __init__(self, type=None, owner=None, index=None, name=None, 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) From 0ba59bee3e0fdae265837e710904c62943dace9b Mon Sep 17 00:00:00 2001 From: Maxim Kochurov Date: Sat, 14 Jan 2017 13:41:18 +0300 Subject: [PATCH 15/17] fix typo in `as_tensor` function --- pymc3/model.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pymc3/model.py b/pymc3/model.py index bdfaaf0c15a..56241971fe2 100644 --- a/pymc3/model.py +++ b/pymc3/model.py @@ -772,6 +772,8 @@ def pandas_to_array(data): return data elif isinstance(data, theano.gof.graph.Variable): return data + elif hasattr(data, '__next__'): + return generator(data) else: return np.asarray(data) @@ -793,8 +795,6 @@ def as_tensor(data, name, model, distribution): constant[data.mask.nonzero()], missing_values) dataTensor.missing_values = missing_values return dataTensor - elif hasattr(data, '__next__'): - return generator(data) else: data = tt.as_tensor_variable(data, name=name) data.missing_values = None From 1630a7602bac59c24d749e3db513d29e1dc83054 Mon Sep 17 00:00:00 2001 From: Maxim Kochurov Date: Fri, 20 Jan 2017 01:03:38 +0300 Subject: [PATCH 16/17] add test for density rescaling --- pymc3/model.py | 12 +++++++++--- pymc3/tests/test_model.py | 12 +++++++++++- 2 files changed, 20 insertions(+), 4 deletions(-) diff --git a/pymc3/model.py b/pymc3/model.py index 56241971fe2..46d598b7b3d 100644 --- a/pymc3/model.py +++ b/pymc3/model.py @@ -746,8 +746,10 @@ def __init__(self, type=None, owner=None, index=None, name=None, distribution.shape, distribution.dtype) * distribution.default() logp_elemwiset = distribution.logp(self) if total_size is None: - coef = tt.as_tensor(1) + 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 @@ -833,8 +835,10 @@ def __init__(self, type=None, owner=None, index=None, name=None, data=None, logp_elemwiset = distribution.logp(data) if total_size is None: - coef = tt.as_tensor(1) + 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 @@ -877,8 +881,10 @@ def __init__(self, name, data, distribution, model, total_size=None): if datum.missing_values is not None] logp_elemwiset = distribution.logp(**self.data) if total_size is None: - coef = tt.as_tensor(1) + 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 diff --git a/pymc3/tests/test_model.py b/pymc3/tests/test_model.py index 4eb96e6789f..dd86c293a9f 100644 --- a/pymc3/tests/test_model.py +++ b/pymc3/tests/test_model.py @@ -1,5 +1,5 @@ import unittest -import theano.tensor as tt +from theano import theano, tensor as tt import pymc3 as pm from pymc3.distributions import HalfCauchy, Normal from pymc3 import Potential, Deterministic @@ -118,3 +118,13 @@ 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()) From 91011b5870b9d27886f330118739a1e543a6e552 Mon Sep 17 00:00:00 2001 From: Maxim Kochurov Date: Fri, 20 Jan 2017 01:15:33 +0300 Subject: [PATCH 17/17] add test for density with generators --- pymc3/tests/test_model.py | 31 ++++++++++++++++++++++++++++++- 1 file changed, 30 insertions(+), 1 deletion(-) diff --git a/pymc3/tests/test_model.py b/pymc3/tests/test_model.py index dd86c293a9f..42d56e667e6 100644 --- a/pymc3/tests/test_model.py +++ b/pymc3/tests/test_model.py @@ -1,5 +1,6 @@ import unittest 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 @@ -127,4 +128,32 @@ def test_density_scaling(self): with pm.Model() as model2: Normal('n', observed=[[1]], total_size=2) p2 = theano.function([], model2.logpt) - self.assertEqual(p1() * 2, p2()) + 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