Skip to content

fix autoreparam because dims are no longer static #363

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 12 commits into from
Jul 27, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
22 changes: 14 additions & 8 deletions pymc_experimental/model/transforms/autoreparam.py
Original file line number Diff line number Diff line change
@@ -1,3 +1,4 @@
import logging
from dataclasses import dataclass
from functools import singledispatch
from typing import Dict, List, Optional, Sequence, Tuple, Union
Expand All @@ -8,7 +9,6 @@
import pytensor.tensor as pt
import scipy.special
from pymc.distributions import SymbolicRandomVariable
from pymc.exceptions import NotConstantValueError
from pymc.logprob.transforms import Transform
from pymc.model.fgraph import (
ModelDeterministic,
Expand All @@ -19,10 +19,12 @@
model_from_fgraph,
model_named,
)
from pymc.pytensorf import constant_fold, toposort_replace
from pymc.pytensorf import toposort_replace
from pytensor.graph.basic import Apply, Variable
from pytensor.tensor.random.op import RandomVariable

_log = logging.getLogger("pmx")


@dataclass
class VIP:
Expand Down Expand Up @@ -174,15 +176,19 @@ def vip_reparam_node(
) -> Tuple[ModelDeterministic, ModelNamed]:
if not isinstance(node.op, RandomVariable | SymbolicRandomVariable):
raise TypeError("Op should be RandomVariable type")
Copy link
Member

@ricardoV94 ricardoV94 Jul 26, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Do you want to raise NotImplementedError for op.ndim_supp>0?

Not sure if for those you would want one lambda per shape or size, mentioning because those will always be different for multivariate RVs

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

it will not pass the dispatch anyway, it is not needed

Copy link
Member Author

@ferrine ferrine Jul 26, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

ndim_supp>0, will never be supported I think, but I'm not sure, maybe it will, e.g. Dirichlet, there are some reparameterizations with Gamma https://stats.stackexchange.com/questions/548620/reparameterization-trick-for-the-dirichlet-distribution

Copy link
Member Author

@ferrine ferrine Jul 26, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

one lambda per shape or size

Per size for sure 100%

I believe that if I work with multivariate rv, it will be the only way to have one lambda per size (one per independent draw)

rv = node.default_output()
try:
[rv_shape] = constant_fold([rv.shape])
except NotConstantValueError:
raise ValueError("Size should be static for autoreparametrization.")
# FIXME: This is wrong when size is None
_, size, *_ = node.inputs
eval_size = size.eval(mode="FAST_COMPILE")
if eval_size is not None:
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

eval_size will never be None

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It can
image

Copy link
Member

@ricardoV94 ricardoV94 Jul 23, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Wait why are you evaling size instead of rv.shape?

rv.shape.eval()

That's the one you should use

Copy link
Member Author

@ferrine ferrine Jul 23, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I was afraid I get the dependence on the original RV, ah, I do not seem to do

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If you use eval you can only get a numerical output

Copy link
Member

@ricardoV94 ricardoV94 Jul 23, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Do you need shape of transformed or original? If it is the transformed, you could request model.initial_point before starting the dispatch process?

Alternatively we should be able to use shape inference but may need to implement the infer_shape of the ModelVariable. Then the helper infer_static_shape, possibly followed by an eval, should do the job and not complain about value variables because those shouldn't be in the graph after infer_shape

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Can you point me to which test failed when you tried rv.shape.eval()?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Do you need shape of transformed or original? If it is the transformed, you could request model.initial_point before starting the dispatch process?

This looks like an interesting solution

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This looks like an interesting solution

model.initial_point()

Seems to give a point in the sampling space, I need the untransformed one

rv_shape = tuple(eval_size)
else:
rv_shape = ()
lam_name = f"{name}::lam_logit__"
_log.debug(f"Creating {lam_name} with shape of {rv_shape}")
logit_lam_ = pytensor.shared(
np.zeros(rv_shape),
shape=rv_shape,
name=f"{name}::lam_logit__",
name=lam_name,
)
logit_lam = model_named(logit_lam_, *dims)
lam = pt.sigmoid(logit_lam)
Expand Down
35 changes: 31 additions & 4 deletions tests/model/transforms/test_autoreparam.py
Original file line number Diff line number Diff line change
Expand Up @@ -7,20 +7,21 @@

@pytest.fixture
def model_c():
with pm.Model() as mod:
# TODO: Restructure tests so they check one dist at a time
with pm.Model(coords=dict(a=range(5))) as mod:
m = pm.Normal("m")
s = pm.LogNormal("s")
pm.Normal("g", m, s, shape=5)
pm.Normal("g", m, s, dims="a")
pm.Exponential("e", scale=s, shape=7)
return mod


@pytest.fixture
def model_nc():
with pm.Model() as mod:
with pm.Model(coords=dict(a=range(5))) as mod:
m = pm.Normal("m")
s = pm.LogNormal("s")
pm.Deterministic("g", pm.Normal("z", shape=5) * s + m)
pm.Deterministic("g", pm.Normal("z", dims="a") * s + m)
pm.Deterministic("e", pm.Exponential("z_e", 1, shape=7) * s)
return mod

Expand Down Expand Up @@ -102,3 +103,29 @@ def test_set_truncate(model_c: pm.Model):
vip.truncate_lambda(g=0.2)
np.testing.assert_allclose(vip.get_lambda()["g"], 1)
np.testing.assert_allclose(vip.get_lambda()["m"], 0.9)


@pytest.mark.xfail(reason="FIX shape computation for lambda")
def test_lambda_shape():
with pm.Model(coords=dict(a=[1, 2])) as model:
b1 = pm.Normal("b1", dims="a")
b2 = pm.Normal("b2", shape=2)
b3 = pm.Normal("b3", size=2)
b4 = pm.Normal("b4", np.asarray([1, 2]))
model_v, vip = vip_reparametrize(model, ["b1", "b2", "b3", "b4"])
lams = vip.get_lambda()
for v in ["b1", "b2", "b3", "b4"]:
assert lams[v].shape == (2,), v


@pytest.mark.xfail(reason="FIX shape computation for lambda")
def test_lambda_shape_transformed_1d():
with pm.Model(coords=dict(a=[1, 2])) as model:
b1 = pm.Exponential("b1", 1, dims="a")
b2 = pm.Exponential("b2", 1, shape=2)
b3 = pm.Exponential("b3", 1, size=2)
b4 = pm.Exponential("b4", np.asarray([1, 2]))
model_v, vip = vip_reparametrize(model, ["b1", "b2", "b3", "b4"])
lams = vip.get_lambda()
for v in ["b1", "b2", "b3", "b4"]:
assert lams[v].shape == (2,), v
Loading