Fix different / alternative parametrizations between RandomOps and Logp methods #4548
Conversation
ricardoV94
force-pushed
the
fix_halfcauchy
branch
from
8413dd5
to
105ed4a
Compare
There was a problem hiding this comment.
Looks like you have the right idea, and it should be as simple as that!
We should add a simple test for this, though. For now, we can add a single, simple beta parameter test to the/an existing HalfCauchy test method in test_distributions.
Otherwise, we should consider moving this new HalfCauchyRV and the MultinomialRV to Aesara. These changes shouldn't break anything over there, and it would be better to have a single, good implementation of these distributions.
pymc3/distributions/continuous.py
Outdated
| class HalfCauchyRV(RandomVariable): | ||
| name = "halfcauchy" | ||
| ndim_supp = 0 | ||
| ndims_params = [0] | ||
| dtype = "floatX" | ||
| _print_name = ("C**+", "\\operatorname{C^{+}}") |
There was a problem hiding this comment.
It seems like you should be able to subclass the existing HalfCauchyRV, just like MultinomialRV does. That should prevent the need to duplicate these class-level variables, and preserve our ability to do things like isinstance(..., HalfCauchyRV) (using aesara.tensor.random.basic.HalfCauchyRV) on instances of this new class.
There was a problem hiding this comment.
I tried to do that first but I didn't find a way to deal with the __call__ methods as that invokes super.__call__, which I think is strongly linked to the two parameter logic of the original HalfCauchyRV?
The current test_distributions.py::test_half_cauchy already tests for the beta parameter. In terms of logp everything is fine. I am checking how to refactor the test_distributions_random.py to test the random method there.
Sure, that's also fine. |
Before doing that, take a look at the conversation about |
|
I was thinking this one is still useful: https://github.com/pymc-devs/pymc3/blob/9ce64f5a50498fd8db7f2065d2b2427dfc01bc09/pymc3/tests/test_distributions_random.py#L477 It tests that the pymc3 paramaters and scipy random reference give similar values, and it's not concerned with shapes. I refactored the main method to something like this: def pymc3_random(
dist,
paramdomains,
ref_rand,
valuedomain=Domain([0]),
size=10000,
alpha=0.05,
fails=10,
extra_args=None,
model_args=None,
):
if model_args is None:
model_args = {}
model, param_vars = build_model(dist, valuedomain, paramdomains, extra_args)
model_dist = change_rv_size(model.named_vars["value"], size, expand=True)
pymc_rand = aesara.function([], model_dist)
domains = paramdomains.copy()
for pt in product(domains, n_samples=100):
pt = pm.Point(pt, model=model)
pt.update(model_args)
# Update the shared parameter variables in `param_vars`
for k, v in pt.items():
nv = param_vars.get(k, model.named_vars.get(k))
if nv.name in param_vars:
param_vars[nv.name].set_value(v)
p = alpha
# Allow KS test to fail (i.e., the samples be different)
# a certain number of times. Crude, but necessary.
f = fails
while p <= alpha and f > 0:
s0 = pymc_rand()
s1 = ref_rand(size=size, **pt)
_, p = st.ks_2samp(np.atleast_1d(s0).flatten(), np.atleast_1d(s1).flatten())
f -= 1
assert p > alpha, str(pt) |
|
That approach to testing only makes sense when you've implemented a sampler from scratch; otherwise, you're testing the numerical correctness of the underlying NumPy/SciPy samplers. Worse yet, it might actually be confirming that the NumPy/SciPy samplers are correct according to the same NumPy/SciPy samplers. When it comes to checking the mappings between our parameterizations and NumPy/SciPy's, it's a very costly overkill. |
I agree it is an overkill. What should we do instead? |
Standard unit testing. We just need to test the logic we've implemented (e.g. the conversion/remapping of parameters and the calls to SciPy/NumPy). If there is no conversion/remapping, then we really only need to make sure that a single call to the underlying |
|
I found other issues in the mapping of arguments between pymc3 and aesara random ops in the The following already refactored distributions still have issues: I did not refactor the test for discrete variables. We also need to establish a method for mapping optional parametrizations (such as tau/sigma in the HalfNormal). The changes here were just for unused / different arguments. But maybe both can be done in a single swipe. |
ricardoV94
changed the title
ricardoV94
force-pushed
the
fix_halfcauchy
branch
2 times, most recently
from
77259a3
to
e77bed8
Compare
|
This PR will probably just be discarded in the end but hopefully it will help in figuring out how to refactor the old distributions to V4... |
pymc3/tests/test_distributions.py
Outdated
| @@ -1354,7 +1352,7 @@ def test_fun(value, mu, sigma): | |||
| Rplus, | |||
| {"mu": Rplus, "sigma": Rplus}, | |||
| test_fun, | |||
| decimal=select_by_precision(float64=5, float32=3), | |||
| decimal=select_by_precision(float64=3, float32=3), | |||
There was a problem hiding this comment.
We seem to have lost some (test) precision here by using InverseGamma.get_alpha_beta
|
HalfNormal is now passing both logp/cdf and random tests. On On |
|
Okay so every "already-refactored" continuous distributions seems to be in order except for the @brandonwillard do the failing tests say anything to you? It seems the automatic distribution refactoring is not working for this one... in test_distributions.py::TestMatchesScipy::test_mvnormal[1] FAILED [ 33%]
pymc3/tests/test_distributions.py:1631 (TestMatchesScipy.test_mvnormal[1])
self = <pymc3.tests.test_distributions.TestMatchesScipy object at 0x7f750b257850>
n = 1
@pytest.mark.parametrize("n", [1, 2, 3])
def test_mvnormal(self, n):
> self.check_logp(
MvNormal,
RealMatrix(5, n),
{"mu": Vector(R, n), "tau": PdMatrix(n)},
normal_logpdf_tau,
)
test_distributions.py:1634:
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
test_distributions.py:604: in check_logp
model, param_vars = build_model(pymc3_dist, domain, paramdomains, extra_args)
test_distributions.py:232: in build_model
distfam("value", **param_vars, transform=None)
../distributions/distribution.py:166: in __new__
rv_out = cls.dist(*args, rng=rng, **kwargs)
../distributions/multivariate.py:225: in dist
return super().__init__([mu, cov], **kwargs)
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
self = [mu, dot.0], shape = (), dtype = 'float64'
defaults = ('median', 'mean', 'mode'), args = ()
kwargs = {'rng': RandomStateSharedVariable(<RandomState(MT19937) at 0x7F750B304D40>)}
def __init__(self, shape=(), dtype=None, defaults=("median", "mean", "mode"), *args, **kwargs):
if dtype is None:
dtype = aesara.config.floatX
> super().__init__(shape, dtype, defaults=defaults, *args, **kwargs)
E TypeError: super(type, obj): obj must be an instance or subtype of type |
ricardoV94
force-pushed
the
fix_halfcauchy
branch
from
2fc525b
to
1996ed3
Compare
|
On the discrete distributions this part of the It currently raises the error: > return aet.switch(aet.gt(alpha, 1e10), Poisson.dist(mu).logp(value), negbinom)
E AttributeError: 'TensorVariable' object has no attribute 'logp' |
That's odd, because I just re-enabled tests for those distributions (e.g. in |
Yes, that's just a syntax change. We no longer have |
Changes needed to be made in order to handle all the variations of It looks like those changes might need to be further refined and tested. |
Just to be clear, this is for the random part not the logp. It was just that the aesara randomOp works on |
|
The random test for the Error details /home/ricardo/Documents/Projects/pymc3-venv/bin/python3 /opt/pycharm-2020.2.3/plugins/python/helpers/pycharm/_jb_pytest_runner.py --target test_distributions_random.py::TestScalarParameterSamples.test_categorical_random
Testing started at 18:35 ...
Launching pytest with arguments test_distributions_random.py::TestScalarParameterSamples::test_categorical_random in /home/ricardo/Documents/Projects/pymc3/pymc3/tests
============================= test session starts ==============================
platform linux -- Python 3.8.5, pytest-6.1.2, py-1.9.0, pluggy-0.13.1 -- /home/ricardo/Documents/Projects/pymc3-venv/bin/python3
cachedir: .pytest_cache
rootdir: /home/ricardo/Documents/Projects/pymc3, configfile: setup.cfg
plugins: cov-2.10.1
collecting ... collected 1 item
test_distributions_random.py::TestScalarParameterSamples::test_categorical_random[2] FAILED [100%]
pymc3/tests/test_distributions_random.py:789 (TestScalarParameterSamples.test_categorical_random[2])
self = <aesara.compile.function.types.Function object at 0x7fd8277c8a90>
args = (), kwargs = {}
restore_defaults = <function Function.__call__.<locals>.restore_defaults at 0x7fd8278cff70>
profile = None, t0 = 1616002542.7315042, output_subset = None, i = 0
c = <array([0.001, 0.999])>
def __call__(self, *args, **kwargs):
"""
Evaluates value of a function on given arguments.
Parameters
----------
args : list
List of inputs to the function. All inputs are required, even when
some of them are not necessary to calculate requested subset of
outputs.
kwargs : dict
The function inputs can be passed as keyword argument. For this, use
the name of the input or the input instance as the key.
Keyword argument ``output_subset`` is a list of either indices of the
function's outputs or the keys belonging to the `output_keys` dict
and represent outputs that are requested to be calculated. Regardless
of the presence of ``output_subset``, the updates are always calculated
and processed. To disable the updates, you should use the ``copy``
method with ``delete_updates=True``.
Returns
-------
list
List of outputs on indices/keys from ``output_subset`` or all of them,
if ``output_subset`` is not passed.
"""
def restore_defaults():
for i, (required, refeed, value) in enumerate(self.defaults):
if refeed:
if isinstance(value, Container):
value = value.storage[0]
self[i] = value
profile = self.profile
t0 = time.time()
output_subset = kwargs.pop("output_subset", None)
if output_subset is not None and self.output_keys is not None:
output_subset = [self.output_keys.index(key) for key in output_subset]
# Reinitialize each container's 'provided' counter
if self.trust_input:
i = 0
for arg in args:
s = self.input_storage[i]
s.storage[0] = arg
i += 1
else:
for c in self.input_storage:
c.provided = 0
if len(args) + len(kwargs) > len(self.input_storage):
raise TypeError("Too many parameter passed to aesara function")
# Set positional arguments
i = 0
for arg in args:
# TODO: provide a Param option for skipping the filter if we
# really want speed.
s = self.input_storage[i]
# see this emails for a discuation about None as input
# https://groups.google.com/group/theano-dev/browse_thread/thread/920a5e904e8a8525/4f1b311a28fc27e5
if arg is None:
s.storage[0] = arg
else:
try:
s.storage[0] = s.type.filter(
arg, strict=s.strict, allow_downcast=s.allow_downcast
)
except Exception as e:
function_name = "aesara function"
argument_name = "argument"
if self.name:
function_name += ' with name "' + self.name + '"'
if hasattr(arg, "name") and arg.name:
argument_name += ' with name "' + arg.name + '"'
where = get_variable_trace_string(self.maker.inputs[i].variable)
if len(e.args) == 1:
e.args = (
"Bad input "
+ argument_name
+ " to "
+ function_name
+ f" at index {int(i)} (0-based). {where}"
+ e.args[0],
)
else:
e.args = (
"Bad input "
+ argument_name
+ " to "
+ function_name
+ f" at index {int(i)} (0-based). {where}"
) + e.args
restore_defaults()
raise
s.provided += 1
i += 1
# Set keyword arguments
if kwargs: # for speed, skip the items for empty kwargs
for k, arg in kwargs.items():
self[k] = arg
if (
not self.trust_input
and
# The getattr is only needed for old pickle
getattr(self, "_check_for_aliased_inputs", True)
):
# Collect aliased inputs among the storage space
args_share_memory = []
for i in range(len(self.input_storage)):
i_var = self.maker.inputs[i].variable
i_val = self.input_storage[i].storage[0]
if hasattr(i_var.type, "may_share_memory"):
is_aliased = False
for j in range(len(args_share_memory)):
group_j = zip(
[
self.maker.inputs[k].variable
for k in args_share_memory[j]
],
[
self.input_storage[k].storage[0]
for k in args_share_memory[j]
],
)
if any(
[
(
var.type is i_var.type
and var.type.may_share_memory(val, i_val)
)
for (var, val) in group_j
]
):
is_aliased = True
args_share_memory[j].append(i)
break
if not is_aliased:
args_share_memory.append([i])
# Check for groups of more than one argument that share memory
for group in args_share_memory:
if len(group) > 1:
# copy all but the first
for j in group[1:]:
self.input_storage[j].storage[0] = copy.copy(
self.input_storage[j].storage[0]
)
# Check if inputs are missing, or if inputs were set more than once, or
# if we tried to provide inputs that are supposed to be implicit.
if not self.trust_input:
for c in self.input_storage:
if c.required and not c.provided:
restore_defaults()
raise TypeError(
f"Missing required input: {getattr(self.inv_finder[c], 'variable', self.inv_finder[c])}"
)
if c.provided > 1:
restore_defaults()
raise TypeError(
f"Multiple values for input: {getattr(self.inv_finder[c], 'variable', self.inv_finder[c])}"
)
if c.implicit and c.provided > 0:
restore_defaults()
raise TypeError(
f"Tried to provide value for implicit input: {getattr(self.inv_finder[c], 'variable', self.inv_finder[c])}"
)
# Do the actual work
t0_fn = time.time()
try:
outputs = (
> self.fn()
if output_subset is None
else self.fn(output_subset=output_subset)
)
../../../pymc3-venv/lib/python3.8/site-packages/aesara/compile/function/types.py:974:
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
p = <bound method RandomVariable.perform of <aesara.tensor.random.basic.CategoricalRV object at 0x7fd827c32fd0>>
i = [[RandomState(MT19937) at 0x7FD82755F040], [array([100000])], [array(4)], [array([0.001, 0.999])]]
o = [[RandomState(MT19937) at 0x7FD82755F040], [None]]
n = categorical_rv(RandomStateSharedVariable(<RandomState(MT19937) at 0x7FD82755F040>), TensorConstant{(1,) of 100000}, TensorConstant{4}, p)
def rval(p=p, i=node_input_storage, o=node_output_storage, n=node):
> r = p(n, [x[0] for x in i], o)
../../../pymc3-venv/lib/python3.8/site-packages/aesara/graph/op.py:473:
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
self = <aesara.tensor.random.basic.CategoricalRV object at 0x7fd827c32fd0>
node = categorical_rv(RandomStateSharedVariable(<RandomState(MT19937) at 0x7FD82755F040>), TensorConstant{(1,) of 100000}, TensorConstant{4}, p)
inputs = [RandomState(MT19937) at 0x7FD82755F040, array([100000]), array(4), array([0.001, 0.999])]
outputs = [[RandomState(MT19937) at 0x7FD82755F040], [None]]
def perform(self, node, inputs, outputs):
rng_var_out, smpl_out = outputs
rng, size, dtype, *args = inputs
out_var = node.outputs[1]
# If `size == []`, that means no size is enforced, and NumPy is trusted
# to draw the appropriate number of samples, NumPy uses `size=None` to
# represent that. Otherwise, NumPy expects a tuple.
if np.size(size) == 0:
size = None
else:
size = tuple(size)
# Draw from `rng` if `self.inplace` is `True`, and from a copy of `rng`
# otherwise.
if not self.inplace:
rng = copy(rng)
rng_var_out[0] = rng
> smpl_val = self.rng_fn(rng, *(args + [size]))
../../../pymc3-venv/lib/python3.8/site-packages/aesara/tensor/random/op.py:417:
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
cls = <class 'aesara.tensor.random.basic.CategoricalRV'>
rng = RandomState(MT19937) at 0x7FD82755F040, p = array([0.001, 0.999])
size = (100000,)
@classmethod
def rng_fn(cls, rng, p, size):
if size is None:
size = ()
size = tuple(np.atleast_1d(size))
ind_shape = p.shape[:-1]
if len(size) > 0 and size[-len(ind_shape) :] != ind_shape:
> raise ValueError("Parameters shape and size do not match.")
E ValueError: Parameters shape and size do not match.
../../../pymc3-venv/lib/python3.8/site-packages/aesara/tensor/random/basic.py:375: ValueError
During handling of the above exception, another exception occurred:
self = <pymc3.tests.test_distributions_random.TestScalarParameterSamples object at 0x7fd8277c1d60>
s = 2
@pytest.mark.parametrize("s", [2,])
def test_categorical_random(self, s):
def ref_rand(size, p):
return nr.choice(np.arange(p.shape[0]), p=p, size=size)
> pymc3_random_discrete(pm.Categorical, {"p": Simplex(s)}, ref_rand=ref_rand)
test_distributions_random.py:795:
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
test_distributions_random.py:132: in pymc3_random_discrete
o = pymc_rand()
../../../pymc3-venv/lib/python3.8/site-packages/aesara/compile/function/types.py:987: in __call__
raise_with_op(
../../../pymc3-venv/lib/python3.8/site-packages/aesara/link/utils.py:508: in raise_with_op
raise exc_value.with_traceback(exc_trace)
../../../pymc3-venv/lib/python3.8/site-packages/aesara/compile/function/types.py:974: in __call__
self.fn()
../../../pymc3-venv/lib/python3.8/site-packages/aesara/graph/op.py:473: in rval
r = p(n, [x[0] for x in i], o)
../../../pymc3-venv/lib/python3.8/site-packages/aesara/tensor/random/op.py:417: in perform
smpl_val = self.rng_fn(rng, *(args + [size]))
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
cls = <class 'aesara.tensor.random.basic.CategoricalRV'>
rng = RandomState(MT19937) at 0x7FD82755F040, p = array([0.001, 0.999])
size = (100000,)
@classmethod
def rng_fn(cls, rng, p, size):
if size is None:
size = ()
size = tuple(np.atleast_1d(size))
ind_shape = p.shape[:-1]
if len(size) > 0 and size[-len(ind_shape) :] != ind_shape:
> raise ValueError("Parameters shape and size do not match.")
E ValueError: Parameters shape and size do not match.
E Apply node that caused the error: categorical_rv(RandomStateSharedVariable(<RandomState(MT19937) at 0x7FD82755F040>), TensorConstant{(1,) of 100000}, TensorConstant{4}, p)
E Toposort index: 0
E Inputs types: [RandomStateType, TensorType(int64, (True,)), TensorType(int64, scalar), TensorType(float64, vector)]
E Inputs shapes: ['No shapes', (1,), (), (2,)]
E Inputs strides: ['No strides', (8,), (), (8,)]
E Inputs values: [RandomState(MT19937) at 0x7FD82755F040, array([100000]), array(4), array([0.001, 0.999])]
E Inputs type_num: ['', 7, 7, 12]
E Outputs clients: [[], ['output']]
E
E Backtrace when the node is created(use Aesara flag traceback__limit=N to make it longer):
E File "/home/ricardo/Documents/Projects/pymc3-venv/lib/python3.8/site-packages/pluggy/hooks.py", line 286, in __call__
E return self._hookexec(self, self.get_hookimpls(), kwargs)
E File "/home/ricardo/Documents/Projects/pymc3-venv/lib/python3.8/site-packages/pluggy/manager.py", line 93, in _hookexec
E return self._inner_hookexec(hook, methods, kwargs)
E File "/home/ricardo/Documents/Projects/pymc3-venv/lib/python3.8/site-packages/pluggy/manager.py", line 84, in <lambda>
E self._inner_hookexec = lambda hook, methods, kwargs: hook.multicall(
E File "/home/ricardo/Documents/Projects/pymc3-venv/lib/python3.8/site-packages/pluggy/callers.py", line 187, in _multicall
E res = hook_impl.function(*args)
E File "/home/ricardo/Documents/Projects/pymc3-venv/lib/python3.8/site-packages/_pytest/python.py", line 184, in pytest_pyfunc_call
E result = testfunction(**testargs)
E File "/home/ricardo/Documents/Projects/pymc3/pymc3/tests/test_distributions_random.py", line 795, in test_categorical_random
E pymc3_random_discrete(pm.Categorical, {"p": Simplex(s)}, ref_rand=ref_rand)
E File "/home/ricardo/Documents/Projects/pymc3/pymc3/tests/test_distributions_random.py", line 114, in pymc3_random_discrete
E model_dist = change_rv_size(model.named_vars["value"], size, expand=True)
E File "/home/ricardo/Documents/Projects/pymc3/pymc3/distributions/__init__.py", line 130, in change_rv_size
E new_rv_node = rv_node.op.make_node(rng, new_size, dtype, *dist_params)
E
E Debugprint of the apply node:
E categorical_rv.0 [id A] <RandomStateType> ''
E |RandomStateSharedVariable(<RandomState(MT19937) at 0x7FD82755F040>) [id B] <RandomStateType>
E |TensorConstant{(1,) of 100000} [id C] <TensorType(int64, (True,))>
E |TensorConstant{4} [id D] <TensorType(int64, scalar)>
E |p [id E] <TensorType(float64, vector)>
E categorical_rv.1 [id A] <TensorType(int64, vector)> 'value'
E
E Storage map footprint:
E - p, Shared Input, Shape: (2,), ElemSize: 8 Byte(s), TotalSize: 16 Byte(s)
E - TensorConstant{(1,) of 100000}, Shape: (1,), ElemSize: 8 Byte(s), TotalSize: 8 Byte(s)
E - TensorConstant{4}, Shape: (), ElemSize: 8 Byte(s), TotalSize: 8.0 Byte(s)
E TotalSize: 32.0 Byte(s) 0.000 GB
E TotalSize inputs: 32.0 Byte(s) 0.000 GB
../../../pymc3-venv/lib/python3.8/site-packages/aesara/tensor/random/basic.py:375: ValueError
Assertion failed
Assertion failed
Assertion failed
Assertion failed
=================================== FAILURES ===================================
____________ TestScalarParameterSamples.test_categorical_random[2] _____________
self = <aesara.compile.function.types.Function object at 0x7fd8277c8a90>
args = (), kwargs = {}
restore_defaults = <function Function.__call__.<locals>.restore_defaults at 0x7fd8278cff70>
profile = None, t0 = 1616002542.7315042, output_subset = None, i = 0
c = <array([0.001, 0.999])>
def __call__(self, *args, **kwargs):
"""
Evaluates value of a function on given arguments.
Parameters
----------
args : list
List of inputs to the function. All inputs are required, even when
some of them are not necessary to calculate requested subset of
outputs.
kwargs : dict
The function inputs can be passed as keyword argument. For this, use
the name of the input or the input instance as the key.
Keyword argument ``output_subset`` is a list of either indices of the
function's outputs or the keys belonging to the `output_keys` dict
and represent outputs that are requested to be calculated. Regardless
of the presence of ``output_subset``, the updates are always calculated
and processed. To disable the updates, you should use the ``copy``
method with ``delete_updates=True``.
Returns
-------
list
List of outputs on indices/keys from ``output_subset`` or all of them,
if ``output_subset`` is not passed.
"""
def restore_defaults():
for i, (required, refeed, value) in enumerate(self.defaults):
if refeed:
if isinstance(value, Container):
value = value.storage[0]
self[i] = value
profile = self.profile
t0 = time.time()
output_subset = kwargs.pop("output_subset", None)
if output_subset is not None and self.output_keys is not None:
output_subset = [self.output_keys.index(key) for key in output_subset]
# Reinitialize each container's 'provided' counter
if self.trust_input:
i = 0
for arg in args:
s = self.input_storage[i]
s.storage[0] = arg
i += 1
else:
for c in self.input_storage:
c.provided = 0
if len(args) + len(kwargs) > len(self.input_storage):
raise TypeError("Too many parameter passed to aesara function")
# Set positional arguments
i = 0
for arg in args:
# TODO: provide a Param option for skipping the filter if we
# really want speed.
s = self.input_storage[i]
# see this emails for a discuation about None as input
# https://groups.google.com/group/theano-dev/browse_thread/thread/920a5e904e8a8525/4f1b311a28fc27e5
if arg is None:
s.storage[0] = arg
else:
try:
s.storage[0] = s.type.filter(
arg, strict=s.strict, allow_downcast=s.allow_downcast
)
except Exception as e:
function_name = "aesara function"
argument_name = "argument"
if self.name:
function_name += ' with name "' + self.name + '"'
if hasattr(arg, "name") and arg.name:
argument_name += ' with name "' + arg.name + '"'
where = get_variable_trace_string(self.maker.inputs[i].variable)
if len(e.args) == 1:
e.args = (
"Bad input "
+ argument_name
+ " to "
+ function_name
+ f" at index {int(i)} (0-based). {where}"
+ e.args[0],
)
else:
e.args = (
"Bad input "
+ argument_name
+ " to "
+ function_name
+ f" at index {int(i)} (0-based). {where}"
) + e.args
restore_defaults()
raise
s.provided += 1
i += 1
# Set keyword arguments
if kwargs: # for speed, skip the items for empty kwargs
for k, arg in kwargs.items():
self[k] = arg
if (
not self.trust_input
and
# The getattr is only needed for old pickle
getattr(self, "_check_for_aliased_inputs", True)
):
# Collect aliased inputs among the storage space
args_share_memory = []
for i in range(len(self.input_storage)):
i_var = self.maker.inputs[i].variable
i_val = self.input_storage[i].storage[0]
if hasattr(i_var.type, "may_share_memory"):
is_aliased = False
for j in range(len(args_share_memory)):
group_j = zip(
[
self.maker.inputs[k].variable
for k in args_share_memory[j]
],
[
self.input_storage[k].storage[0]
for k in args_share_memory[j]
],
)
if any(
[
(
var.type is i_var.type
and var.type.may_share_memory(val, i_val)
)
for (var, val) in group_j
]
):
is_aliased = True
args_share_memory[j].append(i)
break
if not is_aliased:
args_share_memory.append([i])
# Check for groups of more than one argument that share memory
for group in args_share_memory:
if len(group) > 1:
# copy all but the first
for j in group[1:]:
self.input_storage[j].storage[0] = copy.copy(
self.input_storage[j].storage[0]
)
# Check if inputs are missing, or if inputs were set more than once, or
# if we tried to provide inputs that are supposed to be implicit.
if not self.trust_input:
for c in self.input_storage:
if c.required and not c.provided:
restore_defaults()
raise TypeError(
f"Missing required input: {getattr(self.inv_finder[c], 'variable', self.inv_finder[c])}"
)
if c.provided > 1:
restore_defaults()
raise TypeError(
f"Multiple values for input: {getattr(self.inv_finder[c], 'variable', self.inv_finder[c])}"
)
if c.implicit and c.provided > 0:
restore_defaults()
raise TypeError(
f"Tried to provide value for implicit input: {getattr(self.inv_finder[c], 'variable', self.inv_finder[c])}"
)
# Do the actual work
t0_fn = time.time()
try:
outputs = (
> self.fn()
if output_subset is None
else self.fn(output_subset=output_subset)
)
../../../pymc3-venv/lib/python3.8/site-packages/aesara/compile/function/types.py:974:
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
p = <bound method RandomVariable.perform of <aesara.tensor.random.basic.CategoricalRV object at 0x7fd827c32fd0>>
i = [[RandomState(MT19937) at 0x7FD82755F040], [array([100000])], [array(4)], [array([0.001, 0.999])]]
o = [[RandomState(MT19937) at 0x7FD82755F040], [None]]
n = categorical_rv(RandomStateSharedVariable(<RandomState(MT19937) at 0x7FD82755F040>), TensorConstant{(1,) of 100000}, TensorConstant{4}, p)
def rval(p=p, i=node_input_storage, o=node_output_storage, n=node):
> r = p(n, [x[0] for x in i], o)
../../../pymc3-venv/lib/python3.8/site-packages/aesara/graph/op.py:473:
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
self = <aesara.tensor.random.basic.CategoricalRV object at 0x7fd827c32fd0>
node = categorical_rv(RandomStateSharedVariable(<RandomState(MT19937) at 0x7FD82755F040>), TensorConstant{(1,) of 100000}, TensorConstant{4}, p)
inputs = [RandomState(MT19937) at 0x7FD82755F040, array([100000]), array(4), array([0.001, 0.999])]
outputs = [[RandomState(MT19937) at 0x7FD82755F040], [None]]
def perform(self, node, inputs, outputs):
rng_var_out, smpl_out = outputs
rng, size, dtype, *args = inputs
out_var = node.outputs[1]
# If `size == []`, that means no size is enforced, and NumPy is trusted
# to draw the appropriate number of samples, NumPy uses `size=None` to
# represent that. Otherwise, NumPy expects a tuple.
if np.size(size) == 0:
size = None
else:
size = tuple(size)
# Draw from `rng` if `self.inplace` is `True`, and from a copy of `rng`
# otherwise.
if not self.inplace:
rng = copy(rng)
rng_var_out[0] = rng
> smpl_val = self.rng_fn(rng, *(args + [size]))
../../../pymc3-venv/lib/python3.8/site-packages/aesara/tensor/random/op.py:417:
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
cls = <class 'aesara.tensor.random.basic.CategoricalRV'>
rng = RandomState(MT19937) at 0x7FD82755F040, p = array([0.001, 0.999])
size = (100000,)
@classmethod
def rng_fn(cls, rng, p, size):
if size is None:
size = ()
size = tuple(np.atleast_1d(size))
ind_shape = p.shape[:-1]
if len(size) > 0 and size[-len(ind_shape) :] != ind_shape:
> raise ValueError("Parameters shape and size do not match.")
E ValueError: Parameters shape and size do not match.
../../../pymc3-venv/lib/python3.8/site-packages/aesara/tensor/random/basic.py:375: ValueError
During handling of the above exception, another exception occurred:
self = <pymc3.tests.test_distributions_random.TestScalarParameterSamples object at 0x7fd8277c1d60>
s = 2
@pytest.mark.parametrize("s", [2,])
def test_categorical_random(self, s):
def ref_rand(size, p):
return nr.choice(np.arange(p.shape[0]), p=p, size=size)
> pymc3_random_discrete(pm.Categorical, {"p": Simplex(s)}, ref_rand=ref_rand)
test_distributions_random.py:795:
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
test_distributions_random.py:132: in pymc3_random_discrete
o = pymc_rand()
../../../pymc3-venv/lib/python3.8/site-packages/aesara/compile/function/types.py:987: in __call__
raise_with_op(
../../../pymc3-venv/lib/python3.8/site-packages/aesara/link/utils.py:508: in raise_with_op
raise exc_value.with_traceback(exc_trace)
../../../pymc3-venv/lib/python3.8/site-packages/aesara/compile/function/types.py:974: in __call__
self.fn()
../../../pymc3-venv/lib/python3.8/site-packages/aesara/graph/op.py:473: in rval
r = p(n, [x[0] for x in i], o)
../../../pymc3-venv/lib/python3.8/site-packages/aesara/tensor/random/op.py:417: in perform
smpl_val = self.rng_fn(rng, *(args + [size]))
_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
cls = <class 'aesara.tensor.random.basic.CategoricalRV'>
rng = RandomState(MT19937) at 0x7FD82755F040, p = array([0.001, 0.999])
size = (100000,)
@classmethod
def rng_fn(cls, rng, p, size):
if size is None:
size = ()
size = tuple(np.atleast_1d(size))
ind_shape = p.shape[:-1]
if len(size) > 0 and size[-len(ind_shape) :] != ind_shape:
> raise ValueError("Parameters shape and size do not match.")
E ValueError: Parameters shape and size do not match.
E Apply node that caused the error: categorical_rv(RandomStateSharedVariable(<RandomState(MT19937) at 0x7FD82755F040>), TensorConstant{(1,) of 100000}, TensorConstant{4}, p)
E Toposort index: 0
E Inputs types: [RandomStateType, TensorType(int64, (True,)), TensorType(int64, scalar), TensorType(float64, vector)]
E Inputs shapes: ['No shapes', (1,), (), (2,)]
E Inputs strides: ['No strides', (8,), (), (8,)]
E Inputs values: [RandomState(MT19937) at 0x7FD82755F040, array([100000]), array(4), array([0.001, 0.999])]
E Inputs type_num: ['', 7, 7, 12]
E Outputs clients: [[], ['output']]
E
E Backtrace when the node is created(use Aesara flag traceback__limit=N to make it longer):
E File "/home/ricardo/Documents/Projects/pymc3-venv/lib/python3.8/site-packages/pluggy/hooks.py", line 286, in __call__
E return self._hookexec(self, self.get_hookimpls(), kwargs)
E File "/home/ricardo/Documents/Projects/pymc3-venv/lib/python3.8/site-packages/pluggy/manager.py", line 93, in _hookexec
E return self._inner_hookexec(hook, methods, kwargs)
E File "/home/ricardo/Documents/Projects/pymc3-venv/lib/python3.8/site-packages/pluggy/manager.py", line 84, in <lambda>
E self._inner_hookexec = lambda hook, methods, kwargs: hook.multicall(
E File "/home/ricardo/Documents/Projects/pymc3-venv/lib/python3.8/site-packages/pluggy/callers.py", line 187, in _multicall
E res = hook_impl.function(*args)
E File "/home/ricardo/Documents/Projects/pymc3-venv/lib/python3.8/site-packages/_pytest/python.py", line 184, in pytest_pyfunc_call
E result = testfunction(**testargs)
E File "/home/ricardo/Documents/Projects/pymc3/pymc3/tests/test_distributions_random.py", line 795, in test_categorical_random
E pymc3_random_discrete(pm.Categorical, {"p": Simplex(s)}, ref_rand=ref_rand)
E File "/home/ricardo/Documents/Projects/pymc3/pymc3/tests/test_distributions_random.py", line 114, in pymc3_random_discrete
E model_dist = change_rv_size(model.named_vars["value"], size, expand=True)
E File "/home/ricardo/Documents/Projects/pymc3/pymc3/distributions/__init__.py", line 130, in change_rv_size
E new_rv_node = rv_node.op.make_node(rng, new_size, dtype, *dist_params)
E
E Debugprint of the apply node:
E categorical_rv.0 [id A] <RandomStateType> ''
E |RandomStateSharedVariable(<RandomState(MT19937) at 0x7FD82755F040>) [id B] <RandomStateType>
E |TensorConstant{(1,) of 100000} [id C] <TensorType(int64, (True,))>
E |TensorConstant{4} [id D] <TensorType(int64, scalar)>
E |p [id E] <TensorType(float64, vector)>
E categorical_rv.1 [id A] <TensorType(int64, vector)> 'value'
E
E Storage map footprint:
E - p, Shared Input, Shape: (2,), ElemSize: 8 Byte(s), TotalSize: 16 Byte(s)
E - TensorConstant{(1,) of 100000}, Shape: (1,), ElemSize: 8 Byte(s), TotalSize: 8 Byte(s)
E - TensorConstant{4}, Shape: (), ElemSize: 8 Byte(s), TotalSize: 8.0 Byte(s)
E TotalSize: 32.0 Byte(s) 0.000 GB
E TotalSize inputs: 32.0 Byte(s) 0.000 GB
../../../pymc3-venv/lib/python3.8/site-packages/aesara/tensor/random/basic.py:375: ValueError
=========================== short test summary info ============================
FAILED test_distributions_random.py::TestScalarParameterSamples::test_categorical_random[2]
============================== 1 failed in 0.52s ===============================
Process finished with exit code 1
Assertion failed
Assertion failed |
pymc3/distributions/continuous.py
Outdated
| exponential, | ||
| gamma, | ||
| halfcauchy, | ||
| halfnormal, |
There was a problem hiding this comment.
I have a hard time believing that we need to recreate all of these distributions. Something's wrong here.
There was a problem hiding this comment.
The tests pass now (which was my goal). We can now rework from here and see what is needed or not. The original code before this PR was certainly not working yet.
There was a problem hiding this comment.
Yes, we're building toward complete conversion of the code. It isn't going to happen in a single PR/commit.
It's very important that we not introduce new design issues just because we want to make things "work"/tests pass. Remember, a big part of these v4 changes involves fixing old design issues that most likely arose under similar circumstances.
pymc3/distributions/continuous.py
Outdated
| @@ -722,6 +722,24 @@ def _distr_parameters_for_repr(self): | |||
| return ["mu", "sigma", "lower", "upper"] | |||
|
|
|||
|
|
|||
| class PyMC3HalfNormalRV(RandomVariable): | |||
There was a problem hiding this comment.
This PyMC3* naming scheme isn't good. The RandomVariables defined in PyMC3 should be generally applicable. Furthermore, we should expect that other libraries will use these as the canonical Aesara RandomVariables—when these variables aren't already in Aesara.
There was a problem hiding this comment.
Agree. This was just to distinguish the Aesara / PyMC variations when they were incompatible / extended
pymc3/distributions/continuous.py
Outdated
| class PyMC3HalfCauchyRV(RandomVariable): | ||
| name = "halfcauchy" | ||
| ndim_supp = 0 | ||
| ndims_params = [0] | ||
| dtype = "floatX" | ||
| _print_name = ("C**+", "\\operatorname{C^{+}}") | ||
|
|
||
| def __call__(self, beta=1.0, size=None, **kwargs): | ||
| return super().__call__(beta, size=size, **kwargs) | ||
|
|
||
| @classmethod | ||
| def rng_fn(cls, rng, beta, size): | ||
| return stats.halfcauchy.rvs(scale=beta, random_state=rng, size=size) | ||
|
|
||
|
|
||
| halfcauchy = PyMC3HalfCauchyRV() |
There was a problem hiding this comment.
Again, creating a new class for a reduced parameter space doesn't make sense unless it's absolutely necessary.
Here, it looks like a default loc of zero will work, and not require a new class.
There was a problem hiding this comment.
The problem is that the original Aesara class expects two parameters (or at least that's what I figured out from trying to just inherit from it). Am I mistaken?
There was a problem hiding this comment.
Then we make the HalfCauchy.dist interface take only the used parameter, scale/beta, and set the unused loc parameter to zero.
There was a problem hiding this comment.
Don't we still need to at least invert the order of the arguments in rng_fn. Sounds like a messier solution, but I'm happy to try
Original Aesara HalfCauchyRV:
class HalfCauchyRV(RandomVariable):
name = "cauchy"
ndim_supp = 0
ndims_params = [0, 0]
dtype = "floatX"
_print_name = ("C**+", "\\operatorname{C^{+}}")
def __call__(self, loc=0.0, scale=1.0, size=None, **kwargs):
return super().__call__(loc, scale, size=size, **kwargs)
@classmethod
def rng_fn(cls, rng, loc, scale, size):
return stats.halfcauchy.rvs(loc=loc, scale=scale, random_state=rng, size=size)There was a problem hiding this comment.
Remember, users are not creating RandomVariables directly; they're using Distribution.dist, so that's where we need to focus.
There was a problem hiding this comment.
I think I just need to se an example where the loc parameter is ignored, set to zero by default and the relevant scale argument is used properly. Before this refactoring the scale argument was mistaken as the loc during random sampling (with the scale defaulting to 1). That's why it failed in the reintroduced test_distributions_random.py test
There was a problem hiding this comment.
Yeah, it sounds like there was a bug underlying that issue. I'll try to push the changes I'm alluding to in a minute.
|
Will need to also refactor NegativeBinomial to be in terms of |
Again, we only need to change |
Yes that's exactly what I had in mind. My assumption was that since most randomOps will be created in PyMC3 anyway, we might as well have all of them work with the old pymc3 parametrizations and just tweak those that we want to recycle from Aesara. If the final split is like 80% brand new RandomOps created in Pymc3 and 20% adapted from Asera with the parameters converted to match them, is it worth it to add the extra boilerplate to the logp/logcdf/dist methods instead of writing those 20% from scratch (inheriting from RandomOp)? Whatever the final solution, I got the hang of how to code the parameter conversions. I am still unsure how to do it for parameter reduction, but I'll try to figure it out. |
pymc3/distributions/continuous.py
Outdated
| """ | ||
| rv_op = exponential | ||
|
|
||
| @classmethod | ||
| def dist(cls, lam, *args, **kwargs): | ||
| lam = aet.as_tensor_variable(floatX(lam)) | ||
| def dist(cls, lam=None, mu=None, *args, **kwargs): |
There was a problem hiding this comment.
Introducing new parametrization mu. This is not strictly needed, but since we have to reparametrize behind the scenes for the exponential Op anway, we might as well add it as an option for the logp / logcdf.
|
Reintroduced xfails for:
There are still temporary hacks in place for:
Still to be done:
After this (and review) I think we will be in a position to quickly refactor the remaining standard distributions. |
ricardoV94
force-pushed
the
fix_halfcauchy
branch
from
2f6c0a9
to
3ca2162
Compare
ricardoV94
force-pushed
the
fix_halfcauchy
branch
2 times, most recently
from
2a76568
to
fbca286
Compare
|
Sometimes getting a failing test related to the https://github.com/pymc-devs/pymc3/runs/2147720067?check_suite_focus=true alpha = array([1., 1., 1.])
x = array([0.83469003, 0.02537253, 0.13993742], dtype=float32)
def _dirichlet_check_input(alpha, x):
x = np.asarray(x)
if x.shape[0] + 1 != alpha.shape[0] and x.shape[0] != alpha.shape[0]:
raise ValueError("Vector 'x' must have either the same number "
"of entries as, or one entry fewer than, "
"parameter vector 'a', but alpha.shape = %s "
"and x.shape = %s." % (alpha.shape, x.shape))
if x.shape[0] != alpha.shape[0]:
xk = np.array([1 - np.sum(x, 0)])
if xk.ndim == 1:
x = np.append(x, xk)
elif xk.ndim == 2:
x = np.vstack((x, xk))
else:
raise ValueError("The input must be one dimensional or a two "
"dimensional matrix containing the entries.")
if np.min(x) < 0:
raise ValueError("Each entry in 'x' must be greater than or equal "
"to zero.")
if np.max(x) > 1:
raise ValueError("Each entry in 'x' must be smaller or equal one.")
# Check x_i > 0 or alpha_i > 1
xeq0 = (x == 0)
alphalt1 = (alpha < 1)
if x.shape != alpha.shape:
alphalt1 = np.repeat(alphalt1, x.shape[-1], axis=-1).reshape(x.shape)
chk = np.logical_and(xeq0, alphalt1)
if np.sum(chk):
raise ValueError("Each entry in 'x' must be greater than zero if its "
"alpha is less than one.")
if (np.abs(np.sum(x, 0) - 1.0) > 10e-10).any():
raise ValueError("The input vector 'x' must lie within the normal "
> "simplex. but np.sum(x, 0) = %s." % np.sum(x, 0))
E ValueError: The input vector 'x' must lie within the normal simplex. but np.sum(x, 0) = 0.99999994. |
ricardoV94
force-pushed
the
fix_halfcauchy
branch
from
fbca286
to
6ba2c19
Compare
ricardoV94
force-pushed
the
fix_halfcauchy
branch
4 times, most recently
from
aebdfcf
to
136dc07
Compare
brandonwillard
force-pushed
the
fix_halfcauchy
branch
from
136dc07
to
1c5bc07
Compare
brandonwillard
force-pushed
the
fix_halfcauchy
branch
from
1c5bc07
to
da457f2
Compare
brandonwillard
force-pushed
the
fix_halfcauchy
branch
from
94c5980
to
49f6d6f
Compare
Updated state of this PR in #4548 (comment)