From 5ec2041b6afd03638d4e0c0238ff340e7f479092 Mon Sep 17 00:00:00 2001 From: Ricardo Vieira Date: Fri, 21 Oct 2022 15:45:13 +0200 Subject: [PATCH 1/4] Remove unnecessary call to `inputvars` in `SMC_KERNEL` --- pymc/smc/kernels.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/pymc/smc/kernels.py b/pymc/smc/kernels.py index 43d060da50..dea2606699 100644 --- a/pymc/smc/kernels.py +++ b/pymc/smc/kernels.py @@ -28,7 +28,6 @@ from pymc.aesaraf import ( compile_pymc, floatX, - inputvars, join_nonshared_inputs, make_shared_replacements, ) @@ -160,7 +159,7 @@ def __init__( self.rng = np.random.default_rng(seed=random_seed) self.model = modelcontext(model) - self.variables = inputvars(self.model.value_vars) + self.variables = self.model.value_vars self.var_info = {} self.tempered_posterior = None From abbc5b8a235f803599a21983b24346aec7d51aae Mon Sep 17 00:00:00 2001 From: Ricardo Vieira Date: Thu, 20 Oct 2022 16:14:27 +0200 Subject: [PATCH 2/4] Implement `get_value_vars_from_user_vars` util to check sampler vars are valid Deprecates old `allinmodel` --- pymc/model.py | 2 ++ pymc/step_methods/hmc/base_hmc.py | 4 +-- pymc/step_methods/metropolis.py | 20 +++++++-------- pymc/step_methods/slicer.py | 5 ++-- pymc/tests/test_util.py | 30 ++++++++++++++++++++++ pymc/tests/tuning/test_starting.py | 27 +------------------- pymc/tuning/starting.py | 14 ++--------- pymc/util.py | 40 ++++++++++++++++++++++++++++++ 8 files changed, 88 insertions(+), 54 deletions(-) diff --git a/pymc/model.py b/pymc/model.py index 2c640b1a69..248c7be462 100644 --- a/pymc/model.py +++ b/pymc/model.py @@ -69,6 +69,7 @@ UNSET, WithMemoization, get_transformed_name, + get_value_vars_from_user_vars, get_var_name, treedict, treelist, @@ -617,6 +618,7 @@ def logp_dlogp_function(self, grad_vars=None, tempered=False, **kwargs): if grad_vars is None: grad_vars = self.continuous_value_vars else: + grad_vars = get_value_vars_from_user_vars(grad_vars, self) for i, var in enumerate(grad_vars): if var.dtype not in continuous_types: raise ValueError(f"Can only compute the gradient of continuous types: {var}") diff --git a/pymc/step_methods/hmc/base_hmc.py b/pymc/step_methods/hmc/base_hmc.py index a08dd86071..6ccf063b68 100644 --- a/pymc/step_methods/hmc/base_hmc.py +++ b/pymc/step_methods/hmc/base_hmc.py @@ -31,6 +31,7 @@ from pymc.step_methods.hmc import integration from pymc.step_methods.hmc.quadpotential import QuadPotentialDiagAdapt, quad_potential from pymc.tuning import guess_scaling +from pymc.util import get_value_vars_from_user_vars logger = logging.getLogger("pymc") @@ -91,8 +92,7 @@ def __init__( if vars is None: vars = self._model.continuous_value_vars else: - vars = [self._model.rvs_to_values.get(var, var) for var in vars] - + vars = get_value_vars_from_user_vars(vars, self._model) super().__init__(vars, blocked=blocked, model=self._model, dtype=dtype, **aesara_kwargs) self.adapt_step_size = adapt_step_size diff --git a/pymc/step_methods/metropolis.py b/pymc/step_methods/metropolis.py index 00b9ecd752..62f556251b 100644 --- a/pymc/step_methods/metropolis.py +++ b/pymc/step_methods/metropolis.py @@ -56,6 +56,8 @@ "MultivariateNormalProposal", ] +from pymc.util import get_value_vars_from_user_vars + # Available proposal distributions for Metropolis @@ -176,9 +178,7 @@ def __init__( if vars is None: vars = model.value_vars else: - vars = [model.rvs_to_values.get(var, var) for var in vars] - - vars = pm.inputvars(vars) + vars = get_value_vars_from_user_vars(vars, model) initial_values_shape = [initial_values[v.name].shape for v in vars] if S is None: @@ -394,7 +394,7 @@ def __init__(self, vars, scaling=1.0, tune=True, tune_interval=100, model=None): self.steps_until_tune = tune_interval self.accepted = 0 - vars = [model.rvs_to_values.get(var, var) for var in vars] + vars = get_value_vars_from_user_vars(vars, model) if not all([v.dtype in pm.discrete_types for v in vars]): raise ValueError("All variables must be Bernoulli for BinaryMetropolis") @@ -484,8 +484,9 @@ def __init__(self, vars, order="random", transit_p=0.8, model=None): # transition probabilities self.transit_p = transit_p + vars = get_value_vars_from_user_vars(vars, model) + initial_point = model.initial_point() - vars = [model.rvs_to_values.get(var, var) for var in vars] self.dim = sum(initial_point[v.name].size for v in vars) if order == "random": @@ -566,8 +567,7 @@ def __init__(self, vars, proposal="uniform", order="random", model=None): model = pm.modelcontext(model) - vars = [model.rvs_to_values.get(var, var) for var in vars] - vars = pm.inputvars(vars) + vars = get_value_vars_from_user_vars(vars, model) initial_point = model.initial_point() @@ -777,8 +777,7 @@ def __init__( if vars is None: vars = model.continuous_value_vars else: - vars = [model.rvs_to_values.get(var, var) for var in vars] - vars = pm.inputvars(vars) + vars = get_value_vars_from_user_vars(vars, model) if S is None: S = np.ones(initial_values_size) @@ -928,8 +927,7 @@ def __init__( if vars is None: vars = model.continuous_value_vars else: - vars = [model.rvs_to_values.get(var, var) for var in vars] - vars = pm.inputvars(vars) + vars = get_value_vars_from_user_vars(vars, model) if S is None: S = np.ones(initial_values_size) diff --git a/pymc/step_methods/slicer.py b/pymc/step_methods/slicer.py index 09d8935ad4..1344b4166d 100644 --- a/pymc/step_methods/slicer.py +++ b/pymc/step_methods/slicer.py @@ -17,10 +17,10 @@ import numpy as np import numpy.random as nr -from pymc.aesaraf import inputvars from pymc.blocking import RaveledVars from pymc.model import modelcontext from pymc.step_methods.arraystep import ArrayStep, Competence +from pymc.util import get_value_vars_from_user_vars from pymc.vartypes import continuous_types __all__ = ["Slice"] @@ -65,8 +65,7 @@ def __init__(self, vars=None, w=1.0, tune=True, model=None, iter_limit=np.inf, * if vars is None: vars = self.model.continuous_value_vars else: - vars = [self.model.rvs_to_values.get(var, var) for var in vars] - vars = inputvars(vars) + vars = get_value_vars_from_user_vars(vars, self.model) super().__init__(vars, [self.model.compile_logp()], **kwargs) diff --git a/pymc/tests/test_util.py b/pymc/tests/test_util.py index 0e2e2c2a4f..320de7b54d 100644 --- a/pymc/tests/test_util.py +++ b/pymc/tests/test_util.py @@ -28,6 +28,7 @@ _get_seeds_per_chain, dataset_to_point_list, drop_warning_stat, + get_value_vars_from_user_vars, hash_key, hashable, locally_cachedmethod, @@ -236,3 +237,32 @@ def test_get_seeds_per_chain(): with pytest.raises(ValueError, match=re.escape("The `seeds` must be array-like")): _get_seeds_per_chain({1: 1, 2: 2}, 2) + + +def test_get_value_vars_from_user_vars(): + with pm.Model() as model1: + x1 = pm.Normal("x1", mu=0, sigma=1) + y1 = pm.Normal("y1", mu=0, sigma=1) + + x1_value = model1.rvs_to_values[x1] + y1_value = model1.rvs_to_values[y1] + assert get_value_vars_from_user_vars([x1, y1], model1) == [x1_value, y1_value] + assert get_value_vars_from_user_vars([x1], model1) == [x1_value] + # The next line does not wrap the variable in a list on purpose, to test the + # utility function can handle those as promised + assert get_value_vars_from_user_vars(x1_value, model1) == [x1_value] + + with pm.Model() as model2: + x2 = pm.Normal("x2", mu=0, sigma=1) + y2 = pm.Normal("y2", mu=0, sigma=1) + det2 = pm.Deterministic("det2", x2 + y2) + + prefix = "The following variables are not random variables in the model:" + with pytest.raises(ValueError, match=rf"{prefix} \['x2', 'y2'\]"): + get_value_vars_from_user_vars([x2, y2], model1) + with pytest.raises(ValueError, match=rf"{prefix} \['x2'\]"): + get_value_vars_from_user_vars([x2, y1], model1) + with pytest.raises(ValueError, match=rf"{prefix} \['x2'\]"): + get_value_vars_from_user_vars([x2], model1) + with pytest.raises(ValueError, match=rf"{prefix} \['det2'\]"): + get_value_vars_from_user_vars([det2], model2) diff --git a/pymc/tests/tuning/test_starting.py b/pymc/tests/tuning/test_starting.py index 528e1fc386..3afd57c3fb 100644 --- a/pymc/tests/tuning/test_starting.py +++ b/pymc/tests/tuning/test_starting.py @@ -22,7 +22,7 @@ from pymc.tests.checks import close_to from pymc.tests.helpers import select_by_precision from pymc.tests.models import non_normal, simple_arbitrary_det, simple_model -from pymc.tuning import find_MAP, starting +from pymc.tuning import find_MAP @pytest.mark.parametrize("bounded", [False, True]) @@ -147,28 +147,3 @@ def test_find_MAP_issue_4488(): assert not set.difference({"x_missing", "x_missing_log__", "y"}, set(map_estimate.keys())) np.testing.assert_allclose(map_estimate["x_missing"], 0.2, rtol=1e-4, atol=1e-4) np.testing.assert_allclose(map_estimate["y"], [2.0, map_estimate["x_missing"][0] + 1]) - - -def test_allinmodel(): - model1 = pm.Model() - model2 = pm.Model() - with model1: - x1 = pm.Normal("x1", mu=0, sigma=1) - y1 = pm.Normal("y1", mu=0, sigma=1) - with model2: - x2 = pm.Normal("x2", mu=0, sigma=1) - y2 = pm.Normal("y2", mu=0, sigma=1) - - x1 = model1.rvs_to_values[x1] - y1 = model1.rvs_to_values[y1] - x2 = model2.rvs_to_values[x2] - y2 = model2.rvs_to_values[y2] - - starting.allinmodel([x1, y1], model1) - starting.allinmodel([x1], model1) - with pytest.raises(ValueError, match=r"Some variables not in the model: \['x2', 'y2'\]"): - starting.allinmodel([x2, y2], model1) - with pytest.raises(ValueError, match=r"Some variables not in the model: \['x2'\]"): - starting.allinmodel([x2, y1], model1) - with pytest.raises(ValueError, match=r"Some variables not in the model: \['x2'\]"): - starting.allinmodel([x2], model1) diff --git a/pymc/tuning/starting.py b/pymc/tuning/starting.py index 70139778ca..b14d0c30c9 100644 --- a/pymc/tuning/starting.py +++ b/pymc/tuning/starting.py @@ -30,11 +30,10 @@ import pymc as pm -from pymc.aesaraf import inputvars from pymc.blocking import DictToArrayBijection, RaveledVars from pymc.initial_point import make_initial_point_fn from pymc.model import modelcontext -from pymc.util import get_default_varnames, get_var_name +from pymc.util import get_default_varnames, get_value_vars_from_user_vars from pymc.vartypes import discrete_types, typefilter __all__ = ["find_MAP"] @@ -96,11 +95,9 @@ def find_MAP( if not vars: raise ValueError("Model has no unobserved continuous variables.") else: - vars = [model.rvs_to_values.get(var, var) for var in vars] + vars = get_value_vars_from_user_vars(vars, model) - vars = inputvars(vars) disc_vars = list(typefilter(vars, discrete_types)) - allinmodel(vars, model) ipfn = make_initial_point_fn( model=model, jitter_rvs=set(), @@ -182,13 +179,6 @@ def allfinite(x): return np.all(isfinite(x)) -def allinmodel(vars, model): - notin = [v for v in vars if v not in model.value_vars] - if notin: - notin = list(map(get_var_name, notin)) - raise ValueError("Some variables not in the model: " + str(notin)) - - class CostFuncWrapper: def __init__(self, maxeval=5000, progressbar=True, logp_func=None, dlogp_func=None): self.n_eval = 0 diff --git a/pymc/util.py b/pymc/util.py index 2ff465cea9..64b3a376fa 100644 --- a/pymc/util.py +++ b/pymc/util.py @@ -21,6 +21,7 @@ import numpy as np import xarray +from aesara import Variable from aesara.compile import SharedVariable from cachetools import LRUCache, cachedmethod @@ -441,3 +442,42 @@ def _get_unique_seeds_per_chain(integers_fn): ) return random_state + + +def get_value_vars_from_user_vars( + vars: Union[Variable, Sequence[Variable]], model +) -> List[Variable]: + """This function converts user "vars" input into value variables + + More often than not, users will pass random variables, and we will extract the + respective value variables, but we also allow for the input to already be value + variables, in case the function is called internally or by a "super-user" + + Returns + ------- + value_vars: list of TensorVariable + List of model value variables that correspond to the input vars + + Raises + ------ + ValueError: + If any of the provided variables do not correspond to any model value variable + """ + if not isinstance(vars, Sequence): + # Single var was passed + value_vars = [model.rvs_to_values.get(vars, vars)] + else: + value_vars = [model.rvs_to_values.get(var, var) for var in vars] + + # Check that we only have value vars from the model + model_value_vars = model.value_vars + notin = [v for v in value_vars if v not in model_value_vars] + if notin: + notin = list(map(get_var_name, notin)) + # We mention random variables, even though the input may be a wrong value variable + # because most users don't know about that duality + raise ValueError( + "The following variables are not random variables in the model: " + str(notin) + ) + + return value_vars From 059ef47f9ae848d9370c7d60db9276893d84246c Mon Sep 17 00:00:00 2001 From: Ricardo Vieira Date: Thu, 20 Oct 2022 17:58:01 +0200 Subject: [PATCH 3/4] Emit warning when passing non free_RVs to `find_MAP` Also improve docstrings --- pymc/aesaraf.py | 6 ++-- pymc/tests/tuning/test_starting.py | 15 +++++++++ pymc/tuning/starting.py | 49 +++++++++++++++++++++--------- 3 files changed, 52 insertions(+), 18 deletions(-) diff --git a/pymc/aesaraf.py b/pymc/aesaraf.py index a22fa5abc9..a10301fd8f 100644 --- a/pymc/aesaraf.py +++ b/pymc/aesaraf.py @@ -261,11 +261,11 @@ def expand_replace(var): def rvs_to_value_vars( - graphs: Iterable[TensorVariable], + graphs: Iterable[Variable], apply_transforms: bool = True, - initial_replacements: Optional[Dict[TensorVariable, TensorVariable]] = None, + initial_replacements: Optional[Dict[Variable, Variable]] = None, **kwargs, -) -> List[TensorVariable]: +) -> List[Variable]: """Clone and replace random variables in graphs with their value variables. This will *not* recompute test values in the resulting graphs. diff --git a/pymc/tests/tuning/test_starting.py b/pymc/tests/tuning/test_starting.py index 3afd57c3fb..e4902e8ced 100644 --- a/pymc/tests/tuning/test_starting.py +++ b/pymc/tests/tuning/test_starting.py @@ -11,6 +11,8 @@ # WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. # See the License for the specific language governing permissions and # limitations under the License. +import re + import numpy as np import pytest @@ -147,3 +149,16 @@ def test_find_MAP_issue_4488(): assert not set.difference({"x_missing", "x_missing_log__", "y"}, set(map_estimate.keys())) np.testing.assert_allclose(map_estimate["x_missing"], 0.2, rtol=1e-4, atol=1e-4) np.testing.assert_allclose(map_estimate["y"], [2.0, map_estimate["x_missing"][0] + 1]) + + +def test_find_MAP_warning_non_free_RVs(): + with pm.Model() as m: + x = pm.Normal("x") + y = pm.Normal("y") + det = pm.Deterministic("det", x + y) + pm.Normal("z", det, 1e-5, observed=100) + + msg = "Intermediate variables (such as Deterministic or Potential) were passed" + with pytest.warns(UserWarning, match=re.escape(msg)): + r = pm.find_MAP(vars=[det]) + np.testing.assert_allclose([r["x"], r["y"], r["det"]], [50, 50, 100]) diff --git a/pymc/tuning/starting.py b/pymc/tuning/starting.py index b14d0c30c9..f9050e3dcf 100644 --- a/pymc/tuning/starting.py +++ b/pymc/tuning/starting.py @@ -18,12 +18,14 @@ @author: johnsalvatier """ import sys +import warnings -from typing import Optional +from typing import Optional, Sequence import aesara.gradient as tg import numpy as np +from aesara import Variable from fastprogress.fastprogress import ProgressBar, progress_bar from numpy import isfinite from scipy.optimize import minimize @@ -41,7 +43,7 @@ def find_MAP( start=None, - vars=None, + vars: Optional[Sequence[Variable]] = None, method="L-BFGS-B", return_raw=False, include_transformed=True, @@ -61,20 +63,23 @@ def find_MAP( Parameters ---------- start: `dict` of parameter values (Defaults to `model.initial_point`) - vars: list - List of variables to optimize and set to optimum (Defaults to all continuous). - method: string or callable - Optimization algorithm (Defaults to 'L-BFGS-B' unless - discrete variables are specified in `vars`, then - `Powell` which will perform better). For instructions on use of a callable, - refer to SciPy's documentation of `optimize.minimize`. - return_raw: bool - Whether to return the full output of scipy.optimize.minimize (Defaults to `False`) + These values will be fixed and used for any free RandomVariables that are + not being optimized. + vars: list of TensorVariable + List of free RandomVariables to optimize the posterior with respect to. + Defaults to all continuous RVs in a model. The respective value variables + may also be passed instead. + method: string or callable, optional + Optimization algorithm. Defaults to 'L-BFGS-B' unless discrete variables are + specified in `vars`, then `Powell` which will perform better. For instructions + on use of a callable, refer to SciPy's documentation of `optimize.minimize`. + return_raw: bool, optional defaults to False + Whether to return the full output of scipy.optimize.minimize include_transformed: bool, optional defaults to True - Flag for reporting automatically transformed variables in addition - to original variables. + Flag for reporting automatically unconstrained transformed values in addition + to the constrained values progressbar: bool, optional defaults to True - Whether or not to display a progress bar in the command line. + Whether to display a progress bar in the command line. maxeval: int, optional, defaults to 5000 The maximum number of times the posterior distribution is evaluated. model: Model (optional if in `with` context) @@ -95,7 +100,21 @@ def find_MAP( if not vars: raise ValueError("Model has no unobserved continuous variables.") else: - vars = get_value_vars_from_user_vars(vars, model) + try: + vars = get_value_vars_from_user_vars(vars, model) + except ValueError as exc: + # Accomodate case where user passed non-pure RV nodes + vars = pm.inputvars(pm.aesaraf.rvs_to_value_vars(vars)) + if vars: + # Make sure they belong to current model again... + vars = get_value_vars_from_user_vars(vars, model) + warnings.warn( + "Intermediate variables (such as Deterministic or Potential) were passed. " + "find_MAP will optimize the underlying free_RVs instead.", + UserWarning, + ) + else: + raise exc disc_vars = list(typefilter(vars, discrete_types)) ipfn = make_initial_point_fn( From f66f320c796386fd1637556a0c12b7e093e63d33 Mon Sep 17 00:00:00 2001 From: Ricardo Vieira Date: Wed, 9 Nov 2022 12:10:52 +0100 Subject: [PATCH 4/4] Temporarily pin mypy version --- conda-envs/environment-dev.yml | 2 +- conda-envs/environment-test.yml | 2 +- conda-envs/windows-environment-dev.yml | 2 +- conda-envs/windows-environment-test.yml | 2 +- requirements-dev.txt | 2 +- 5 files changed, 5 insertions(+), 5 deletions(-) diff --git a/conda-envs/environment-dev.yml b/conda-envs/environment-dev.yml index d6a1895a34..5ceef53dfd 100644 --- a/conda-envs/environment-dev.yml +++ b/conda-envs/environment-dev.yml @@ -38,5 +38,5 @@ dependencies: - watermark - polyagamma - sphinx-remove-toctrees -- mypy +- mypy=0.982 - types-cachetools diff --git a/conda-envs/environment-test.yml b/conda-envs/environment-test.yml index 2dea663920..19fcd7ecd4 100644 --- a/conda-envs/environment-test.yml +++ b/conda-envs/environment-test.yml @@ -27,5 +27,5 @@ dependencies: - pre-commit>=2.8.0 - pytest-cov>=2.5 - pytest>=3.0 -- mypy +- mypy=0.982 - types-cachetools diff --git a/conda-envs/windows-environment-dev.yml b/conda-envs/windows-environment-dev.yml index 6a7012eea7..c5aea815e5 100644 --- a/conda-envs/windows-environment-dev.yml +++ b/conda-envs/windows-environment-dev.yml @@ -35,5 +35,5 @@ dependencies: - sphinx>=1.5 - watermark - sphinx-remove-toctrees -- mypy +- mypy=0.982 - types-cachetools diff --git a/conda-envs/windows-environment-test.yml b/conda-envs/windows-environment-test.yml index a510c4d4aa..b27c453284 100644 --- a/conda-envs/windows-environment-test.yml +++ b/conda-envs/windows-environment-test.yml @@ -28,5 +28,5 @@ dependencies: - pre-commit>=2.8.0 - pytest-cov>=2.5 - pytest>=3.0 -- mypy +- mypy=0.982 - types-cachetools diff --git a/requirements-dev.txt b/requirements-dev.txt index 9475a128f4..09819e6317 100644 --- a/requirements-dev.txt +++ b/requirements-dev.txt @@ -10,7 +10,7 @@ fastprogress>=0.2.0 h5py>=2.7 ipython>=7.16 jupyter-sphinx -mypy +mypy==0.982 myst-nb numpy>=1.15.0 numpydoc