diff --git a/docs/api/hssm.md b/docs/api/hssm.md index 765c1ab2c..081035036 100644 --- a/docs/api/hssm.md +++ b/docs/api/hssm.md @@ -17,6 +17,9 @@ Use `hssm.HSSM` class to construct an HSSM model. - sample_prior_predictive - vi - find_MAP + - find_MLE + - map + - mle - log_likelihood - graph - plot_predictive diff --git a/docs/api/point_estimate.md b/docs/api/point_estimate.md new file mode 100644 index 000000000..18221166e --- /dev/null +++ b/docs/api/point_estimate.md @@ -0,0 +1,14 @@ +`hssm.PointEstimate` is what [`HSSM.find_MAP`](hssm.md) and +[`HSSM.find_MLE`](hssm.md) return, and what the `model.map` and `model.mle` +properties hold. It is a `dict` subclass, so it can be passed anywhere a plain +point dictionary was accepted before — notably `model.sample(initvals=...)` — +while also carrying the optimizer metadata and the ArviZ-friendly exporters +documented below. See the +[point estimation tutorial](../tutorials/map_mle.ipynb) for worked examples. + +::: hssm.PointEstimate + options: + # The project default suppresses attribute docs, but for this class the + # attributes *are* the API — everything the optimizer reports back lives + # there rather than on a method. + show_docstring_attributes: true diff --git a/docs/api/rl.md b/docs/api/rl.md index 26236aa1c..2280d7df1 100644 --- a/docs/api/rl.md +++ b/docs/api/rl.md @@ -25,6 +25,9 @@ sampling model from a named model string. - sample_prior_predictive - vi - find_MAP + - find_MLE + - map + - mle - log_likelihood - plot_predictive - plot_quantile_probability diff --git a/docs/changelog.md b/docs/changelog.md index 350bf8c76..17072aaa3 100644 --- a/docs/changelog.md +++ b/docs/changelog.md @@ -6,7 +6,9 @@ 2. **`plot_model_cartoon` correctness follow-up** (#1125). Each posterior draw's geometry now derives from **one coherent θ vector** — the across-trials mean by default, or one trial via the new `obs=` parameter (which then conditions every simulated layer on that trial; the observed histogram stays pooled, with a warning). New `random_state=` (int or `np.random.Generator`) makes the whole figure reproducible: which posterior draws are displayed (previously chosen by the unseeded global RNG before any seed applied), every simulator seed, and the trajectories — facets consume successive segments of one stream. Geometry and trajectory simulations now receive a `max_t` horizon derived from the x-limits instead of the simulator's 20 s default (previously ~75% of every boundary polyline lay outside the axes); the noisy RT simulations deliberately keep the long horizon so defective densities are not silently re-normalized by censoring. Trajectories are noisy realizations of the reduced reference θ, so their crossing markers land on the drawn boundary by construction. `n_trajectories`, `xlims`, and `ylims` are promoted to documented parameters (None keeps each renderer's defaults). Bug fix: the >2-choice renderer simulated each draw's RT histogram with `n_samples=1` (silently ignoring `n_reps`) and one seed shared across draws, so its uncertainty bands summarized degenerate, noise-correlated histograms. Deliberate value changes for a fixed `random_state`: a new seeding protocol, trial-mean instead of trial-0 geometry for regression models (intercept-only models unaffected), and reference-θ trajectories. Documented caveat: geometry is nonlinear in θ, so the trial-mean curve is not the mean of per-trial curves and need not sit mid-band — the reduction is a display convention. -3. **`plot_model_cartoon` accepts `hist_height="auto"`** (#1127): fits the tallest RT-histogram curve to 90% of the vertical headroom between the histogram baseline (ribbon-aware — expanding bounds are cleared too) and the upper y-limit, so histograms can never overrun the axes; resolves per facet. Its mirror image, `ylims="auto"`, keeps the raw density scale and grows the frame around the content instead (never below the default limits; the two spellings are mutually exclusive). The defaults remain unchanged. The >2-choice renderer now also draws per-choice drift-uncertainty cones in band mode (graded fills in each accumulator's color; previously slope uncertainty appeared only as per-draw spaghetti in the samples display). +3. **Point estimation: `find_MAP` is now trustworthy, and `find_MLE` is new** (#1102). `find_MAP` previously delegated to `pm.find_MAP` with no starting point, so PyMC rebuilt its own (`t=2.0`, `a=2.0`) and ignored HSSM's processed initial values. For many models — including a plain hierarchical DDM on `cavanagh_theta` in the default float64 configuration — the *gradient* of the log-density is non-finite there while the logp is finite, so PyMC's start check passed, L-BFGS-B aborted at `nit=0`, and the untouched start point was stored as `model.map` with no warning. `find_MAP` now defaults `start` to `model.initvals`, audits the optimizer's result and **never caches a failed estimate** (it warns and returns `None`, or raises under `strict=True`); `sample(initvals="map")` raises instead of silently falling back when the MAP does not converge. It also refreshes RL/aDDM extra fields the way `sample()` does, warns and switches to the gradient-free `"Powell"` method under `hssm.set_floatX("float32")` (where L-BFGS-B stops early *and reports success*), and gains `n_starts`, `strict`, `se`, `seed` and `method` parameters. New `find_MLE()` maximizes the observed log-likelihood only; it raises on hierarchical models (the group-level scale is unidentified without the priors under the non-centered parameterization — use `find_MAP`) and on models with `pm.Potential` terms, with an `allow_unidentified=True` escape hatch. Both return a `PointEstimate` — a `dict` subclass, so existing `initvals=` code keeps working — carrying `.params`, `.success`, `.logp`, `.se`, `.to_dataframe()` and `.to_datatree()` (a 1-chain/1-draw posterior that feeds `az.summary`, `sample_posterior_predictive` and the plotting functions; `sd` is 0 and `ess`/`r_hat` are `NaN` by construction). New `model.mle` property, `hssm.PointEstimate` exported at the top level (with a `.copy()` that keeps the metadata a plain `dict.copy()` would drop), and `model.map`'s guard changed from a falsy check to `is None`. Both methods take `progressbar=`, and both record the optimizer that *actually* ran — a likelihood without a gradient is silently switched to `"Powell"` by SciPy's caller, so the requested method would otherwise be misreported. `sample(sampler="laplace")` now warns: bambi's laplace path provides no initval plumbing, so it hits the same bad start and HSSM cannot reach it. See the new "MAP and MLE point estimation" tutorial. + +4. **`plot_model_cartoon` accepts `hist_height="auto"`** (#1127): fits the tallest RT-histogram curve to 90% of the vertical headroom between the histogram baseline (ribbon-aware — expanding bounds are cleared too) and the upper y-limit, so histograms can never overrun the axes; resolves per facet. Its mirror image, `ylims="auto"`, keeps the raw density scale and grows the frame around the content instead (never below the default limits; the two spellings are mutually exclusive). The defaults remain unchanged. The >2-choice renderer now also draws per-choice drift-uncertainty cones in band mode (graded fills in each accumulator's color; previously slope uncertainty appeared only as per-draw spaghetti in the samples display). ### 0.4.0 diff --git a/docs/tutorials/map_mle.ipynb b/docs/tutorials/map_mle.ipynb new file mode 100644 index 000000000..afe4cb9fc --- /dev/null +++ b/docs/tutorials/map_mle.ipynb @@ -0,0 +1,990 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# MAP and MLE point estimation\n", + "\n", + "HSSM is a Bayesian package, but you do not always want a full posterior. Two\n", + "point estimates are available:\n", + "\n", + "- **`model.find_MAP()`** — the *maximum a posteriori* estimate: the mode of the\n", + " joint density, likelihood and priors together.\n", + "- **`model.find_MLE()`** — the *maximum likelihood* estimate: the mode of the\n", + " observed likelihood alone, with the priors dropped.\n", + "\n", + "Both are far cheaper than sampling — seconds on the simple models below,\n", + "though a large hierarchy (or `se=True` on one) takes longer — which makes them\n", + "useful for:\n", + "\n", + "- **debugging a model specification** before committing to a sampling run,\n", + "- **cheap estimation** when you do not need the full posterior,\n", + "- **comparing** a frequentist point estimate against the Bayesian one,\n", + "- **initializing the sampler** via `model.sample(initvals=\"map\")`.\n", + "\n", + "A point estimate is not a posterior. It carries no uncertainty of its own —\n", + "`se=True` below adds asymptotic standard errors, but those are a normal\n", + "approximation around the mode, not the real posterior shape. For inference,\n", + "sample." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Setting PyTensor floatX type to float64.\n", + "Setting \"jax_enable_x64\" to True. If this is not intended, please set `jax` to False.\n" + ] + } + ], + "source": [ + "import arviz as az\n", + "\n", + "import hssm\n", + "\n", + "hssm.set_floatX(\"float64\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Point estimation is gradient-based, so it is sensitive to precision. Under\n", + "`hssm.set_floatX(\"float32\")` the gradient noise exceeds the optimizer's default\n", + "tolerances and the optimizer can stop far from the optimum *while reporting\n", + "success*. HSSM warns and switches to a gradient-free method in that case, but\n", + "`float64` (the default) is what you want here." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Model initialized successfully.\n" + ] + }, + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
rtresponse
01.2591331.0
14.1828881.0
23.3357141.0
33.4901311.0
45.707030-1.0
\n", + "
" + ], + "text/plain": [ + " rt response\n", + "0 1.259133 1.0\n", + "1 4.182888 1.0\n", + "2 3.335714 1.0\n", + "3 3.490131 1.0\n", + "4 5.707030 -1.0" + ] + }, + "execution_count": 2, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "dataset = hssm.simulate_data(\n", + " model=\"ddm\", theta=[0.5, 1.5, 0.5, 0.5], size=500, random_state=42\n", + ")\n", + "model = hssm.HSSM(data=dataset, model=\"ddm\", loglik_kind=\"analytical\")\n", + "dataset.head()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Maximum a posteriori" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "MAP estimate (L-BFGS-B, converged, logp=-1027)\n", + " estimate\n", + "parameter \n", + "a 1.459973\n", + "t 0.512108\n", + "z 0.497889\n", + "v 0.481167" + ] + }, + "execution_count": 3, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "estimate = model.find_MAP(progressbar=False)\n", + "estimate" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The returned object is a `PointEstimate`. It is a `dict` subclass — so it can\n", + "go anywhere a point dictionary was accepted before — carrying the raw optimizer\n", + "point (constrained values, transformed value-variable names such as `t_log__`,\n", + "and any deterministics) plus the optimizer metadata:" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "converged: True\n", + "method: L-BFGS-B\n", + "logp: -1027.110353543079\n", + "message: CONVERGENCE: RELATIVE REDUCTION OF F <= FACTR*EPSMCH\n" + ] + }, + { + "data": { + "text/plain": [ + "{'a': array(1.45997307),\n", + " 't': array(0.51210766),\n", + " 'z': array(0.49788874),\n", + " 'v': array(0.48116658)}" + ] + }, + "execution_count": 4, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "print(\"converged:\", estimate.success)\n", + "print(\"method: \", estimate.method)\n", + "print(\"logp: \", estimate.logp)\n", + "print(\"message: \", estimate.message)\n", + "\n", + "estimate.params" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "`.params` holds only the constrained free parameters — the estimates you would\n", + "report. `.to_dataframe()` tidies them, one row per scalar entry.\n", + "\n", + "### Standard errors\n", + "\n", + "`se=True` adds the inverse Hessian at the optimum, computed on the untransformed\n", + "(constrained) scale. What that quantity *is* depends on which estimate you asked\n", + "for: for an MLE these are standard errors in the usual sense (the observed\n", + "information), while for a MAP they are posterior standard deviations under a\n", + "Laplace — normal — approximation around the mode. Same arithmetic, different\n", + "object.\n", + "\n", + "Expect `NaN` when the optimum sits on a parameter bound or a parameter is weakly\n", + "identified; HSSM warns rather than raising in that case. The cost scales with\n", + "the number of parameters — `2n` gradient evaluations plus an `n × n` inverse —\n", + "so on a large hierarchy this step can outlast the optimization itself." + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
estimatese
parameter
a1.4599730.036789
t0.5121080.027102
z0.4978890.018053
v0.4811670.046027
\n", + "
" + ], + "text/plain": [ + " estimate se\n", + "parameter \n", + "a 1.459973 0.036789\n", + "t 0.512108 0.027102\n", + "z 0.497889 0.018053\n", + "v 0.481167 0.046027" + ] + }, + "execution_count": 5, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "model.find_MAP(se=True, progressbar=False).to_dataframe()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Starting points matter\n", + "\n", + "By default `find_MAP` starts from `model.initvals`, HSSM's own processed initial\n", + "values. This is not cosmetic: PyMC's default start (`t=2.0`) puts many models\n", + "where the *gradient* of the log-density is non-finite while the logp is still\n", + "finite, so the start check passes and the optimizer aborts without moving. HSSM\n", + "detects that and refuses to hand back the start point as if it were an\n", + "estimate." + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "{'a': array(1.5), 't': array(0.025), 'z': array(0.5), 'v': array(0.)}" + ] + }, + "execution_count": 6, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "model.initvals" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "If the optimizer does fail, `find_MAP` warns and returns `None`, and **nothing\n", + "is cached** — so a later `model.map` or `sample(initvals=\"map\")` cannot silently\n", + "pick up a bogus point. Pass `strict=True` to turn the warning into an error.\n", + "\n", + "`n_starts` runs the optimization from several jittered starting points and keeps\n", + "the best one, which helps on multimodal or badly conditioned problems:" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "find_MAP: 3 of 3 starts converged.\n", + "3 of 3 starts converged\n", + "best logp: -1027.1103533898117\n" + ] + } + ], + "source": [ + "multi = model.find_MAP(n_starts=3, seed=0, progressbar=False)\n", + "print(f\"{multi.n_converged} of {multi.n_starts} starts converged\")\n", + "print(\"best logp:\", multi.logp)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Maximum likelihood" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
estimate
parameter
a1.460432
t0.511961
z0.497842
v0.481510
\n", + "
" + ], + "text/plain": [ + " estimate\n", + "parameter \n", + "a 1.460432\n", + "t 0.511961\n", + "z 0.497842\n", + "v 0.481510" + ] + }, + "execution_count": 8, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "mle = model.find_MLE(progressbar=False)\n", + "mle.to_dataframe()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "To compare MAP and MLE, both points must be scored on the **same** objective.\n", + "The MAP maximizes the joint density and the MLE the observed likelihood, so\n", + "their `.logp` values are not comparable to each other. Scoring both on the\n", + "observed log-likelihood is the meaningful comparison — and the MLE should win by\n", + "construction:" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "observed logp at MAP: -1023.3322\n", + "observed logp at MLE: -1023.3321\n", + "difference: 9.71e-05\n" + ] + } + ], + "source": [ + "from hssm.optimize import score_point\n", + "\n", + "map_score = score_point(model.pymc_model, estimate, observed_only=True)\n", + "mle_score = score_point(model.pymc_model, mle, observed_only=True)\n", + "\n", + "print(f\"observed logp at MAP: {map_score:.4f}\")\n", + "print(f\"observed logp at MLE: {mle_score:.4f}\")\n", + "print(f\"difference: {mle_score - map_score:.2e}\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The gap is small here because the priors are weak relative to 500 trials. With\n", + "a small dataset or informative priors, the MAP is pulled further toward the\n", + "prior and the gap widens.\n", + "\n", + "### MLE is refused on hierarchical models\n", + "\n", + "Dropping the priors leaves the group-level scale unidentified, in either\n", + "parameterization:\n", + "\n", + "- **Non-centered** (HSSM's default), where a group-level effect is\n", + " `offset * sigma`: the likelihood is invariant under rescaling `offset` against\n", + " `sigma`, so the optimum is a flat ridge.\n", + "- **Centered** (`noncentered=False`): `sigma` does not enter the likelihood at\n", + " all, so its gradient is exactly zero and the optimizer simply leaves it where\n", + " it started.\n", + "\n", + "Either way there is no maximum to find, so `find_MLE` raises rather than\n", + "returning a number off the ridge.\n", + "\n", + "MAP has no such problem — the priors pin the scale down — so use `find_MAP` for\n", + "hierarchical models. (`allow_unidentified=True` overrides the check if you know\n", + "what you are doing, e.g. for penalized-likelihood work; the estimate you get\n", + "back reports the group-level scale, but you should not read it as identified.)" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Model initialized successfully.\n", + "find_MLE is not well posed for hierarchical models. This model has group-specific terms ['1|participant_id']. Dropping the priors leaves the group-level scale unidentified: under the non-centered parameterization (HSSM's default) the likelihood is invariant under rescaling the offsets against sigma, and under the centered parameterization sigma does not enter the likelihood at all. Either way the optimum is a flat ridge rather than a point. Use `find_MAP()` instead, which is well behaved for these models, or pass `allow_unidentified=True` to proceed anyway.\n" + ] + } + ], + "source": [ + "hierarchical = hssm.HSSM(\n", + " data=hssm.load_data(\"cavanagh_theta\"),\n", + " model=\"ddm\",\n", + " loglik_kind=\"analytical\",\n", + " include=[{\"name\": \"v\", \"formula\": \"v ~ 1 + (1|participant_id)\"}],\n", + ")\n", + "\n", + "try:\n", + " hierarchical.find_MLE()\n", + "except ValueError as error:\n", + " print(error)" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
estimate
parameter
a1.081437
t0.345613
z0.477671
v_Intercept0.503523
v_1|participant_id_sigma0.661233
v_1|participant_id_offset[0]-0.561643
\n", + "
" + ], + "text/plain": [ + " estimate\n", + "parameter \n", + "a 1.081437\n", + "t 0.345613\n", + "z 0.477671\n", + "v_Intercept 0.503523\n", + "v_1|participant_id_sigma 0.661233\n", + "v_1|participant_id_offset[0] -0.561643" + ] + }, + "execution_count": 11, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "hierarchical_map = hierarchical.find_MAP(progressbar=False)\n", + "hierarchical_map.to_dataframe().head(6)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Note that a MAP over a hierarchy is the mode of the *joint* posterior over\n", + "population parameters and group-level offsets, not of the marginal posterior of\n", + "the population parameters.\n", + "\n", + "## Using a point estimate downstream\n", + "\n", + "### As a sampler starting point" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "Initializing NUTS using adapt_diag...\n", + "Sequential sampling (1 chains in 1 job)\n", + "NUTS: [a, t, z, v]\n", + "Sampling 1 chain for 100 tune and 100 draw iterations (100 + 100 draws total) took 1 seconds.\n", + "Only one chain was sampled, this makes it impossible to run some convergence checks\n" + ] + }, + { + "data": { + "text/html": [ + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
meansdeti89_lbeti89_ubess_bulkess_tailr_hatmcse_meanmcse_sd
a1.4550.0351.41.58476nan0.00390.0024
v0.480.0520.40.563468nan0.00890.0043
t0.5180.0280.470.564754nan0.00410.0023
z0.5050.0190.470.533839nan0.00310.0022
" + ], + "text/latex": [ + "\\begin{tabular}{llllllllll}\n", + "\\toprule\n", + " & mean & sd & eti89\\_lb & eti89\\_ub & ess\\_bulk & ess\\_tail & r\\_hat & mcse\\_mean & mcse\\_sd \\\\\n", + "\\midrule\n", + "a & 1.455 & 0.035 & 1.4 & 1.5 & 84 & 76 & nan & 0.0039 & 0.0024 \\\\\n", + "v & 0.48 & 0.052 & 0.4 & 0.56 & 34 & 68 & nan & 0.0089 & 0.0043 \\\\\n", + "t & 0.518 & 0.028 & 0.47 & 0.56 & 47 & 54 & nan & 0.0041 & 0.0023 \\\\\n", + "z & 0.505 & 0.019 & 0.47 & 0.53 & 38 & 39 & nan & 0.0031 & 0.0022 \\\\\n", + "\\bottomrule\n", + "\\end{tabular}\n" + ], + "text/plain": [ + " mean sd eti89_lb eti89_ub ess_bulk ess_tail r_hat mcse_mean mcse_sd\n", + "a 1.455 0.035 1.4 1.5 84 76 nan 0.0039 0.0024\n", + "v 0.48 0.052 0.4 0.56 34 68 nan 0.0089 0.0043\n", + "t 0.518 0.028 0.47 0.56 47 54 nan 0.0041 0.0023\n", + "z 0.505 0.019 0.47 0.53 38 39 nan 0.0031 0.0022" + ] + }, + "execution_count": 12, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "model.sample(initvals=\"map\", draws=100, tune=100, chains=1, cores=1, progressbar=False)\n", + "az.summary(model.traces)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### As a one-draw posterior\n", + "\n", + "`to_datatree()` returns the estimate as a 1-chain, 1-draw posterior, which the\n", + "rest of the ArviZ and HSSM toolchain accepts — including posterior predictive\n", + "sampling and the plotting functions.\n", + "\n", + "Read the summary with care: with a single draw there is no sampling\n", + "distribution, so `sd` is `0` and `ess_bulk`, `ess_tail` and `r_hat` come back\n", + "`NaN`. That is expected, not a diagnostic failure." + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
meansdeti89_lbeti89_ubess_bulkess_tailr_hatmcse_meanmcse_sd
a1.501.51.5nannannannannan
t0.5100.510.51nannannannannan
z0.500.50.5nannannannannan
v0.4800.480.48nannannannannan
" + ], + "text/latex": [ + "\\begin{tabular}{llllllllll}\n", + "\\toprule\n", + " & mean & sd & eti89\\_lb & eti89\\_ub & ess\\_bulk & ess\\_tail & r\\_hat & mcse\\_mean & mcse\\_sd \\\\\n", + "\\midrule\n", + "a & 1.5 & 0 & 1.5 & 1.5 & nan & nan & nan & nan & nan \\\\\n", + "t & 0.51 & 0 & 0.51 & 0.51 & nan & nan & nan & nan & nan \\\\\n", + "z & 0.5 & 0 & 0.5 & 0.5 & nan & nan & nan & nan & nan \\\\\n", + "v & 0.48 & 0 & 0.48 & 0.48 & nan & nan & nan & nan & nan \\\\\n", + "\\bottomrule\n", + "\\end{tabular}\n" + ], + "text/plain": [ + " mean sd eti89_lb eti89_ub ess_bulk ess_tail r_hat mcse_mean mcse_sd\n", + "a 1.5 0 1.5 1.5 nan nan nan nan nan\n", + "t 0.51 0 0.51 0.51 nan nan nan nan nan\n", + "z 0.5 0 0.5 0.5 nan nan nan nan nan\n", + "v 0.48 0 0.48 0.48 nan nan nan nan nan" + ] + }, + "execution_count": 13, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "point_fit = estimate.to_datatree()\n", + "az.summary(point_fit)" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "dt=None, we use the traces assigned to the HSSM object as datatree.\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/Users/fmuia/Library/CloudStorage/Dropbox/02-Data-Science/github-projects/Inccbrown/hssm-repo/implement_mle_map/.venv/lib/python3.13/site-packages/pytensor/link/numba/dispatch/basic.py:234: UserWarning: Numba will use object mode to run ddm_RV_rv{\"(),(),(),(),()->(2)\"}'s perform method. Set `pytensor.config.compiler_verbose = True` to see more details.\n", + " warnings.warn(\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "n_samples > n_draws. Using the entire dataset.\n" + ] + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkAAAAHGCAYAAACcmzRuAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjExLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvctoD+AAAAAlwSFlzAAAPYQAAD2EBqD+naQAAWzlJREFUeJzt3QeUFNXWt/FNBslBMiJJEMREFkGCiCCCYAATYgAjKhhRr+EqQRRRREUMKKKiJFEMV+EiIIgkUTEQFCRHyRmmv/U/963+eobpMDM9sZ7fWr2Grq6uPl3dTO3ZZ59zcgUCgYABAAD4SO7MbgAAAEBGIwACAAC+QwAEAAB8hwAIAAD4DgEQAADwHQIgAADgOwRAAADAdwiAAACA7xAAAdnQxIkTbcKECeYHq1evthdffNE2btwY3PbLL7+4bfv27Yvb66THMdNLRrc1uddbtGiR23b48OEMaUO4dgCplYuZoJHTHThwwEaPHh28nytXLitWrJidffbZds4556TLa3700UeWJ08eu+KKK9Ll+K1atbJjx47Zd999Z+nts88+sz///DN4v0CBAnbKKadY69at7aSTTkr31582bZpdeumlNmfOHDv//PPdtlGjRtntt99u69ats8qVK8d8rB9//NFmzZplffr0OaHtqT1mWq1YscK++OKL4P38+fNb8eLF7bTTTrMzzzzTne+k0uP9R5Lc6w0ZMsQGDBhgO3futBIlSsR8rOz4GSFnIgOEHG/Pnj3Wr18/e+edd2zNmjXuYq4MSoMGDaxNmza2e/fuuL/m8OHDbeTIkZZeFFhdddVVlhHeeust69+/vzt3uv3www92yy23WJUqVezLL7+0zKDA4J577rGiRYum6Hm6sOq7oO9EvI6ZVkuWLHFtmj59uju/v//+uwv6evToYeXLl7cHH3zQ9u/fn+7vP5KMPDdZ8TNCzpQ3sxsAZBRlD5Q+93z44Yd2zTXX2GOPPWYvv/xytvog7rrrrgx9vdy5cyc6d7t27bKzzjrLnT/9NV6kSJEMbc95553nbln9mCnRq1evEzKGkyZNcttnz55tM2fOtEKFCmVKWzP73GS1diBnIACCb1199dV222232ddffx3clpCQ4C42f/zxhxUsWNCaNWtmtWvXPuG5y5Yts6VLl9qRI0fs9NNPt6ZNm7quNXnjjTds8+bNljdv3mDQULhwYevdu3eiY/z888+ujuLo0aMumNAxQr3//vvueZdddpn99ttvNm/ePNf1dNFFF7kMltYxvvLKKxM9J5b2RzpurNTl0b17d3vuueds4cKF1qhRI3vzzTddt5j+Sv/vf//r2tC+fXurWbOme87WrVvt22+/te3bt1ulSpWsXbt2yXbDzJ8/353bcuXKWceOHcPWgsyYMcNlopIGXzq+sgjbtm2zGjVq2AUXXOC6lbRNN+8z8rIICjrUnZL0mGqD2nvttdfaySefnOg1/vnnHxs7dqw7dmg3aqzvMVaXX365y1DefPPNLkhXNijS+9+xY0fwveszVbCg7jSJ9v71XVSXqrqYDh486LJ7ep933nlnxPPtnXPtr/8P+g5Ur1490eM6J/qu3XHHHYm2KxurLlYv25XSzyjUypUrbe7cuXbo0CH3nW/ZsqXrhvYkfX+ff/656x5Xd3KtWrVS/Rkh+6ILDL6WL18+V0sj6n5QXZC6lvSLUkXGZ5xxhgtcjh8/7vZR0NGzZ08XWHz66afuYq2LkoIXXXREGREVhuoXsddtpG0e1Ux06NDBHUO/hPVaqnFRsBBa3KngQjUPTz31lN16660uqJg8ebJ7TN1rL730UqL3Ekv7ox03pedOdP6UEVK3xX/+8x/33tS+Tz75xAURXr2ILsi6iKvL5/HHH3eBkbrTPDqOLoS6cOli+vHHH7uL6ZYtW054bV3o9Hp63VCDBg1yXXN6PV3whg0bZueee667OOpirpv3GXmfjT6n5I6pbIvuqwswqXfffdc9psDKE8t7TI3rr7/eXew/+OCDiO9f50uvP2LECHfeX3vtNdfNq+40ifb+1QWnYyrTpPP+1VdfBTOj4c63fP/993bhhRe675G6mVW7pC7gUArYveAtlAIaHVftiKWNybVDQb/+kNEfIjoH2kffI/1f+Ouvv4L7ee9PAb8CU90fM2aMe964ceNS+ekgW1MRNJCTbdq0KaCv+p133plo+8yZM93266+/3t1v0KBBoHLlyoH169cH95k8ebLbZ9CgQe7+rFmz3P2vv/460bEWL14c2LJlS/B+kyZNAhdccEGy7enQoUOgdOnSgVWrVgW3bdiwIVC2bNlAnz59gtvOOuusQPny5QNPPvlkovciOnbz5s0THTeW9kc7bnK6dOkSyJMnT6JtR48eDZxxxhmBggULBrZv3x5Yt26de50qVaoE5s+fn+i47777rnts1KhRwe3Hjh0LXHXVVYEKFSoE9u/f77YNHjzY7ffVV18F91uxYkWgTp06bvucOXOC21977TW3Ta/reeutt9y24cOHJ2qrzvOff/7p/q3HtE9y7ze5Y55//vmB00477YR99d6bNm0avB/re0zOhx9+6J47YcKEsPvotXLnzh04fvx42LZWrFgxcNNNNyV63j///JPo84j0/r3z37lz58CePXvcNm+/5F4vdP+9e/cGtz/yyCOBXLlyBebOnRvcpv97hQsXPuE1p0yZ4o7x/fffx9TG5NoxdOhQt+2jjz4Kbtu4cWOgatWq7ruekJCQqL1XXnllYN++fW6bHrv00kvd/4dDhw4le+6Rc5EBgm/89NNPrktKf52qhqZTp04u9T1w4EBbsGCBLV682P2FqK4LT9euXV3tkFfQrOyNeH+RepRlKFu2bNQ2qOtM2Y0HHnjAdc94Klas6LIxyiyEDitWqj70L2d1EyQn1van9LgeZb507nRT5qhhw4a2fPlyl20oXbp0cD91hTVp0iTRcYcOHer+Gtf786hr4oknnrBNmza5LhDRSD11RygT5tHnc8kll1gslNlSV+K9996baLvOc9IumVip60mjtDQCLfRc63PUY55Y32NqqRtLmY5ww7/1+Sgrou+O/u0pWbJkos8jFqo58rqeon0v5KabbkrUHfXoo4+6+6+//rplBGW69L0LHRRQoUIF939M/+eVEUraXnUBi7qtr7vuOtdlrS5b+As1QPCNvXv3unS6fumVKlXKBRsKgjTMWOlwSW5YvLapS0kXGKXOVeOi+pm2bdu6UWS66RewVwMUiYIUWbt2rQtKdLHyLliqh9AFTG306nYUAHiFr5H8+uuvMbXfG64c63E9aqPXTaFuH9Vg6NydeuqpifbTuQmlQFFt00XYe7/e8VT7JLrwqBZD8/107tz5hNdWUBONXkfHUbviSRdVjTp6++23rUWLFm6b/q0LqGqgUvIe00KBj75f3oU7KT2m4FfBvLr+FDSqPknfzZQWqCf9DKNJ+vmo5knfLwWJ6U2BvL43+n+ZlPd/Qe3wpk8QdQuH0h8fsmHDhpi+a8g5CIDg21FgobwaGRUuJ+Vt0z765a6i3ylTprjiadUQaC4U/bJVrYX3yzQc74KoYtVVq1YlekyFtrrYhl6wQrMrkcTa/pQeN9wosHCSHterr1Jgl/T9it6v6lSU3ZDQotWk7U/t+08Lfd6qJ1HhuLJdap9GDyow8rIksb7H1NK5UQFxnTp1kj0/nmeeecbVkqn2SgGv6nf0XVJxerdu3WJ+vZR+N8J9ZiqIDt3H+4yTBjBpkdLvvSQNCL39QtsLfyAAAv4vIyKag8X7S9+ji48uCsoaeRkQ/fXvZQA0akVFo+oG8YIEBQzJqVu3rvupbp4bb7wxU9qfkXSxUVGyugejBVDqbkkuU6L3FI0yIyoAVpdHJOE+l0iUVVL33Pjx4122UPPThHZ/peQ9poaK7dX1GsvUB8pCeV1eGpmlLKVGXnkBUGrefzT6zPT+PQoIFQjqtT36w0DBjjJZoQFIcp93StqoY6m7K7nviL73wggvhEMNEPB/2SHVimhkVWidhUbwKNNzww03uG4GjSZKOhKmcePG7sIYWheki6GyPElp5JfqZzRayasnCqXui/Rsf2bo27ev62IMnW7Ao1E63qgfjXbSKDJND+DROdTInlhouLZGJCUd0abPS3U44tVpJffZhKPuTXULqetLN2Vimjdvnqr3mFL6PmiEk7pE77vvvrD7KbjwukE9ZcqUcW0NrSlLzfuPRqPkQrMsChb1fjVa0uMFZQrmPArQklvOJaVt1OtomL1Gd3k0ceQLL7xgVatWdX+cAMkhAwT8X4pey1do3hld8DTviAIUDevVL9Cnn37anSf9Zashv7oA6uKiOg91hynDotmSQ4uPVUyqi7KGQusvVQ1HVxCiC7RqiDRcWFkkzW+iIb8qtFUhreaXSa/2Z4b777/fvT+1Te9bdRa6QCnQUXDwzTffuOyUho3rHGgYvM6dup80tFoBgOpbYnkdBaiaG0nz5yhoWb9+vRvWrQutMgWqiVFBseaCUXeRhvJ7c8xEooyPurK8YuvUvsdINCWC2qtgQhd/Ba86HzqeurO8+XzCda1qUkp9DxVg66eWlNB3LbQAPrXvPxJNe+DVwqmOTd2F+t6rRsyj4vaLL77YLW+hKQLU7aSARfslzWyltI0qNFdtneqA9L3Re9cEkpqWQufUm64BSIoACDmeukd08UraNZSU6jR0AdUvT6XmvV+kmiDQy57ol71S66r30U8FHo888oh16dIl0ZpNyrioS0cjUFTwHDp1v7oLVEekCd00j5AyFPXr13dZBAVVHk3ApzXLkqMLQtKailjaH+24ydFFKHTEWnL0/nSOFcAlpddW/YwudppbRsWmqnfSaC0Fk16Xh4JEXfAVrKgrSxkM1bMom6Bjh45uS25JBB1Hk+fpdZRJUgCoYFDZNi/4UHZBWRUFBhr5o2DDy9xFWmZBI4W8OWVCMxspfY/JUXbHC65UaO6tBaYRZZqfRoFbUknbqs9T50yTYCpw0ntXIKl5kELPW6T3r3OlYya39lhy58bbX/VQGgWpPwTUDanvtQKepJT9UTeiipLVpqlTp7qAT8cIfY8p/YxUzK/MmzJwqn3SvvpjRP9HQoPGcO9Pr63t0b7jyHlYDBUAAPgONUAAAMB3CIAAAIDvEAABAADfIQACAAC+QwAEAAB8hwAIAAD4DgFQGJrgTlPeh66sDAAAcgYCoAgrh2sSLf0EAAA5CwEQAADwHQIgAADgOwRAAADAdwiAAACA7xAAAQAA3yEAAgAAvkMABAAAfIcACAAA+A4BEAAA8B0CIAAA4DsEQAAAwHcIgAAASGdjxoyxt956K0ed58mTJ9uLL75o2RUBEAAA6WzWrFk2c+bMHHWeFyxYYF999ZVlVwRAAADAd/JmdgMAAMju/vrrLxs3bpxt2bLFatWqZTfeeKMVL178hP2+//57lzXZu3evXXnlldasWbPgY4cPH7YPP/zQfvrpJytdurRdfvnldvrppwcf37Fjh40dO9b+/PNPq1KlivXo0cOqVq2aqJstISHB6tWr57qn8uXLZ5deeqlr14gRIyx37v+f8/jyyy9dVmrIkCExHVt++eUXd6w8efJY69atLbsjAAKATDbs+ak283BC1P1aF8ht993fJUPahNgtXrzYLrjgAuvYsaM1adLEPv74YxdwLFmyxEqUKBHcb/r06bZo0SIXHG3bts1atmxpEydOtC5d/veZKuBZv3699ezZ0/bt22fXXXedvfLKK9a0aVNbvny5Czp002v8/vvvdvbZZ9t//vMfa9y4sXu+AhrdL1eunHuuApjq1avb66+/bldccYW1atUq2JYhQ4bYGWec4f4dy7EVuOlxBW1nnXWWPfjgg7Z161arX79+tv2q5AoEAoHMbkRWtGfPHhe9796924oVK5bZzQGQg3UaOMW2FipoZQ8eCruP9/i0R7tmaNsQnQKL8uXL2/jx4939I0eOuMyNggUvw9KrVy/3uDJFFStWdNseeeQRmzBhgq1YscJlf0466SRbuHChNWjQIHicXbt2WdmyZa1du3Z25pln2rBhw4KvO2DAAJs/f36wtkiv8emnn9rq1asTZZ86dOhglStXtjfeeMPdX7dunQuO5s6d6zJQsRxbwVqNGjVclkmUwapWrZo1bNgw29YBkQECgCwgWnCjIMmP+kxYYFv2hg8M00O5ogVt9JX/y3xEc+zYMRdIKJPjyZ8/vwt+lJEJdd555wWDH+nevbsNHjzYNmzY4AIUBRRDhw61hx9+2GVZdBwFPwcOHHCBiO7fdtttpryFbgp0li5desJrJO16u/baa61v3742cuRIK1CggH3wwQcuM6TgJ5ZjHz161ObNm+fa5SlatKgLrJTJyq4oggYAIJX++ecfFwSVKVMm0XbdVz1QKNX1JN1H1JXkdZEpeLnsssvs5JNPdgGJMkA7d+6048ePuy4rZVwaNWrkuqYUQIVmbSS0y83TtWtXF8R88cUX7v7777/vgiKJ5dh6j9onXPuzKzJAAIAsK9ZMTGZRoFKwYEH7+++/rXnz5sHtun/KKack2lddT6G0j6joWJQBGj16tPv3b7/95mqCnnzySXv22WddhkaP33LLLSluY+HChV2dkQKfmjVrumLmSZMmBYOYaMdWFkqZI7VfNUJJ259dkQECACCVcuXKZd26dXPFyqrjkY0bN7puJgUwSefNUV2NqJvppZdesvPPP98FUcqyKAPkqVu3rtWpU8eNzlLwoS615557LlFWSYXSGs0VCxVFT5s2zbVTGZ5atWq57bEcW+9RWalXX33VZZJEdUteRim7IgMEAEAaqG6nbdu2rpBYtTuq/VGmRF1YoRTUqGtJRc5r1651dTZe0KMh6zpO//793cgqjQZTFuibb75xj6t+RwGVjqERZ4cOHXKjtVRIHQsVOqt7TRkmBV6hYjm2apVUCK3RYdpPgZw3iiy7YhRYGIwCA5BRvALnWIqgGQWWNWnE1rfffmubN2+20047zQ1dDzV79mw3R4+Cmx9++MGNorrwwgtPqKtZtmyZ66IqWbKkCzg0MiyUhtYrMFLXlYIs7Zf0NUKHu4eaMWOGC7oU7JQMeV4sx/ayQhoar2BNtUIK0pS5at++vWVHBEBhEAAByCgEQIAPu8DU56i+Uv1UZKz0YN68kZulWTK//vprF41qxktVuCsi9WheBQ1LDFWhQgV76KGH0u19AACA7CNTi6BXrlzpgh5NoqQA5oknnrCLL77YDbcL5/777w8ODZTHH3/cpeKUsQlN8ykVeeqppwZvlSpVypD3BAAAsr5MzQApI6Mqd1Waa42S3r17uyF6yghdf/31yT5Hw/Sef/754H1N7qRpvz///HO7+uqrg9t1nHvvvTdD3gcAAMheMi0D5E3KdM011wQXaNOcCSre+uSTT8I+TwFTqIMHD7rhhEmLtbSgmyrYNX+CZrAEAADI9ABIQwA1Z4LWFgml++oai+SPP/5w2Z2bbrrJDe0bOHCg6zrzaM6CUqVKucmptO6KhicqUxSJ2qJutNAbAADImTKtC0zrj3jriYTSwqPeY+EUKlTI1fVs377dDflTwbMKoosUKeIe13olWujN06NHD2vTpo2bCVPDDpOjOQ6eeuqpOLwzAACQ1WVaBsgLfLxiZo/WJUkaFCWl4EYZoGeeecatnKvbiy++mOjxUK1bt3ZTjX/33Xdhj6mVb7Xyu3dLOmU5AADIOTItAFK9jzI26s4KpfuaZTJWmtlSdUFJj5PcJFXeFN7J0XTgyj6F3gAAQM6UaQGQCp81G+U777zjpt2Wn3/+2XVnXXXVVcH9PvroI1fI7AUxSbM4mtXyxx9/dFOLi1bl1TTkobx5hi666KIMeGcAAKQvzYWnGlePBhVl1uKkXydpS3aRqfMADRkyxI3i0jw+PXv2dHU6WrBNtToerYPy3nvvBYubVaejxeP69OnjAiWtvaIC6DvuuCO4z6BBg+y8885zw+o1RffNN99sTz/9tFvjBACA7O7uu+92gYdH18SZM2fG/Pxp06a5wUjp0ZbsIlPnASpfvrwtXbrUzQOkDI0+QAU3oVTA7K1roskSFRBpvRJlfQoXLmz//ve/Ew2Nz5Mnj1urxNunQ4cO9tZbb1nlypUz/P0BAJARLrnkEjc4KFa3/N+ceko6+FWmL4WhEV3dunUL+3hyo7bOPfdcd4skln0AAEirRYsWufnotDSTlmrSQqf6Yz50IVNvH9W4ajFU1aR6i4jq3wsWLHDP0+OqkU1q27Ztbh+tapDcKuzqOalWrVqibSoJ0etqsFHDhg3dIqeiRMLhw4fdACJNF6Plpy677LK4tSW7yPQACACA7GzkyJGuhlUjiBUwbNiwwU3nMn369GAPhbfP/v373ajk008/3QVA6gXRepYaeFOxYkUXHGklhJdeeilRd5V6QxSQqH5W07/otUKpB0Ujo3v16uXuqwdEdbbat3bt2rZ8+XJXT6v1NtXrcvDgQfv+++9dW5WIUAAUr7ZkGwEka/fu3QGdHv0EgPR0yTOT3S2t+yBz3HDDDe568c0337j7x44dC3Tq1CnQrl27RPvkyZMnsGjRouC2w4cPB6pWrRoYPnx4cNvGjRsDZcuWDUyaNMndP3DgQKBixYqBJ598MrjPM888417vtddeC26rVKlSYMyYMe7fBw8eDFSuXDnQq1cv1xbZv39/4MsvvwzuX65cucB7772XLm3JLsgAAQCyrN3dn7GETTsy9DVzVyhtxT96LEXP0Uhkr2RDtagPPPCAG3izY8cOK126tNveokWL4Ihlb+FuzTmnriQtAaUuMt20lqUeU3nI7NmzbfPmzW4hcE///v1d/Ws4yjytX7/ennvuOdcWUXdc6IoJSaVXW7IyAiAAANIoaf1N9erV3U+NtPICIHUrJZ3GJX/+/G66l1AVKlQIPl9D28uWLesG/XjUZaV9wtFztByUV/MTi9Xp1JasjAAIAJBlpTQTk1mSrmrg3Vcg4tE0LaG06oGKjjXViwKJ5Oj5ydXYJH29pBMEaz1LFUGrwDkWRdOpLVlZps4DBABATjBv3jy3PqVnypQprthZt3A0950KiTVVS6jjx4+7kVbSuHFjF5iEzrOj+X4iFR57U8eMHz8+0XZ1x3mKFCniRoKld1uyMjJAAACkkZZTUg2QJuVds2aNm2NnzJgxLqgIR/U2qtNRHY2Wc9KkwKrDmTBhglugu2PHjm5U2V133WXXXHONPfLII+54L7zwQqIh9klp3jtNCKy5fjTyTCPOtIqCusS8lRUaNGjggh0dxxsFlh5tycrIAAEAkEYKELS6wbJly2zTpk02depUu/baa4OPK6BQBiWpe+65xwUnCqA0P48yLKrD0fE8w4YNc8fW0HYd+7PPPnPBTY0aNcJOhKgibE0KrOH4Kl5WAbaO4Rk5cqS1a9fOHWvSpElxbUt2kUtDwTK7EVmR+k/Vj6rUHgujAkhPnQZOcT+nPdo1Tfsgc2juHdXbjBs3jo8gGyEDBAAAfIcaIAAA0kDdWyoWRvZCAAQAQBrceeednL9siC4wAADgOwRAAADAdwiAAACA7xAAAQAA3yEAAgAAvkMABABADnHkyBE7dOiQZXdHMuB9EAABABAnmhE6Fvv377eDBw8m+1hCQkKy27V4qZaniERrefXo0SOmNoQeN6sFTal5HylFAAQAyNKOb9yRobeU0npbumCXLVvWChYs6NbF+vjjj8Pu/8EHH7jV2Fu3bp1o+9y5c61atWpWtGhRt6hq6EpVe/futbp169ry5csjtkXreGlx05S455577LrrrjO/IQACACAN+vXrZ5MnT7YZM2a4rpuXX37ZbrjhBrcIaVKrVq2yBx980C666KITHrv55pvtiSeesPXr17vFSD/99NPgYw8//LBbXPWMM86I2Bat3K5V6EOzO7qFyy4dOXLEZa00k/W+ffvcLems1uGyTqHHVrCmjJb2VUCYlF5DWS/veOFeKyMRAAEAkAYTJkywu+66y+rXr2+5c+d2q6drdfYXXnjhhGBD3ToKUpTpCbVz505bsWKFC3JKlixpXbp0se+//949Nm/ePPv222/t0UcfTXHX0a233mpXXHGFW7C1TJkydtJJJ1m3bt3cgt/y7LPP2rvvvutWdS9fvry7zZkzxz32yiuvWNWqVV1GqUKFCvb4448n6uLTsS+//HLXZmW0mjZtajNnzrQSJUrYjh2JM2kK+tq0aeP+PXz48OBr6XlaSkTZr4xGAAQAQBrkyZPnhCyJ7itwCaUsTs2aNe36669P9hjKongZEQUa2qagqU+fPjZq1CjXvRWuPiiSadOmWePGjW3jxo0uyFq6dKk999xz7rF//etfLvN02WWXBbMyrVq1csHPsGHDbOLEie69KABToKeAKdTnn39uDRo0cAHPTz/9ZG3btrXSpUu7fT1q8/jx44PdbAqGvNfatWuXC9D0+vp3RiIAAgAgDa6++mp76aWXXBfY1q1bXUblq6++sm3btgUDli+++MImTZrkApnkFCtWzM466yzXfbZy5UrXpXbBBRfYoEGDrEWLFi5Tctppp7mfN954Y4q6jpo1a+ZqivLnz2+nnHKKde/ePZhdCkdZqscee8xltVQgXaVKFbv33ntt7NixifY788wzXdZJtU+ioE0ZqPfffz+4z3//+193LpIratb50Vpqel5GZ4EIgAAASIOhQ4da7969rW/fvq5QWcHLI4884i7quXLlsu3bt7ugRVmVvHnzusyHMjy6+OvfXpCkwEkZF3UVaf/KlSu7gGPIkCEuC6Ri5S1bttjixYvto48+irl9SbvbihcvHjHboiBuw4YNrltP3WYnn3yyK/C+77773HsJVadOnROery4xBTN///23u69gSDVPOo5XB3XppZe6Ym+1RV1hes1169ZZRmI1eAAA0kDZj6efftrdPMpqVK9e3QVAurCrADg0A6LiYQU+uvir4FlZGmWAfvjhB/e4usPOP/98V0ekGhwFPdOnT3dBg2p41L12zTXXxNQ+tSE1Jk2aZB06dIi4T758+U7Y1rBhQ5et0mg3r0A8NPOl83D66ae77jjVFokyU7FOIRAvZIAAAIgjjYbSRV8FwnLOOecEa168m+puVPyrfyv4SerVV191wZFqYxTA6OYFCF59ULzkz58/UZeasj0KSFTfk1rKAinzo+JqBXoq6ha9zo8//ujevxf8KEBUximjEQABAJAGypQ8//zzbvj6b7/95gIf1fQMGDAgVcfTcVRsPHLkyGCWRUGS7v/111+u+0v1QfFSrVo1++WXX2zt2rXBoekDBw60119/3RVLa7uyNa+99prLbMUaAP36669uWL8yVhp9Jgrc1E04evRoV5StwmllhFJT3J1WBEAAAKRBp06d3DB2dVmpy0jZEw0lVxAUqdvMCwqS0igpDTn3MiTy5ptvuq6yli1buoCia9euMU2EqNfxCpRDMz4nhby26o3U/aZh7N4weI3YUmCnDM65555rnTt3dsGdRo1FOrZH3X+a6FHBU9JRb+oa27Rpkyugvuqqq1w90Nlnn+3aFe59pIdcgdCpJhGkORJUnLV79+6IX2IASKtOA6e4n9Me7ZqmfQDEjgwQAADwHQIgAADgOwRAAADAdwiAAACA72T6RIh79+61Tz75xM1uqSm327dvH/U5ml1SU2truF69evVcpXnSiZ5Sc1wAAOAPmZoB0lwHCk609smff/7phuJdeeWVbgbMcDTTplbanT9/vv3xxx9ueJ0WbtPEU2k5LgAA8I9MzQA99NBDbm0QrRmiiZ600JoyOpp7QKvDJkfriWiBNi/jo/VWtEib5irQfAKpPS4AAPCPTMsAaaZJdVH17NkzuJZI7dq13aq3WgwunCZNmiTq7tLESbqvSZPSclwAAOAfmZYB0uyQBw4ccAumhdJ9bzG4cDQVuGaS1GSF//nPf1yGR7NUpuW4WphON4+ODQAAcqZMywCpgFk023KoEiVKBB8LR1meQ4cOuSDln3/+sR07dgSDl9Qed/Dgwe453k3dagAAZCfDhw93pSHZ3fAMeB+ZFgAVLlw42UyLlp7wHgunVq1a9swzz9ioUaNs8eLFrv7HWzQutcfVonXax7tpdVoAAKL58ssv3TpgSW+hg3NkzZo1dvvtt1vbtm3d4BwtBBpq+/bt7nH1aKiUI5QG8ai0QyupR6KBP1qzKyWef/75RGt8ZQWpeR/ZpgtMi8VpobNVq1a5wmbPypUrXc1OrMqWLetGfC1dujRNx1UNkVdHBADIOo7OT98LYVL5mtZN0f6abkWrpU+ePDnR9tBrihb/1GKjqkdV2ca0adOsefPmbkTzGWec4fbRAqSVKlVyK6n36tXLvv32W7dIqGj1dP1xfs4550RsS//+/ROVc8Ri1apVLvjym0wLgPLmzetWgH3vvfesT58+7r4ivtmzZ7ttns8//9w2bNjg9jl69Kj7kmlEl2fz5s0uivbm+Yn1uAAAxIsG5CjrE84LL7zgVmAfP3685cmTx12nli1bZv/+97/t448/dtki1bRqVXmVbMyYMcOmTp3qAqCNGzfawIED7fvvv4/aDj1n27ZtNmjQIHdfP1U2op4TPaZSEL1279693QCiN954w23X9dVr/4svvmgNGzZ07fOmk1FZyM0335zoPXrHPvXUU23ChAmul0XXXbX1iy++SLS6+5gxY2z69On2/vvv29ixY11AJyo50Wrz999//wmlKzl6GPyzzz7rImBNZNioUSM3Suviiy8ODmeXKVOmuAhZJ1UpwJtuuskqVKhgdevWtV27drmh7YqI+/btm6LjAgAQL7oederUyQUVClr69etnpUqVCj7+zTffuDnsFPx4unTp4upPRYFJ7ty5rVixYu6+nuuVctx1112uTEPZoWgUrGguPI+SBgqwdA3U9VPBlLJE+fLlc91wSh4o0aDskubZEwVL3333neuK0/u4/PLL7ffff3fvT4kEBVChx27Tpo3ddttt7tqsAUfz5s2zr776Kjg4SUaMGOHaILo2V69e3f1bdbyvvfaam89PJS06B74IgBQ1KsJU5KgUok6QTljoCdAJP+uss9y/FU0qGNIXSf2gFStWdKlCBTspPS4AAPGgoKdr167WrVs3O3LkiLvmvPnmm653QmUa3goG3bt3T/Q8BTQKnLRygeauK1OmjAs8lGWZNWuW+8Nf3Wpbt261q6++2u6++25bvXq1u+716NEj5vbpWqnroRd86fo5efJkFwCpbKR8+fKutyQ0u6NuugcffNAefvhhd18lJepaU5DkBUBeoKZjhWZ79LgyPV4ApFoelamMGzfO3Vc2KXSgkYIwnSf11CgQ8s1SGCVLlnQfcjiXXXbZCV80fRCh9T2pOS4AAPGgwOaGG24I3teFv06dOi67o9FMoi6mggULJnqe6lVFQZNXjKyskAIhBQTKGKluSN1Jup55WZtbb73VBS7nnXdeTO1r0KBBosyTnrto0aKw+yvzpGyM2qxMjnpfdNOIawVyodTLEhr8eLVMOic6jjJaCnyUFfPKVzSK+6233nI1Tgru1I2m11L2ylcBEAAA2VlygY2yKd7gHC9TogAilO4rMPG6vbS0k7qJFBQogLrjjjtcwKMBPKrT0Rx4yhpplQOVh8QaAHmTAocmEgIRlobypoxRzY/qc5I+N1Ryo6v1HlTvpMyQAsMPP/zQdeN5dNyff/7ZdcUpGFOxuAKmpKPm0hsBEAAAcabyCwUBHtWqKqsSasGCBW4EWGiAoq4w3dQdpFoaPUcZE2WJVBzt9XDEc6qW3EnKQ5R9KlKkiOvyilTYHY7ej2pu1Q1Ws2ZNN0GxuvAkISHBBUaffvqptWvXzm1T4JM0OMwIFMUAAJAGmodOo7c8Kgz+73//m6jmRwXIGtmloEaWL1/u6nK0PSkFHioq1kgpdS8pkKpWrZqrC1LmRl1H3tD5eChbtqwbbe1RPZC62YYOHZooi6UMlEaNxUJ1SjoH6tZTkbTqkLxgS8HVL7/8EgyItH5nSofuxwMZIAAA0qBo0aJusI66srQUk7IZuvBr4kKPantUUKxgoEaNGq7eRbUyd9555wnHU6GxJkts1qxZcNtzzz3nggovO6OusXjp0aOHvfLKK250tbrqNAxeQ9yVmVEbFHwpC6XMjh6LhQYnVa1a1XXdvfPOO4keU12Uapq0XedKo85iGeEWb7kCkToCfUzFW5qTQEMDvf5ZAEgPnQZOcT+nPdo1Tfsg86iQV8PCVdOjgCFp3Y1HF3xvXh0NG0+OusZUMJy0vkbP1YSKqg9SliY5ytIom3L66acHJwHWZT50fUwNk9+xY0dwhLUo2NG+uvZpcmFvTh7d1/tSYbYCmtAaoOSOHUpZLs1JpCJsr+Dbo9FvmoCxdOnS7nwtWbLEnQ/vnCR9H+mBACgMAiAAGYUACMh41AABAADfoQYIALKJrYUKWpe3Z0fcp22tcnZ3i9jXUwT8igwQAGQDLddvtrIHD0XcZ/PeQzZj5ZYMaxOQnZEBAoBsoPevK9yt5Oz/zSycnGjZIQD/HxkgAADgOwRAAADAdwiAAACA7xAAAQAA3yEAAgAAvkMABAAAfIcACAAA+A4BEAAA8B0CIAAA4DsEQAAAwHcIgAAAgO8QAAEAAN8hAAIAAL5DAAQAAHyHAAgAAPgOARAAAPAdAiAAAOA7BEAAAMB3CIAAAIDvEAABAADfIQACAAC+QwAEAAB8hwAIAAD4DgEQAADwnbyZ3YAjR47YjBkzbMuWLVa/fn1r0KBB1OesXbvWfvjhB8ubN681btzYKlWqlOjxOXPm2K+//ppoW6lSpeyqq66Ke/sBAED2k6kB0LZt26xNmzYuCDrzzDOtX79+1r17dxs1alSy+wcCAbvmmmtc8NOwYUM7cOCAXXvttTZo0CC79957g/u9//779s0331i7du2C25IGSQCQEfYPGW9Hvvgh8k5tm5vlzcMHAvglABowYID7+eOPP9pJJ51kixcvtkaNGtmll15ql1xySbIBUJcuXVyAkzv3/3rv3nvvPevVq5d7To0aNYL7nnPOOWEDKQDIKK/8vdNmRwlwthYqaOUswIcC+KEGKCEhwT7++GO78cYbXfAj6v4677zzbPz48ck+R0FPjx49gsGPtG/f3h1r+fLlifbdunWrjR071qZOnWobN25M53cDAMmbXbm8bS1a2PJUKB32VqFEYWvX6P//AQcgB2eA1q1bZ3v37rW6desm2q77ixYtivk406ZNszx58rgutFB///23ffXVV7ZhwwZbsGCBPfvss3b33XeHPc7hw4fdzbNnz54UvR8ACKfswUM2tX8HThCQhWRaBsgLMEqUKJFoe8mSJWMOPn7//Xfr37+/u1WuXDm4XVmlP//80z744AObNWuWvfrqq66+SF1s4QwePNiKFy8evFWpUiXV7w0AAGRtmRYAFSpUyP1UFiiUgh+vSyySv/76yy666CLr2LGjDRkyJNFjTZo0cSPEQgOiMmXK2Ndffx2xHmn37t3BmzJUAAAgZ8q0LrCqVata/vz5bfXq1Ym2637NmjUjPlf7tGrVytULqQg6tCYonIIFC9rOnTvDPl6gQAF3AwAAOV+mZYDy5ctnF198sSt41uguUbHyzJkzrXPnzsH91IWlYmnPmjVrXPDTrFkzNxpM9T+hjh8/7vYJ9d1337m5gxQwAQAAZOoweBUmKyjR0PamTZu6UVvqvrruuuuC+yjDM3/+fDeJ4aFDh6x169buZ8uWLe3NN98M7qegqE6dOi6Y0pB4DaevV6+eC3zeeustN1+QXgcAACBTAyAFLD///LMLfDQT9P333289e/ZMVL+jwEbdZV52R8Pe5Zdffkl0LM0iLXrukiVLbMKECW5+oXLlytmXX35pLVq0yND3BgAAsq5cAa//CScUY2s0mAqiixUrxtkBkCqdBk5xP6c92jVNZ3Bny37uZ8nZw8Pu0+Xt2e7n1Jtapum1AD/I9LXAAACxSdi4IxgIJed4+5aWqxCDOYBYsBo8AGQD+Ts2sdwVS0fe6dhxCxz8/xO6AgiPDBAAZAOFH+7hbhH9X3cbgOjIAAEAAN8hAAIAAL5DAAQAAHyHAAgAAPgOARAAAPAdAiAAAOA7BEAAAMB3CIAAAIDvEAABAADfIQACAAC+QwAEAAB8hwAIAAD4DgEQAADwHQIgAADgOwRAAADAdwiAAACA7xAAAQAA3yEAAgAAvkMABAAAfIcACAAA+A4BEAAA8B0CIAAA4DsEQAAAwHcIgAAAgO8QAAEAAN8hAAIAAL6TqgBo0KBBtmHDhvi3BgAAIKsGQG+//bZVrVrVOnbsaBMnTrQjR47Ev2UAAADpJG9qnrRy5UqbPXu2C4RuuOEGO+mkk+y6666zm266yerXrx//VgIAYnPsuO1s2S/iLvk7NrHCD/fgjMLXUpUBypUrl11wwQX27rvv2ubNm12X2Pz58+3MM8+0Ro0a2ahRo2z//v3xby0AIPzv5kIFzPLmiXiGEjbusCNf/MBZhO+lKgMUqkiRIlatWjV3W7x4se3bt8+efPJJe+yxx1yAdMkll/j+JANARshdoohZiSJWcvbwsPtEyw4BfpHqUWBr1qxxgY4Cny5dulj+/Pntv//9r/3++++2bt06GzBggN16661RjxMIBGzhwoU2bdo0d8xY7Nq1y2bOnGlz5syx3bt3x+24AADAH1IVAF144YVWvXp1++yzz+yhhx6yTZs22TvvvGPnn3++ezxfvnzWr1+/qCPF9uzZYy1atLDOnTvbsGHDrF69ei5zFI6Cmr59+1rdunXt6aeftgceeMBOOeUUl2lKy3EBAIC/pKoL7LTTTrPnnnvOzjnnnLD75M6d2wVGkSgo2bJli8salShRwmbNmmWtWrWyNm3auFtyAVDt2rVt9erVVqBAAbdt5MiR1rt3b7d/lSpVUnVcAADgL6nKAI0dOzZs8KOaIE/58uXDHkPBzLhx4+zmm292QYqosFpF1NqebGNz57a77rorGPzIVVddZUePHrVffvkl1ccFAAD+kqoMULgRXgpEdIvF+vXrbefOnW7kWCjdX7p0acxtUd2RRqWdfvrpaTru4cOH3S20Gw0AAORMKQqAQjMoSbMpCQkJ9sMPP1itWrViOpZXvFyyZMlE20uXLu2KnGPx999/2z333GN9+vRxxdhpOe7gwYPtqaeeiul1AQCAjwKg+++/P9l/e4XPp556qr366qsxHcvrxjpw4ECi7RpGX7BgwajPV33RRRddZA0aNLARI0ak+bgatda/f/9EGSCvpggAAPg4ANKkh3LGGWfYsmXL0vTCCi7y5s1ra9euTbRd9zXCLFo7VMysrM/kyZPdEPy0HleBU2htEQAAyLlSVQOU1uBHlI1p3bq1W0vsxhtvdNt27NhhM2bMcEPXPZpc8Z9//rF27dq5+xrdpeBHa5F98sknJ2R1Yj0uAOREm/cesi5vzw77+PH2La3l+s32YIa2Csh6cgU0bCoGmudHevXqFfx3ONonFkuWLHHz9WgkV7NmzeyNN95wRdSqJfKyMbfccotbZkNBl4qU1eW1bds2F8yEBj8NGzZ0XXCxHjcadYEVL17c1RQVK1YspucAQFKdBk5xP6c92jXdT86IOcttxsotEffZtGu/lT14KEPaA+SIAKhy5crBUVbev8PRPrH6448/bPTo0S6zo4VU77zzTitatGjwca0rtmLFCnvhhRds7969waxOUrfffru1bds25uNGQwAEILsFQNmxPUCWD4D8hgAIQE4MOLJae4BstxZY0mBh/PjxtmDBgngcDgAAIOsFQCow7tGjR3D+HxUdq1bnvPPOY7ZlAACQMwOgZ555xv71r3+5f3/33XeuKFm1NlOmTLGhQ4fGu40AAACZHwCpKLlGjRrBpSi6du1qhQsXdkPV//zzz/i2EAAAICsEQBUrVrS5c+fasWPHXHeYN/pq3bp1VqlSpXi3EQAAIPMnQtT6W5dccomVKlXKzZHTvn17t12F0FdffXV8WwgAAJAVAqC+ffta48aN3WKkWo/Lm1xQ2Z8rr7wy3m0EAADI/ABImjRp4m6hbrrppni0CQAAIGsGQEuXLrV58+a5dbqSeuyxx9LaLgAAgKwVAI0YMcLuvfdeq1WrlpUsWfKExwmAAABAjguAnnvuOVfwrMVGAQAAfDEMXouSdurUKf6tAQAAyKoBkEaAzZ8/P/6tAQAAyKpdYC1atHBrgfXr189q1qxpuXLlSvT4FVdcEa/2AQAAZI0AaNiwYe7ns88+m+zjBEAAACDHBUC7du2Kf0sAAACycg0QAACALwOg3377zQYMGGDdu3cPbvv444/t4MGD8WobAABA1gmAZsyYYQ0bNrRff/3VBT2en3/+2UaOHBnP9gEAAGSNAOiRRx6x0aNH26effppo+7XXXmuvv/56vNoGAACQdQKgZcuWWbdu3dy/Q4fAn3LKKbZ27dr4tQ4AACCrBEBFixa1TZs2nRAALViwwCpWrBi/1gEAAGSVAEhrgGkxVG8l+ISEBJs1a5b17t3bTZAIAACQ4wKgwYMHu6Dn5JNPdj+VEWrVqpXVqVPHnnzyyfi3EgAAILMnQixcuLB9/vnn9uOPP9qiRYtcEHTuuedao0aN4tk2AACArBMAec455xx3AwAAyPEBkOYBmjx5sq1evdoVQVerVs0uv/xya926dfxbCAAAkNk1QH369LELL7zQPvvsMzty5IgdPnzYpk6dam3atLE77rgj3u0DAADI3AzQpEmTbPz48W4CxEsvvTTRYwqCNBFi+/btrUuXLvFuJwAAQOZkgN555x3797//fULwIwp6nnrqKXv77bfj1zoAAIDMDoCWLFliXbt2Dfu4ZofWPgAAADkmANq+fbtVqlQp7OOVK1e2bdu2xaNdAAAAWSMAUtFz3rzhy4by5cvniqIBAABy1DD4+++/P31aAgAAkBUDoHr16tlXX30VdZ+U0nxCW7Zssdq1a1vJkiVjes769ett1apVdvbZZ1uJEiUSPbZixQrbuHHjCbNXM1M1AABIcQC0bNmyuJ61gwcPusVTNbGiJlNUQDNo0CDr169f2OdoxXntM3/+fBc0zZw5061DFuqFF15wQ/ZDgzEdf8yYMXFtPwAA8OFSGGmlhVO1ntiff/5p5cqVc5Mrdu7c2Zo1a2ZNmzZN9jnLly+3G264wV566SU79dRTwx77ggsusIkTJ6Zj6wEAgK9Wg48XzSt0yy23uOBHNL/QmWeeGTFTc/3117uh+Hny5ImaXVK2SAHTsWPH4t52AACQfWVaBmjDhg22detWa9CgQaLtuq+sUFpNnz7dNm3aZJs3b3brlY0aNSrZCRw9Gr0WOoJtz549aW4DAADImjItA7Rz5073s1SpUom2lylTJvhYanXs2NEVQWtSRhVLq8use/futnLlyrDPGTx4sBUvXjx4q1KlSpraAAAAsq5MC4Dy588f7KoKdeDAgeBjqaU6otKlS7t/586d255++mkrUKCAqzEKZ8CAAbZ79+7gbd26dWlqAwAAyLoyrQtMs0YrOFFXWCjdr1q1alxfS/VCyjQlHRofSgGSbgAAIOfLtAzQSSedZM2bN3cry3v27dvnanfatWsX3KYi5oULF8Z83EAgYPv370+0TcdYs2aN1a9fP06tBwAA2VmmDoMfOHCgtW3b1h544AE39H3kyJFWvnx56927d3Cf5557zs35481BpKLmP/74I7jm2NKlS91PDYnX7ejRo27Cw169erl5gNauXWtDhgyxxo0b29VXX51J7xQAAGQlmToMvkWLFjZ79mw3oeHrr79u5557rs2dO9eKFCkS3EezQ4fO4KyAR/MHvfLKK26un08++cTd//bbb93jqh/S5IjKJmnk1/fff2+PPfaYzZkzJ821RQAAIGfIFVCfEU6gYfAaDaaC6GLFinGGAKRKp4FT3M9pj3bNOu05dtzemzE34n75Ozaxwg/3yLB2Ab7KAAEAMlauQgXM8kaeSDZh4w478sUPGdYmwHc1QACAjJW7RBGzEkWs5OzhYffZ2TL8eoxATkEGCAAA+A4BEAAA8B0CIAAA4DsEQAAAwHcIgAAAgO8QAAEAAN8hAAIAAL5DAAQAAHyHAAgAAPgOARAAAPAdAiAAAOA7rAUGAD6zee8h6/L27LCPH2/f0lqu32wPZmirgIxFBggAfKRtrXJWvmjBiPtsLVTQZlcun2FtAjIDGSAA8JG7W9R2t0g6DZySYe0BMgsZIAAA4DsEQAAAwHcIgAAAgO8QAAEAAN8hAAIAAL5DAAQAAHyHAAgAAPgOARAAAPAdAiAAAOA7BEAAAMB3CIAAAIDvEAABAADfIQACAAC+QwAEAAB8hwAIAAD4Tt7MbgAAIOvZWqigdXl7dsR92tYqZ3e3qJ1hbQLiiQwQACCRlus3W9mDhyKelc17D9mMlVs4c8i2yAAByDb2DxlvR774Iep++Ts2scIP98iQNuVEvX9d4W4lZw8Pu0+07BCQ1WWJDNCuXbts1apVduTIkZifs3v3blu2bJnt378/rscFkHUp+EnYuCPiPno8liAJgL9lagB07Ngxu/nmm61cuXLWsmVLK1u2rI0dOzbic3777Te75ZZbrEaNGla/fn1buHBhXI4LIHvIXbG0y0yEu+lxAMjSAdCgQYNs2rRpLqjZuHGjjRgxwm688UZbunRp2OfMmjXLmjRpYnPnzo3rcQEAgH9kag3Q6NGjrXfv3i6bIz179rTBgwfbm2++aSNHjkz2Obfffrv7uX79+rgeFwCQuCtxZ8t+YU/J8fYtLVehApwyZFuZlgHavHmzbdiwwWVzQjVr1swWL16c5Y4LAH6hIvKoXYnHjlvg4OGMahKQczJAO3b8r5CxdOnE/8l0f/v27Rl+3MOHD7ubZ8+ePaluAwBkZxpBF3UU3cApGdUcIGdlgPLm/V/slXSEloKQfPnyZfhx1UVWvHjx4K1KlSqpbgMAAMjaMi0AqlSpkuXKlcs2bdqUaLvupyX4SO1xBwwY4IbWe7d169alug0AACBry7QAqEiRItaoUSP74osvEmVppk+fbq1btw5u0yguzeUT7+MmVaBAAStWrFiiGwAAyJkydRTYU089ZZ06dbIzzjjDFSkPHz7cChcubLfddltwn8cff9zmz5/vJj30JjfUCLAtW/43Bfvq1autTJkybq4f3WI9LgAA8K9MnQfo4osvdvP1aG6ffv36udqb7777zkqUKJGoS6tWrVrB+9q3R48eds8991i9evVs2LBh7v7EiRNTdFwAAOBfmb4WmIIV3cJRNidUly5d3C2txwUAAP6VJdYCAwAAyEgEQAAAwHcIgAAAgO8QAAEAAN8hAAIAAL5DAAQAAHyHAAgAAPgOARAAAPAdAiAAAOA7BEAAAMB3CIAAAIDvZPpaYAAgI+Ystxkrt0Q8Gcfbt7SW6zfbgxH2eaPeaTa7cnnL8/bsiMdqW6uc3d2iNicf8CkyQACyBAU/m/ceirjP1kIFXXATiR7XfpHodaIFWwByNjJAALKM8kUL2tSbWoZ9vNPAKTEdp+zBQza1f4ewj3eJkh0CkPORAQIAAL5DBghAtqLurUgZHD2uDBAAREIGCEC2oQLoaMGNHtd+ABAJGSAAWULCrn0WOHjYdrbsF3afmzfusN4VS1vJ2cPD7hPp+QDgIQMEIEtQ8GPHjkfcJ3fF0pa/Y5MMaxOAnIsMEICsI2+eiNkdAIgXMkAAAMB3CIAAAIDv0AUGAMjSS6AIS5cg3sgAAQCy9BIoLF2C9EAGCACQpZdAufSFL+z4ph1RpzjQCMHCD/dIhxYiJyIDBADI9lMkJGzcYUe++CHD2oTsjwwQACDbT5HABJhIKTJAAADAdwiAAACA7xAAAQAA3yEAAgAAvkMABAAAfIdRYAByHA2JjjQq6Hj7lparUIEMbROArIUMEIAcRZPh5a5YOvJOx47/b24ZAL6VJTJAR44csT179ljp0qUtV65caX7Orl27bN++fYm25cuXz8qVKxfXdgPIejQTcNTZgAdOyajmAMiiMjUDlJCQYPfdd5+VKFHCTj31VKtcubJNmTIlzc95+OGH7bTTTrOmTZsGb9dee206vxsAAJBdZGoANGzYMBszZozNmzfPZXMeeugh6969u/3+++9pfk7Hjh1t/fr1wdv06dMz4B0BAIDsIFMDoFdeecVuueUWO/vssy137tx299132ymnnGKjR4+Oy3P++ecf11UGAACQJQKgrVu32t9//23NmzdPtP3888+3BQsWpPk56harVq2aFSlSxM477zxbsmRJOrwLAACQHWVaEfS2bdvczzJlyiTarvvq3krLcxo2bGi33nqrnXPOObZ792678847rV27dvbrr79a+fLlkz324cOH3c2j7jUAQATHjkddhFSj8qIWpQN+CoDUfSXHjh1LtP3o0aOWJ0+eND1HXWSe4sWL2xtvvGEnn3yyTZw40e66665kjz148GB76qmn0vCOAIQzYs5ym7FyS8QTtLVQQSt78BAnMZvQPEpbLJdd375l+J2OHbeWf2+zBzOyYUBW7wKrVKmS+7l58+ZE27ds2RJ8LB7PkUKFClmFChVszZo1YfcZMGCAyxZ5t3Xr1qXo/QAIT8HP5r2RgxsFPy3XJ/6/jayrXaMaVqFEYctToXTY29aihW125eSz7oBvM0DFihVzhczffPONG8XlZXZmzJhhffv2De63c+dOV8isOXxifU4gEEg0N9CmTZtc7VD16tXDtqdAgQLuBiB9lC9a0KbeFD5bEK0rBVnL3S1qu1sknZhvCVlYpk6E+K9//csFMo0aNbJmzZrZ888/77bffvvtwX0eeOABmz9/vi1btiym56iOp1WrVu559erVs7Vr17rsjuYLuu666zLlfQJ+l7Brn5t5OVKQo+Uros7gjBxXJ8SyJPDlMPhu3brZuHHj7N1337WuXbu6wuNZs2a5eh1PqVKlEhUuR3uOsjivvfaaq/e57LLL7IknnrC2bdva4sWLXQYJQMZzy04cOx5xHwU/KphFzuHWW8ubfE1nEMuSwK9LYSib43VnJWfo0KEpfo66yT744IO4tRFAHOTNYyVnD+dU+kjuEkXMShSJ/LnTTYZMwmKoAADAdwiAAACA72R6FxiA7I05fgBkR2SAAKQJc/wAyI7IAAFIM+b4AZDdEAABSBPm+AGQHdEFBiBNmOMHQHZEBghAHH6TMMcPgOyFDBAAAPAdMkAAgHSzee8h6/L27LCPby1U0MoePBT1ONov0nG0pljL9ZvtwVS3FH5DBggAkC7a1irnRghGouBHgUskejxakKQAaXbl/79uJBANGSAAQLq4u0Vtd4sk0krxnt6/rnC3SGuKdWJNMaQQARAAAP9n/5DxduSLH6Kej/wdm1jhh3tw3rIxusAAAPg/Cn4SNu6IeD70eCxBErI2MkAAAITIXbF0xO62WLrtkPWRAQIAAL5DAAQAAHyHLjAAYY2Ys9yt9h5JrPO4AEBWQgYIQFgKfjSRXVrncQGArIYMEICINJHd1Jtahn2cglAA2REBEJAD5yiJRSzzmCTs2udWe48U5GhIsEbNADmhSzfWJTX0vY8W/DNXUNZGFxiQw+YoiUWs85go+LFjxyPuo+BHv+iBnNClG8uSGvq+Rwv6mSso6yMDBOSwOUpikaJuq7x50vx6QHbp0o1lSQ1lTqNlT+kazvrIAAEAAN8hAwT4tM5B8rw9O+J+DHFHRohWT5PRdWb63neJ8n9DK91HW+gVWRsZIMCHdQ6xYog70lss9TQZWWemAuho81rp/1e0PzKQ9ZEBAnw8dD1abQ91DEhvsdTTZKTev65wt0j/N6Jlh5A9EAABAHKGY8cjBu1e129W6m5D5iEAAnLgHCWA3+QqVOB/0zakEdM6+AcBEJDFanfUfRVtjhICICCx3CWKmJUoErHryiv6Z1oHCAEQkMPmKAEAREcABGRA11Usw2ZjWXbC2ja3rUULRyzCjJZFSslU/jm5HoKhzjmPvvvx+L8RD2/UO81mlz/ZLMofLa0L5Lb77u9iOel3WXbBMHggA4adxzJsNpZlJ1pu3mblLBBxH/2C1y+otA49zsn1EAx1znn0nY8W3MTyfyNe5tSq6v5YiRaEzzycYDntd1l2QQYIyICuq5iHzUZZdkK1Pw/mwKHHGY2hzjmPMhJZKSuhmqQKZja1f4ds1aVdPl6/y7KBLJMBSkhISJfnpOa4AAAgZ8v0DNDgwYPtxRdftO3bt1vdunXtpZdesjZt2qT5Oak5bkbY3f0ZS9gUfTVvdTtklb/Q9w8ZH3XlcNffHWUFZa/rQX99RzqOUsduREc26YOOpXZHw9c1TBfZR0w1WVns/yqyYb1RlLmL4vUdi+X3+PE4/Z6K5bUkd4XSVvyjx8yXGaBRo0bZoEGD7P3337fdu3dbt27drFOnTrZ69eo0PSc1x81KVHgay5cno6gtalMkCn7Unx2JN4Q74nHKn2xbLFe26oOOpXZHj8djjhJkrc81q/1fRfaqN3LBRt48GfIdi+X3uMXp91RMr+X3DNALL7xgN998s1144YXu/pNPPmljxoxxAcyzzz6b6uek5rgZJZZoNysuP6Bi2Gjza0Tr73Z/LZUoHHkOjoFT3Do8UY+T1USp3Yk2EgRZVJTPNSv+X0X2qTeKZe6ieH7Hov0etzj+nor6WllApmWAduzYYStXrrQLLrgguC1Xrlzu/vfff5/q56TmuAAAwF8yLQO0Zcv/ujBOPvnkRNvLli1rCxYsSPVzUnNcOXz4sLt59uzZY+mhz4QFtiXKMEO3Xo1S71kla9C2uftL2JtFNS393dH6zdVNVnbv/qj1NNovq4ygcG2Osnq0t1+kNsd6HMRH1HmQYphzKZb/q3yuSMvvxLhdD2L4Pb41Xr9/Y3gtKVe0oI2+srH5dhRY0lFauq+MTVqfk9Ljqmi6ePHiwVuVKlUss8TSL5yh8uaJWhgXS393LP3mmuNGc92kdQ6XjKS2aDKzSPR4tDbHchzERyzzIMUy51Is/1f5XJGW34lxux7E8Hu8XLx+/8bwWllBrkAgEPl/eDrZuXOnlSpVyiZMmGBXXHFFcPu1115rGzZssG+//TZVz0nNccNlgBQEqYi6WLFicXznAAAgs2Xan5slS5Z0w9NnzpyZKEuj+82bNw9uO3bsmB05ciTm58R63KQKFCjgAp3QGwAAyJkyNd/+0EMP2dtvv22TJk2yjRs3Wv/+/W3fvn12++23B/e57bbb7Nxzz03Rc2LZBwAA+FemDoPv2bOnC0wGDBjgipfr169v33zzjVWuXDm4T758+Vx2JiXPiWUfAADgX5lWA5TVqQZIxdDUAAEAkPMw5AQAAPgOARAAAPAdAiAAAOA7BEAAAMB3CIAAAIDvEAABAADfIQACAAC+QwAEAAB8hwAIAAD4TqYuhZGVeRNka0ZoAACQvRQtWtRy5coV9nECoDD27t3rflapUiV9PhkAAJBuoi1lxVpgYSQkJLiV5KNFkIhMGTQFkevWrYv4RUTacJ4zDuea85yT7MnBv6PJAKVS7ty5WT0+jvQfK6f958qKOM+c65yG7zTnOb1QBA0AAHyHAAgAAPgOARDSVYECBeyJJ55wP8F5zgn4TnOec5ICPv4dTRE0AADwHTJAAADAdwiAAACA7xAAAQAA32EmaKS79evX25o1a+zss8+2IkWKJLuPHt++fbvVqVMn7D6I3U8//RSczdxTrlw5q1WrFqcxDbZs2WJr1661atWqWZkyZTiXcXb48GFbuHDhCdvr1atnJUuW5Hyn0aFDh2zJkiVWsWJFO/XUU8POnrxy5UorX758zp8LLwCkk3nz5gU6d+4cKFOmjBZWCyxcuPCEffbt2xfo0KFDoHDhwoHatWsHTjrppMCbb77JZ5JGDRo0CFStWjXQvHnz4G3gwIGc11RKSEgI3HHHHYECBQoE6tat634+9NBDnM84W716tftdoe9v6Hd37ty5nOs02LZtW6B///6BChUqBAoVKhS47777kt1v+PDhgYIFCwZOP/10t98VV1wROHToUI4992SAkG5++eUXu/HGG23QoEF2xhlnJLvPww8/bMuXL3cZIP1FPW7cOLvhhhusSZMmYZ+D2Nx2223u/CLtRo8ebe+9954tXrzYZSMWLFhgLVq0sHPPPdeuuuoqTnGcjR8/3mrWrMl5jZO///7bZX30O7lt27bJ7jN37lzr37+/TZs2zTp27Ogy940aNbKBAwfav//97xz5WVADhHTTp08fu+yyyyxPnjzJPn7s2DEbO3as3XHHHcHuhOuuu86qVq1q77zzDp9MGimVre4ErWmHtHn77bft8ssvd8GPNG7c2C666CK3HfGnP4h+/PHHE7pxkToNGjSw++67z0qXLh12n7ffftvtp+BH1P3Vq1evHP0dJwBCpvnrr7/cQnz6TxdKf3Xolx/SZsSIEda7d29XV6VzvGzZMk5pKgQCAVdTlfR7qiCI72n66Nmzp/tjSBfsW265xQ4cOJBOrwSPvsvJfcc3bNhg27Zts5yILjDE7I8//nCFypGo6ypfvnwxHe+ff/5xP5P+VaL7K1as4JNJ8stp//79Yc9J3rx5rWnTpsH7d999t1155ZVWqFAh91f01Vdf7bJxSoFrG2Kn867i3OS+p953GPGh7+aUKVPcd1V+/fVXa926tRUuXNheeuklTnM6+ueff5L9jnuPnXzyyTnu/BMAIWbvv/++zZw5M+I+n376qZUqVSqm43mBkkYmhDp48KDlz5+fTybE8OHDXcYsnKJFi9qXX36Z6C/o0Meef/55O/30023RokWudgWx43uacTRS0Qt+RF2OCuZffvllAqAM+J4fSuZ3seTU38cEQIjZ008/HdezpVofUYpV3V4e3T/llFP4ZEKoViqtFxbv3CJltEaSzl/Sc8f3NGPo3G/dutXVDCrTifRRtWrVZL/jquFUAXVORA0QMo0KnzU3kLJGnp07d9rs2bOtXbt2fDKppL/aVLcS6uuvv3Y/GVmXOvo+fvbZZ8H7x48fd6Nl+J7GV3LdvPru1q5dm+AnnbVr186++eabRFmgqVOnWsuWLXPsQqmE00g3mzZtsj///NNNHCcqJNV/Lk0iV6lSJbdNQ+QvvfRSl/E555xzbNiwYVajRo1EXThIea2WRtZpCgJNdqb6IZ1n3ScASp3HHnvMZSlVkNulSxf74IMPXF3EAw88wNczjoYOHep+X1x88cWu7mfSpEnuIjxhwgTOcxooezZ//vxgkKnMznfffee6x88666zgtBmjRo2yrl272u23325z5sxx3erRyh6yM1aDR7qZPHmyvfDCCydsv+uuu6xHjx7B+/oPpv94O3bscEHQQw89xCy7aaRi59dff93N6Kr0tS7aobUVSDkV5CpA15wqCtIffPBB5qqJM2UuP/74Yxf07Nq1y0477TR3MVYGCGmbEuOSSy45YXudOnXszTffTPRH65AhQ9yIUc0Erd/VzZo1y7GnngAIAAD4DjVAAADAdwiAAACA7xAAAQAA3yEAAgAAvkMABAAAfIcACAAA+A4BEAAA8B0CIABZgmYJHz9+vJu0LbNozakvvvgiW7Q1Hvbt2+feR3JLUKREepwPTYT4ySefxO14QFJMhAjEiWawTUhIcP/WFPO1atVyM9n6nc6Jzk0kmlm5SpUqVqFCBTeLdWYt2XH11Ve7Fci19EUkmzdvzvS2xsOqVavc93T16tVu2ZTUSo/zoe/NmWee6ZZx6dy5c1yOCYRiLTAgTq655ho799xzrXr16rZnzx6bNWuWtW/f3q1jpBWV/UoLh4b+Ja/14RYvXmxXXXVVcFubNm1csNi9e3crUaJEprRz6dKlboHT0aNHR923UKFCmdrWrCY9zkfu3LndciMDBgwgAEK6IAAC4qhPnz5uwUz5+eef3dpmY8eOdQuRhq53tGTJElu3bp1bGNZbjDCU1vDSoqZax+vss88OBlAKqooXLx5c5PTo0aNuteaCBQsmer62z5s3zy3Yqb/I9Vd+KO84CtZ0HC2W2KRJEytSpEhM7Yj1fUi+fPlc94hH677pNUO3ed0oWq9M2TOve0YBiRbL3bhxo2tH1apVXVbAW1x3zZo17v0pg5SUumO8BSDV9nLlylkkI0eOtG7dugVf3zuPOoYCWh3DW8RXq2OHttWjNq5YscK1R+tXTZw40QXBJUuWTPR+lDH5/fff3SLA3vtJT3v37g2ei6ZNm57Qbi8wjdQmrRO1YMECF+ycd955ib4r4c5HtNeN9hldfvnlbi2wb7/91lq1apWmcwCcIAAgLvLkyRN44403Em2rUqVK4J577gne37JlS6BZs2aB6tWrBzp37hw49dRTA61btw7s2bMnuE/fvn0DJUqUCFx66aWBJk2aBJo2bRrYvHmze6x9+/buvo6rf1erVi1Qs2bNwLp164LPX7VqldtWo0YNt0/hwoUDt99+e6J2abvaof06duwYqF27dqBq1aqBjRs3xtSOWN5HOK+99po7V0lt2rQpoF9Jv/zyi7u/cuVKd/+CCy4InHnmmYEOHToE8uXLF3jggQcCV155ZeDss89276NAgQKBcePGJTrWRx995Nqu52qf4sWLu9eNpHz58oF33nkneP/vv/9276t+/fruHOi9PvPMM8m2VR5//HHXlgsvvDBQr149117ts3DhwkTvR8fS45dccon7bPr16xexXZMnTw58+OGHYW+ff/55xOePHz8+UKxYscA555wTuOiiiwK1atVKcZuGDx8eKFSokPuMdd5Lly4dmD17dtjPLtrrpuQzatWqVeC+++6L+B6B1CAAAtIpANq9e3egYMGCgYEDBwa3KVjo2bNn4NixY+7+4cOHAy1btgz+gldgoQvJzz//HHzOkiVLAqtXr3b/1oVCQcBPP/3k7h86dCjQvHnzwDXXXBPc/+KLL3YXnCNHjrj7P/74o3vO1KlTg/voOLqIrV+/3t0/evSou9A/8sgjMbUj2vuIZwB02223Bfd5+eWX3bZ77703uG3w4MEuEPToebqIz5kzJ7ht/vz57rNYsWJFsm3asGGDO+6CBQuC2x577DH3njzHjx8PfPLJJ8m2ddmyZYHcuXMHvvzyS3c/ISEh0KNHj2QDoN69e7vHZfr06YFcuXIFP4fkaP/u3buHvfXv3z/sc//44w/32Y8cOTLRedb5iLVNOkbevHkDEydODB5DAbUCbH3uyZ2PWF431s9IgXjo5wDEC11gQBwtXLjQdQ0o9f/OO++4riOvS2zHjh322Wef2dNPP21TpkxxXUi6qTtr5syZbp+8efO6m7rP6tev77apGy2UulS8Lgp1Pdxzzz123XXX2Xvvvee6Wb766iv7+uuvXdeT17VwySWXuC6n0GJS/dvr0tFrNm/e3JYvXx61HbG8j3h3K3qaNWuW7LZHHnnEjhw5Yvnz57cPPvjAypYt67qZVH+ltom6/L777rsTugNl+/bt7qe6qjzq6tF2db/pc1RNSpcuXZJt4+TJk61u3bp28cUXu/u5cuWy+++//4RuPrntttvc46LuS6+r0fsskoqlJikcnYvKlSvbHXfcEdxWvnx5d4u1TerGU1epuqM8jz76qL322mu2aNEi1x2W0tdNyWekz8T7fIB4IgAC4ki1Lapr2LZtmwuGdPHSL3r5+++/3S/6H374wY2WCaX6CClVqpS98cYb9tBDD7kCUNU9XH/99cELqyQdraP6G138dTHxLhS6YIVSTcr333+faJteK5SCKdXhRGtHLO8jnkKDErUxuW1qjxcAqS7o8OHD7sIdSu+hdOnSyb6GV88SOhz8rrvucgFgzZo13ciwdu3aWd++fd1op6RUB5X0cwk3qir0vCtIVWDlnffkKMjU+wmnWLFi1rFjx2QfW7t2rQsmvOAmnEht0ued9PukwEjnXY8lFwBFe92UfEb6TJKrWQLSigAISKciaP1y79GjhyuGVZZCFypRZsD7Kzs5vXr1shtuuMF+++03l2np2rWrvfXWW26UmezcuTPR/rqvC03oRUzFz6GFwbpfpkyZFL2XcO1o3LhxTO8js+g86yKaXPYlHA3BV8ZHw8G9Ym4dR8dQVm3u3Lk2YsQIa9iwoRs6npTOvQqIQyX9nFLryy+/dEXY4SgYCRcAaVRW0iA1pfS9UcF5qAMHDrgAJtx3KtrrpuQz0mei/0NAvDERIpBOrrjiCjfi5+6773YZCmUSdNMoqKQ2bNjgfip7pL94FdAo6/Dwww+7IMMbKSPq3tIFKLT7pUGDBm4kmLpq9Ne6tnm0ryb3O//882Nue6R2xPI+MpOyVLr4KmgJpSBCwUxylPVo0aJFoud470XZIXU7vvjii647TNmepNR9qBFSGinlmTp1alzej7KIChTC3YYNGxb2uRdddJHLSiYNRpShjJW+Nxrtp6yNR91WJ5100gnds7G+bqyfkf7faDTjhRdeGHN7gViRAQLS0eDBg90wbU0EqHlSlEFRPY4Co06dOrlAQ8GJfsFr8j1dQPWY6i3q1KnjhlTPnj3bBSChE8S1bdvWZWd0EXnzzTddlsCjC7Wef/DgQXcM1SKpy0hdOrGK1o5o7yMz6eJ60003WYcOHVyXlQJCDU/XXEQzZsw4Yai/p3fv3nbffffZs88+67qANCxeNS46njIaqrHSBV/Bn2aMDqVzoMyY3r/qXnT+dN4lWvdTetI5UBZS8ywpEFfGRufh5ptvTjQPUyTKLimg0XtTvZmyiUOHDrUnn3wybAYo2uvG+hnpO6dMU2j9ERAvZICAOFGAk3Q+GgUPQ4YMCf4lrCyKupQ05443T4+Kib2gQfurCFQBi375K9hR7U7r1q0TdU0pENFx9Lj2C/0LWUGJjq0Lr56ri4cyN+riCa21SPrXuyZx9Lq0orUj2vuIRAGEzlW0yfRU96H7hQsXDu6jx7Qt9L3o4qptXtG3F6B99NFHrhhd7VPdjs6B5rgJR1186ppRzY0XvOo8K+uj5yqjp/mTFBwlbavOtYrPe/bs6YImPe5l4bz6leTej2hbcnVF8TJu3DiXrVNGS5+Zarq84CfWNumcqB5MmSDVmql794EHHog4EWKk1431M3r55ZetX79+LtsExBtLYQDZiP5yVkbp+eefz+ym5Ejqkpk+fbo98cQTKX6ugsDQOixd/DWLsbp9NKIOKaPRhqoz02izpBN9AvHA/0oACKnl0S01FOxoFJqmHVBB9KuvvuqySAQ/qaMi6TFjxvDdRLohAAKyEXVdacQSsh6NElPdj7pxdPFWsXpKCs8BZCy6wAAAgO9QBA0AAHyHAAgAAPgOARAAAPAdAiAAAOA7BEAAAMB3CIAAAIDvEAABAADfIQACAAC+QwAEAADMb/4fo4Ug3LBCDrQAAAAASUVORK5CYII=", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "model.restore_traces(point_fit)\n", + "model.sample_posterior_predictive()\n", + "model.plot_predictive();" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## A note on `sampler=\"laplace\"`\n", + "\n", + "`model.sample(sampler=\"laplace\")` is provided by bambi, which gives HSSM no way\n", + "to pass initial values into its internal MAP step — so it always starts from\n", + "PyMC's default point, the one described above. When the gradient is non-finite\n", + "there, the \"approximation\" is a Gaussian centred on the *start point*, reported\n", + "without an error. HSSM warns when you request it.\n", + "\n", + "Prefer `find_MAP()` followed by `sample()`, which is what this tutorial does." + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": ".venv", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.13.11" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/mkdocs.yml b/mkdocs.yml index 3b2f8d7ef..97970ecb2 100644 --- a/mkdocs.yml +++ b/mkdocs.yml @@ -53,6 +53,7 @@ nav: - Set initial values: tutorials/initial_values.ipynb - Use alternative samplers (Bayeux): tutorials/tutorial_bayeux.ipynb - Run variational inference: tutorials/variational_inference.ipynb + - Find MAP and MLE point estimates: tutorials/map_mle.ipynb - Working with results: - Plot posteriors and predictions: tutorials/plotting.ipynb - Posterior predictive plot gallery: tutorials/ppc_gallery.ipynb @@ -84,6 +85,7 @@ nav: - hssm.ModelConfig: api/model_config.md - hssm.Prior: api/prior.md - hssm.Link: api/link.md + - hssm.PointEstimate: api/point_estimate.md - hssm.rl: api/rl.md - Useful functions: - hssm.load_data: api/load_data.md @@ -152,6 +154,7 @@ plugins: - tutorials/ppc_gallery.ipynb - tutorials/cartoon_gallery.ipynb - tutorials/scientific_workflow_hssm.ipynb + - tutorials/map_mle.ipynb - archive/pymc_to_hssm.ipynb - tutorials/tutorial_bayeux.ipynb - tutorials/bayesflow_lre_integration.ipynb diff --git a/src/hssm/__init__.py b/src/hssm/__init__.py index e94c6d1f0..5769a2221 100644 --- a/src/hssm/__init__.py +++ b/src/hssm/__init__.py @@ -18,6 +18,7 @@ from .hssm import HSSM from .link import Link from .modelconfig import list_models +from .optimize import PointEstimate from .param import UserParam as Param from .prior import Prior from .register import register_model @@ -42,6 +43,7 @@ "load_data", "ModelConfig", "Param", + "PointEstimate", "Prior", "check_data_for_rl", "register_model", diff --git a/src/hssm/base.py b/src/hssm/base.py index 9d5b509dd..1f72d4601 100644 --- a/src/hssm/base.py +++ b/src/hssm/base.py @@ -55,7 +55,7 @@ _split_array, ) -from . import plotting +from . import optimize, plotting from .config import BaseModelConfig from .modelconfig import list_models from .param import Params @@ -284,7 +284,8 @@ def __init__( self._inference_obj: DataTree | None = None self._inference_obj_vi: DataTree | Approximation | None = None self._vi_approx = None - self._map_dict = None + self._map_dict: optimize.PointEstimate | dict | None = None + self._mle_dict: optimize.PointEstimate | dict | None = None # endregion # ===== Initial Values Configuration ===== @@ -548,16 +549,532 @@ def _store_init_args( result.update(extra_kwargs) return result - def find_MAP(self, **kwargs): + def _point_estimate_setup( + self, + start: dict | None, + method: str | None, + kind: str, + objective_factory: Callable[[], Any] | None = None, + probe_gradient: bool = True, + ) -> tuple[dict, str, Any]: + """Prepare the shared preconditions of ``find_MAP`` and ``find_MLE``. + + Runs the RL/aDDM extra-field refresh, resolves the start point, applies + the ``float32`` safety net and settles which optimizer will really run: + both `pymc.find_MAP` and `optimize_objective` swap in the gradient-free + method behind the caller's back when the likelihood has no gradient, so + the swap has to be resolved here or the estimate would record a method + that never ran. + + The objective is built by ``objective_factory`` rather than passed in, + so that it is constructed *after* the extra-field refresh. Callers that + do not need one (``find_MAP`` differentiates the joint log-density, + which `resolve_gradient_method` builds itself) pass ``None``. + + ``probe_gradient=False`` returns the method with only the ``float32`` + rule applied, leaving the gradient question to the caller. `find_MLE` + uses it: probing here would build the gradient graph that + `optimize.compile_objective` is about to build anyway, and that graph is + the expensive part of both calls. + + Notes + ----- + The extra-field refresh mirrors what ``sample()`` does at the top of a + run, and is kept for parity with it, but note that it currently has no + effect on the likelihood: ``_update_extra_fields`` rebinds an + ``extra_fields`` attribute on the distribution *class*, while the logp + graph reads the ``extra_fields`` list closed over in + ``make_distribution``. The ordering above is therefore a precaution + against that plumbing being repaired, not a live dependency. + """ + if self._check_extra_fields(): + self._update_extra_fields() + + objective = None if objective_factory is None else objective_factory() + resolved_start = self.initvals if start is None else start + resolved_method = optimize.resolve_method(method, kind) + if probe_gradient: + resolved_method = optimize.resolve_gradient_method( + self.pymc_model, resolved_method, objective + ) + return resolved_start, resolved_method, objective + + def find_MAP( + self, + start: dict | None = None, + method: str | None = None, + n_starts: int = 1, + strict: bool = False, + se: bool = False, + seed: int | None = None, + return_raw: bool = False, + **kwargs, + ) -> "optimize.PointEstimate | tuple | None": """Perform Maximum A Posteriori estimation. + Maximizes the joint log-density of the model (likelihood + priors) with + respect to all continuous parameters, delegating the optimization itself + to `pymc.find_MAP`. + + Parameters + ---------- + start : optional + The starting point, as a dictionary of constrained parameter values. + Defaults to `model.initvals` --- HSSM's own processed initial values. + This matters: PyMC's default start (`t=2.0`, `a=2.0`) puts many + models where the gradient of the log-density is non-finite, from + which the optimizer cannot move. + method : optional + Any `scipy.optimize.minimize` method. Defaults to `"L-BFGS-B"`, or + to the gradient-free `"Powell"` when the model runs at `float32` + precision. `"Powell"` is also substituted, whatever was asked for, + when the likelihood exposes no gradient (`blackbox`); the method + that actually ran is recorded on the returned estimate. + n_starts : optional + Number of starting points to try. Values above 1 add jittered copies + of `start` and keep whichever run reaches the highest log-density. + Must be at least 1. Defaults to 1. + strict : optional + If True, raise a `RuntimeError` instead of warning when the + optimizer fails to converge. Defaults to False. + se : optional + If True, also compute the inverse Hessian at the estimate. For a MAP + these are posterior standard deviations under a Laplace (normal) + approximation around the mode, not standard errors in the sampling + sense --- see `find_MLE` for those. Not available for `blackbox` + likelihoods. Costs `2 * n_parameters` gradient evaluations plus an + `n x n` inverse, which is noticeable on large hierarchies. Defaults + to False. + seed : optional + Seed for the jitter applied when `n_starts > 1`, making multi-start + runs reproducible. + return_raw : optional + If True, return a `(point, OptimizeResult)` tuple as + `pymc.find_MAP` does. The raw result is always available as + `model.map.opt_result`. + kwargs + Other arguments passed to `pymc.find_MAP` (`vars`, `maxeval`, + `progressbar`, ...). `include_transformed` is the exception: + HSSM consumes it rather than forwarding it, because scoring and the + convergence audit both need the transformed names. PyMC is always + asked for them and they are dropped afterwards, so the result is + the same as PyMC's. + Returns ------- - dict - A dictionary containing the MAP estimates of the model parameters. + PointEstimate | None + A `dict` subclass holding the MAP estimate along with the optimizer + metadata, or `None` if the optimization failed and `strict` is + False. A failed estimate is never cached on the model. + + Raises + ------ + ValueError + If `n_starts` is less than 1. + + Notes + ----- + On a hierarchical model this is the mode of the *joint* posterior over + parameters and group-level offsets, not of the marginal posterior of the + population-level parameters. + + Build-time jitter (`initval_jitter`) draws from the unseeded global + `numpy` random state, so for models with vector-valued (hierarchical) + parameters the default start --- and hence the MAP --- varies slightly + between model instances. Pass an explicit `start=` for full + reproducibility. """ - self._map_dict = pm.find_MAP(model=self.pymc_model, **kwargs) - return self._map_dict + optimize.validate_n_starts(n_starts) + resolved_start, resolved_method, _ = self._point_estimate_setup( + start, method, "MAP" + ) + rng = np.random.default_rng(seed) + # A concrete seed even when the caller passed none: the expansion below + # and PyMC's own would otherwise each draw fresh entropy, and on a model + # with a stochastic initval strategy the audit would end up comparing + # against a start the optimizer never saw. + run_seed = seed if seed is not None else int(rng.integers(2**32)) + # Scoring and the failure audit both need the transformed value-variable + # names, so always ask PyMC for them and drop them afterwards if the + # caller asked for a constrained-only point. + include_transformed = kwargs.pop("include_transformed", True) + # No scorer is compiled up front. `pm.find_MAP` already reports the + # joint log-density at the point it found, as `opt_result.fun` negated, + # and that agrees with a freshly compiled scorer bit-for-bit under + # L-BFGS-B, Powell, Nelder-Mead and BFGS. Compiling a second evaluator + # for a number PyMC hands back costs a full logp compile --- about a + # tenth of a `find_MAP` on an analytical DDM, and considerably more on + # an ONNX/JAX likelihood. `scorer` below is the fallback for a SciPy + # method that leaves `fun` unset, built at most once and usually never. + scorer: Callable[[dict], float] | None = None + + def candidate_score(point: dict, opt_result: Any) -> float: + """Return the joint log-density at a converged candidate.""" + nonlocal scorer + reported = optimize.objective_value(opt_result) + if reported is not None: + return -reported + if scorer is None: + scorer = optimize.make_scorer(self.pymc_model) + return scorer(point) + + # The (possibly partial) start is expanded once, here, rather than once + # per candidate inside the loop. That gives the audit the *complete* + # start it needs to tell "the optimizer never moved" from "one parameter + # happened not to move", and it is also the point the jitter is applied + # to: transformed space, where no jittered candidate can leave the + # model's support. + base_start = optimize.make_start_point( + self.pymc_model, resolved_start, seed=run_seed, validate=False + ) + + def candidate_overrides(index: int, candidate: dict) -> dict: + """Return the ``start=`` override for this candidate. + + A jittered candidate is handed over in transformed space as it + stands: `pm.find_MAP` routes `start=` through + `convert_str_to_rv_dict`, which recognizes transformed names and + back-transforms them itself. Expanding to constrained space here + instead would buy a compiled pytensor function and, on a regression + model, an evaluation of every trial-wise deterministic per + candidate --- all of it discarded by the same call. + """ + if index == 0: + # The first candidate is the unjittered start, which PyMC + # expands itself from `resolved_start` under the same seed --- + # the very point `base_start` already holds. + return resolved_start + return candidate + + best_point: dict | None = None + best_result: Any = None + best_score = -np.inf + n_converged = 0 + failures: list[list[str]] = [] + + for index, candidate_start in enumerate( + optimize.jittered_starts(base_start, n_starts, rng) + ): + try: + point, opt_result = pm.find_MAP( + start=candidate_overrides(index, candidate_start), + method=resolved_method, + model=self.pymc_model, + return_raw=True, + seed=run_seed, + **kwargs, + ) + # Only `SamplingError` means "PyMC rejected this start". Catching + # `ValueError` here too would relabel an unknown `method=` or a bad + # `options=` -- which SciPy raises from inside this call -- as a bad + # starting point, and turn a usage error into a convergence warning. + except pm.exceptions.SamplingError as error: + failures.append([f"the starting point was rejected by PyMC ({error})"]) + continue + + candidate_reasons = optimize.audit_result( + opt_result, point, candidate_start + ) + if candidate_reasons: + failures.append(candidate_reasons) + continue + + n_converged += 1 + score = candidate_score(point, opt_result) + if score > best_score: + best_point, best_result, best_score = point, opt_result, score + + if best_point is None: + # Drop any estimate an *earlier* successful call cached. Without + # this the "a failed estimate is never cached" contract holds only + # on a first run: `model.map` would keep returning the stale point + # instead of raising, and `sample(initvals="map")` would slip past + # its `_map_dict is None` guard and initialize from it. + self._map_dict = None + optimize.report_failure("MAP", failures, strict) + return None + + if n_starts > 1: + _logger.info("find_MAP: %d of %d starts converged.", n_converged, n_starts) + + estimate = self._make_point_estimate( + best_point, + best_result, + kind="MAP", + method=resolved_method, + logp=best_score, + n_starts=n_starts, + n_converged=n_converged, + se=se, + include_transformed=include_transformed, + ) + self._map_dict = estimate + + if return_raw: + return estimate, best_result + return estimate + + def find_MLE( + self, + start: dict | None = None, + method: str | None = None, + n_starts: int = 1, + strict: bool = False, + se: bool = False, + seed: int | None = None, + return_raw: bool = False, + allow_unidentified: bool = False, + maxeval: int = 5000, + progressbar: bool = True, + **kwargs, + ) -> "optimize.PointEstimate | tuple | None": + """Perform Maximum Likelihood estimation. + + Maximizes the *observed* log-likelihood --- the priors are dropped --- so + the result is a frequentist point estimate rather than a posterior mode. + + Parameters + ---------- + start : optional + The starting point, as a dictionary of constrained parameter values. + Defaults to `model.initvals`. + method : optional + Any `scipy.optimize.minimize` method. Defaults to `"L-BFGS-B"`, or + to the gradient-free `"Powell"` at `float32` precision and for + `blackbox` likelihoods. The method that actually ran is recorded on + the returned estimate. + n_starts : optional + Number of starting points to try. Values above 1 add jittered copies + of `start` and keep the highest-likelihood run. Must be at least 1. + Defaults to 1. + strict : optional + If True, raise a `RuntimeError` instead of warning when the + optimizer fails to converge. Defaults to False. + se : optional + If True, also compute standard errors from the observed information + (the inverse Hessian of the observed log-likelihood). Costs + `2 * n_parameters` gradient evaluations plus an `n x n` inverse. + Defaults to False. + seed : optional + Seed for the jitter applied when `n_starts > 1`. + return_raw : optional + If True, return a `(point, OptimizeResult)` tuple. + allow_unidentified : optional + Bypass the hierarchical-model check. Defaults to False --- see Raises. + maxeval : optional + Maximum number of likelihood evaluations. Defaults to 5000. + Exhausting it is reported as a failure naming the budget, not as a + bad starting point. + progressbar : optional + Whether to display a progress bar, as `find_MAP` does. Defaults to + True. + kwargs + Other arguments passed to `scipy.optimize.minimize`. + + Returns + ------- + PointEstimate | None + A `dict` subclass holding the MLE along with the optimizer metadata, + or `None` if the optimization failed and `strict` is False. A failed + estimate is never cached on the model. + + Raises + ------ + ValueError + If the model has group-specific (hierarchical) terms, unless + `allow_unidentified=True`. Dropping the priors leaves the + group-level scale unidentified either way: under the non-centered + parameterization (HSSM's default) the group effect is + `offset * sigma`, so the likelihood is invariant under + `(offset * c, sigma / c)`; under the centered parameterization + `sigma` does not enter the likelihood at all and its gradient is + exactly zero. Use `find_MAP` for hierarchical models. + ValueError + If the model contains `pm.Potential` terms, which + `model.observedlogp` would silently drop. + ValueError + If `n_starts` is less than 1. + + Notes + ----- + The optimization runs in transformed space. The transforms are + bijective, so no Jacobian correction is needed --- a monotone + reparameterization preserves the argmax --- and the prior-derived + interval transforms usefully keep the parameters inside the range the + likelihood network was trained on. The flip side is that this is a + *bounded* MLE, which may sit on a boundary. + """ + optimize.validate_n_starts(n_starts) + terms = optimize.group_specific_terms(self) + if terms and not allow_unidentified: + raise ValueError( + "find_MLE is not well posed for hierarchical models. This model " + f"has group-specific terms {terms}. Dropping the priors leaves " + "the group-level scale unidentified: under the non-centered " + "parameterization (HSSM's default) the likelihood is invariant " + "under rescaling the offsets against sigma, and under the " + "centered parameterization sigma does not enter the likelihood " + "at all. Either way the optimum is a flat ridge rather than a " + "point. Use `find_MAP()` instead, which is well behaved for " + "these models, or pass `allow_unidentified=True` to proceed " + "anyway." + ) + + if self.pymc_model.potentials: + raise ValueError( + "find_MLE cannot be used on a model with `pm.Potential` terms: " + "`model.observedlogp`, the objective it maximizes, silently " + "excludes them, so the reported optimum would not correspond to " + "the model you specified. Use `find_MAP()` instead." + ) + + resolved_start, resolved_method, objective = self._point_estimate_setup( + start, + method, + "MLE", + objective_factory=lambda: self.pymc_model.observedlogp, + # Resolved from `compile_objective` below instead: probing here + # would build the gradient graph that call is about to build anyway. + probe_gradient=False, + ) + rng = np.random.default_rng(seed) + run_seed = seed if seed is not None else int(rng.integers(2**32)) + # Compiled once, outside the loop: the objective and its gradient do not + # depend on which start the optimizer is given, and recompiling them per + # candidate is the dominant cost of a multi-start run. + compiled = optimize.compile_objective(self.pymc_model, objective) + # Whether the gradient could be built is exactly what `compile_objective` + # just found out, so the method follows from its result. + resolved_method = optimize.method_for_gradient( + resolved_method, compiled[1] is not None + ) + # The scorer evaluates the same objective the optimizer maximizes, so it + # reuses that compiled function rather than compiling `observedlogp` a + # second time. + score_candidate = optimize.make_scorer( + self.pymc_model, observed_only=True, function=compiled[0] + ) + expand = optimize.make_point_expander(self.pymc_model) + # Expanded once and jittered in transformed space, so that no candidate + # can be pushed out of the model's support. See `find_MAP`. + base_start = optimize.make_start_point( + self.pymc_model, resolved_start, seed=run_seed, validate=False + ) + + best_point: dict | None = None + best_result: Any = None + best_score = -np.inf + best_method = resolved_method + n_converged = 0 + failures: list[list[str]] = [] + + for candidate_start in optimize.jittered_starts(base_start, n_starts, rng): + # Still checked per candidate: the jitter cannot leave the support, + # but it can land where the likelihood is -inf anyway -- a `t` above + # the fastest response time, for instance. + try: + self.pymc_model.check_start_vals(candidate_start) + except (ValueError, pm.exceptions.SamplingError) as error: + failures.append([f"the starting point was rejected by PyMC ({error})"]) + continue + + point, opt_result, used_method = optimize.optimize_objective( + self.pymc_model, + objective, + candidate_start, + method=resolved_method, + maxeval=maxeval, + progressbar=progressbar, + compiled=compiled, + expand=expand, + **kwargs, + ) + + candidate_reasons = optimize.audit_result( + opt_result, point, candidate_start + ) + if candidate_reasons: + failures.append(candidate_reasons) + continue + + n_converged += 1 + score = score_candidate(point) + if score > best_score: + best_point, best_result, best_score = point, opt_result, score + best_method = used_method + + if best_point is None: + # See the matching comment in `find_MAP`: a failed run must not + # leave a previous run's estimate reachable through `model.mle`. + self._mle_dict = None + optimize.report_failure("MLE", failures, strict) + return None + + if n_starts > 1: + _logger.info("find_MLE: %d of %d starts converged.", n_converged, n_starts) + + estimate = self._make_point_estimate( + best_point, + best_result, + kind="MLE", + method=best_method, + logp=best_score, + n_starts=n_starts, + n_converged=n_converged, + se=se, + observed_only=True, + ) + + self._mle_dict = estimate + + if return_raw: + return estimate, best_result + return estimate + + def _make_point_estimate( + self, + point: dict, + opt_result: Any, + *, + kind: Literal["MAP", "MLE"], + method: str, + logp: float, + n_starts: int, + n_converged: int, + se: bool, + observed_only: bool = False, + include_transformed: bool = True, + ) -> "optimize.PointEstimate": + """Wrap a converged optimizer point in a `PointEstimate`.""" + params = optimize.constrained_params(self.pymc_model, point) + dims, coords = optimize.dims_and_coords(self.pymc_model, list(params)) + errors = ( + optimize.standard_errors( + self.pymc_model, point, observed_only=observed_only + ) + if se + else None + ) + mapping = ( + point + if include_transformed + else optimize.drop_transformed(self.pymc_model, point) + ) + return optimize.PointEstimate( + mapping, + kind=kind, + params=params, + success=True, + message=str(getattr(opt_result, "message", "")), + logp=logp, + method=method, + n_starts=n_starts, + n_converged=n_converged, + se=errors, + opt_result=opt_result, + dims=dims, + coords=coords, + ) def sample( self, @@ -640,9 +1157,27 @@ def sample( "Running map estimation first..." ) self.find_MAP() - kwargs["initvals"] = self._map_dict - else: - kwargs["initvals"] = self._map_dict + # A failed `find_MAP` caches nothing, so `_map_dict` is + # still None here. Falling through would hand + # `initvals=None` to bambi and start NUTS from PyMC's + # default point -- silently sampling from somewhere the + # user never asked for. Fail loudly instead. + if self._map_dict is None: + raise RuntimeError( + "initvals='map' was requested but MAP " + "estimation did not converge, so there is no " + "MAP estimate to initialize from. Investigate " + "with `model.find_MAP(strict=True)`, or pass " + "`initvals=None` to use the default initial " + "values." + ) + # Hand over only the constrained free parameters. The + # raw optimizer point also carries transformed names and + # trial-wise deterministics; PyMC tolerates both, but + # there is no reason to push them through. + kwargs["initvals"] = dict( + getattr(self._map_dict, "params", self._map_dict) + ) else: raise ValueError( "initvals argument must be a dictionary or 'map'" @@ -661,6 +1196,24 @@ def sample( else: sampler = "pymc" + if sampler == "laplace": + # bambi's `_run_laplace` calls a bare `pm.find_MAP()` and takes the + # Hessian at whatever it returns. It offers no way to pass initial + # values, and HSSM resets `rvs_to_initial_values` to None, so the + # optimization always starts from PyMC's default point -- the same + # start that leaves the gradient non-finite for many HSSM models. + # When that happens the "approximation" is a Gaussian centred on the + # start point, reported without complaint. + _logger.warning( + "The 'laplace' sampler is provided by bambi, which gives HSSM " + "no way to supply initial values, so its internal MAP step " + "starts from PyMC's default point. If the gradient is " + "non-finite there, the result is a Gaussian around that start " + "point rather than around the posterior mode -- with no error " + "raised. Prefer `model.find_MAP()` (which uses HSSM's own " + "initial values) followed by `model.sample()`." + ) + if self.loglik_kind == "blackbox": if sampler in ["blackjax", "numpyro", "nutpie"]: raise ValueError( @@ -1861,24 +2414,49 @@ def vi_approx(self) -> Approximation: return self._vi_approx @property - def map(self) -> dict: - """Return the MAP estimates of the model parameters. + def map(self) -> "optimize.PointEstimate | dict": + """Return the MAP estimate of the model parameters. Raises ------ ValueError - If the model has not been sampled yet. + If `find_MAP` has not been run, or if it ran but did not converge + (a failed estimate is never cached). Returns ------- - dict - A dictionary containing the MAP estimates of the model parameters. + PointEstimate + A `dict` subclass containing the MAP estimate of the model + parameters, along with the optimizer metadata. """ - if not self._map_dict: - raise ValueError("Please compute map first.") + # `is None` rather than a falsy check: None is the only "not computed" + # sentinel, and a computed-but-empty estimate must not read as missing. + if self._map_dict is None: + raise ValueError("Please compute map first, by running `find_MAP()`.") return self._map_dict + @property + def mle(self) -> "optimize.PointEstimate | dict": + """Return the maximum likelihood estimate of the model parameters. + + Raises + ------ + ValueError + If `find_MLE` has not been run, or if it ran but did not converge + (a failed estimate is never cached). + + Returns + ------- + PointEstimate + A `dict` subclass containing the MLE of the model parameters, along + with the optimizer metadata. + """ + if self._mle_dict is None: + raise ValueError("Please compute mle first, by running `find_MLE()`.") + + return self._mle_dict + @property def initvals(self) -> dict: """Return the initial values of the model parameters for sampling. diff --git a/src/hssm/optimize.py b/src/hssm/optimize.py new file mode 100644 index 000000000..d66d34e73 --- /dev/null +++ b/src/hssm/optimize.py @@ -0,0 +1,1366 @@ +"""Point estimation (MAP / MLE) for HSSM models. + +This module holds the machinery shared by ``hssm.HSSM.find_MAP`` and +``hssm.HSSM.find_MLE``: + +* ``PointEstimate`` --- the ``dict`` subclass both methods return. It keeps + the plain-mapping behaviour older code relies on (``sample(initvals="map")`` + checks ``isinstance(initvals, dict)``) while carrying optimizer metadata and + ArviZ-friendly exporters. +* the start-value, failure-auditing, ``float32`` and standard-error helpers. +* ``optimize_objective``, a small ``scipy.optimize.minimize`` driver used by + ``find_MLE``. ``find_MAP`` keeps delegating to ``pymc.find_MAP``, which + already implements the same loop for the *joint* log-density; only the + objective differs for MLE, and PyMC offers no hook to swap it. + +Of these, only ``PointEstimate`` (re-exported as ``hssm.PointEstimate``) and +``score_point`` are meant for users; see ``__all__``. The rest is implementation +shared with ``hssm.base``. + +Notes +----- +Point estimation is gradient-based and therefore precision sensitive. Under +``hssm.set_floatX("float32")`` the gradient noise sits above L-BFGS-B's default +tolerances, so the optimizer stops early *and reports success*. The helpers +below warn unconditionally in that configuration and fall back to the +gradient-free ``"Powell"`` method unless the caller picked one explicitly. +""" + +from __future__ import annotations + +import logging +import warnings +from typing import TYPE_CHECKING, Any, Callable, Literal, cast + +import arviz as az +import numpy as np +import pandas as pd +import pytensor +import pytensor.gradient as tg +import pytensor.tensor as pt +from pymc.blocking import DictToArrayBijection, RaveledVars +from pymc.initial_point import make_initial_point_fn +from pymc.model.transform.conditioning import remove_value_transforms +from pymc.progress_bar import CustomProgress, default_progress_theme +from pymc.pytensorf import rewrite_pregrad +from pymc.util import get_default_varnames +from rich.console import Console +from rich.progress import Progress, TextColumn +from scipy import optimize as sp_optimize + +if TYPE_CHECKING: + from pymc import Model as PyMCModel + from xarray import DataTree + +_logger = logging.getLogger("hssm") + +#: The supported surface of this module. Everything else defined here is +#: machinery shared with ``hssm.base`` --- importable, but not API: it may +#: change shape whenever the two ``find_*`` methods need it to. +__all__ = [ + "PointEstimate", + "score_point", +] + +#: Exceptions PyTensor raises when a graph has no usable gradient. Mirrors the +#: tuple ``pymc.find_MAP`` catches, so blackbox likelihoods take the same +#: gradient-free path here as they do there. +NO_GRADIENT_ERRORS = (AttributeError, NotImplementedError, tg.NullTypeGradError) + +#: PyMC's own default. Kept as a literal because HSSM uses ``method=None`` as +#: its own "user did not choose" sentinel and must never forward that ``None`` +#: to ``scipy.optimize.minimize`` (where it silently auto-selects a *gradient* +#: method while ``jac=True`` is still set). +DEFAULT_METHOD = "L-BFGS-B" + +#: Must match ``pymc.find_MAP``'s case-sensitive ``method != "Powell"`` guard. +GRADIENT_FREE_METHOD = "Powell" + +_FLOAT32_WARNING = ( + "Point estimation is running with `floatX == 'float32'`. Gradient noise at " + "single precision exceeds the optimizer's default tolerances, so the " + "optimizer can stop far from the optimum *and still report success*. The " + "returned estimate should be treated as approximate. Call " + "`hssm.set_floatX('float64')` before building the model for reliable " + "results." +) + +#: Rows of ``PointEstimate.to_dataframe`` shown in ``repr``. A hierarchical +#: model has one row per group-level offset, and an unbounded ``to_string()`` +#: would dump hundreds of them into the REPL. +REPR_MAX_ROWS = 20 + +#: Above this many scalar parameters, ``standard_errors`` logs what it is +#: about to do: the finite-difference Hessian costs ``2 * size`` full gradient +#: evaluations plus a ``size x size`` inverse, which can take minutes. +_SE_SIZE_HINT = 100 + + +# region ===== PointEstimate ===== +def _rebuild_point_estimate(mapping: dict, state: dict) -> "PointEstimate": + """Reconstruct a ``PointEstimate`` during unpickling.""" + obj = PointEstimate.__new__(PointEstimate) + dict.__init__(obj, mapping) + obj.__dict__.update(state) + return obj + + +class PointEstimate(dict): + """A point estimate of the model parameters, with optimizer metadata. + + This is a ``dict`` subclass, so it can be passed anywhere a plain point + dictionary was accepted before (notably ``model.sample(initvals=...)``). The + mapping itself is what the optimizer returned --- constrained values, + transformed value-variable names (``t_log__``, ``z_interval__``) and + deterministics --- while ``params`` holds only the constrained free + parameters. The transformed entries are absent when the estimate was + produced with ``include_transformed=False``; ``params`` is unaffected. + + Attributes + ---------- + kind + ``"MAP"`` or ``"MLE"``. + params + The constrained free parameters only, with transformed duplicates and + deterministics dropped. + success + Always ``True``. A run that did not converge never becomes a + ``PointEstimate`` --- ``find_MAP``/``find_MLE`` warn and return ``None`` + instead --- so this exists for the shape of the record, not as a check + worth branching on. + message + The optimizer's termination message. + logp + The value of the maximized objective at the returned point --- the joint + log-density for ``"MAP"``, the observed log-likelihood for ``"MLE"``. + method + The ``scipy.optimize.minimize`` method that actually produced the + estimate. This is not necessarily the one that was requested: a + likelihood without a gradient (``blackbox``) forces the gradient-free + fallback, and so does ``float32`` precision. + n_starts + How many starting points were tried. + n_converged + How many of those starts converged. + se + Square roots of the diagonal of the inverse negative Hessian at the + estimate, keyed like ``params``, or ``None`` if they were not + requested (or could not be computed). For ``"MLE"`` these are standard + errors in the usual sense --- the observed-information approximation to + the sampling distribution. For ``"MAP"`` they are posterior standard + deviations under a Laplace (normal) approximation around the mode, which + is a different object that happens to be computed the same way. + opt_result + The raw ``scipy.optimize.OptimizeResult``. Never ``None`` on an + estimate that was returned: a missing result is one of the conditions + ``audit_result`` fails a start on, so a run that exhausted ``maxeval`` + or was interrupted yields ``None`` from ``find_MAP``/``find_MLE`` + rather than an estimate carrying no result. + + Notes + ----- + ``dict.copy`` on a ``dict`` subclass returns a plain ``dict``, which would + silently drop every attribute above. ``copy`` is overridden to preserve + them; ``pickle`` and ``cloudpickle`` are handled by ``__reduce__``. + """ + + _META_FIELDS = ( + "kind", + "params", + "success", + "message", + "logp", + "method", + "n_starts", + "n_converged", + "se", + "opt_result", + "dims", + "coords", + ) + + def __init__( + self, + point: dict[str, Any], + *, + kind: Literal["MAP", "MLE"] = "MAP", + params: dict[str, Any] | None = None, + success: bool = True, + message: str = "", + logp: float = np.nan, + method: str | None = None, + n_starts: int = 1, + n_converged: int = 1, + se: dict[str, Any] | None = None, + opt_result: Any = None, + dims: dict[str, list[str]] | None = None, + coords: dict[str, list[Any]] | None = None, + ): + super().__init__(point) + self.kind = kind + self.params = dict(point) if params is None else dict(params) + self.success = success + self.message = message + self.logp = logp + self.method = method + self.n_starts = n_starts + self.n_converged = n_converged + self.se = None if se is None else dict(se) + self.opt_result = opt_result + self.dims = dict(dims) if dims else {} + self.coords = dict(coords) if coords else {} + + def __reduce__(self): + """Preserve the metadata attributes through ``pickle``/``cloudpickle``.""" + state = {field: getattr(self, field) for field in self._META_FIELDS} + return (_rebuild_point_estimate, (dict(self), state)) + + def copy(self) -> "PointEstimate": + """Return a shallow copy that keeps the optimizer metadata. + + ``dict.copy`` would return a plain ``dict``, silently dropping + ``params``, ``se`` and everything else that makes this more than + a mapping. + + Returns + ------- + PointEstimate + A new estimate with the same mapping and metadata. + """ + # Every entry of ``_META_FIELDS`` is also an ``__init__`` keyword, so the + # constructor is the single place that knows how to rebuild the object. + return PointEstimate( + dict(self), **{field: getattr(self, field) for field in self._META_FIELDS} + ) + + __copy__ = copy + + def __repr__(self) -> str: # noqa: D105 + status = "converged" if self.success else "DID NOT CONVERGE" + header = f"{self.kind} estimate ({self.method}, {status}, logp={self.logp:.4g})" + frame = self.to_dataframe() + if len(frame) > REPR_MAX_ROWS: + hidden = len(frame) - REPR_MAX_ROWS + body = ( + f"{frame.head(REPR_MAX_ROWS).to_string()}\n" + f"... and {hidden} more row(s); use `.to_dataframe()` for all." + ) + else: + body = frame.to_string() + return f"{header}\n{body}" + + def to_dataframe(self) -> pd.DataFrame: + """Return the estimate as a tidy DataFrame, one row per scalar entry. + + Vector-valued parameters are expanded to ``name[i]`` rows. An ``se`` + column is included only when standard errors were computed. + + Returns + ------- + pd.DataFrame + Indexed by parameter name, with an ``estimate`` column and, + optionally, an ``se`` column. + """ + names: list[str] = [] + estimates: list[float] = [] + errors: list[float] = [] + + for name, value in self.params.items(): + array = np.atleast_1d(np.asarray(value, dtype=float)) + se_array = ( + np.atleast_1d(np.asarray(self.se[name], dtype=float)).ravel() + if self.se is not None and name in self.se + else None + ) + flat = array.ravel() + scalar = np.ndim(value) == 0 or flat.size == 1 + for i, entry in enumerate(flat): + names.append(name if scalar else f"{name}[{i}]") + estimates.append(float(entry)) + errors.append( + float(se_array[i]) + if se_array is not None and i < se_array.size + else np.nan + ) + + frame = pd.DataFrame( + {"estimate": estimates, "se": errors}, + index=pd.Index(names, name="parameter"), + ) + if frame["se"].isna().all(): + frame = frame.drop(columns=["se"]) + return frame + + def to_datatree(self) -> "DataTree": + """Return the estimate as a 1-chain, 1-draw posterior ``DataTree``. + + This makes the point estimate usable with the rest of the ArviZ/HSSM + toolchain: ``az.summary``, ``model.restore_traces``, + ``model.sample_posterior_predictive`` and the plotting functions all + accept the result. + + Returns + ------- + DataTree + A ``DataTree`` with a single ``posterior`` group of shape + ``(chain=1, draw=1, ...)``. + + Notes + ----- + With a single draw the posterior summaries degenerate: ``az.summary`` + reports a standard deviation of ``0`` and ``NaN`` for ``ess_bulk``, + ``ess_tail`` and ``r_hat``. That is expected --- there is no sampling + distribution to summarize --- not a sign that something went wrong. + """ + data = { + name: np.asarray(value)[np.newaxis, np.newaxis, ...] + for name, value in self.params.items() + } + dims = {name: dim for name, dim in self.dims.items() if name in data} + used = {dim for entry in dims.values() for dim in entry} + coords = {key: value for key, value in self.coords.items() if key in used} + dataset = az.dict_to_dataset( + data, coords=cast("Any", coords), dims=cast("Any", dims) + ) + return az.convert_to_datatree(dataset) + + +# endregion + + +# region ===== shared helpers ===== +def resolve_method(method: str | None, kind: str) -> str: + """Resolve the optimizer method, applying the ``float32`` safety net. + + Parameters + ---------- + method + The user's choice, or ``None`` if they did not make one. HSSM keeps + ``None`` as its own sentinel because PyMC's default is the *string* + ``"L-BFGS-B"``, which makes "user chose L-BFGS-B" indistinguishable from + "user chose nothing". + kind + ``"MAP"`` or ``"MLE"``, used in the warning text. + + Returns + ------- + str + The method to use. + """ + if pytensor.config.floatX != "float32": + return method or DEFAULT_METHOD + + fallback = method is None + warnings.warn( + f"find_{kind} was called with float32 precision. {_FLOAT32_WARNING}" + + ( + f" Falling back to the gradient-free '{GRADIENT_FREE_METHOD}' " + "method, which is more robust at this precision; pass an explicit " + "`method=` to override." + if fallback + else "" + ), + UserWarning, + # resolve_method <- _point_estimate_setup <- find_MAP/find_MLE <- user + stacklevel=4, + ) + # Narrow on `method is None` rather than `fallback`: mypy does not carry the + # narrowing through the intermediate boolean. + return GRADIENT_FREE_METHOD if method is None else method + + +def gradient_available(pymc_model: "PyMCModel", objective: Any = None) -> bool: + """Check whether a symbolic gradient of ``objective`` can be built. + + Only the *graph* is built, not a compiled function, and it raises exactly + where ``pymc.find_MAP`` does --- which is what makes it a faithful predictor + of whether the optimizer will be able to use gradients, rather than a guess + from ``loglik_kind``. + + It is not free: building the log-density and its gradient graph costs on the + order of a tenth of a second for an ONNX/JAX model, a few percent of the + optimization that follows. That is the price of recording the method that + really ran instead of the one that was asked for. Callers that already know + the answer skip it by passing the gradient-free method explicitly. + + Parameters + ---------- + pymc_model + The model whose continuous value variables the gradient is taken over. + objective + The objective to differentiate. Defaults to the joint log-density, the + objective ``find_MAP`` maximizes. + + Returns + ------- + bool + ``True`` unless the likelihood exposes no gradient (``blackbox``). + """ + objective = pymc_model.logp(jacobian=False) if objective is None else objective + try: + flat_gradient(objective, pymc_model.continuous_value_vars) + except NO_GRADIENT_ERRORS: + return False + return True + + +def method_for_gradient(method: str, has_gradient: bool) -> str: + """Return the method the optimizer will *actually* run. + + Split from ``resolve_gradient_method`` so that a caller which has *already* + established gradient availability --- ``find_MLE`` learns it from + ``compile_objective``, which has to build the gradient anyway --- can settle + the method without building the gradient graph a second time. + + Parameters + ---------- + method + The method resolved so far (never ``None``). + has_gradient + Whether a gradient of the objective could be built. + + Returns + ------- + str + ``method``, or ``GRADIENT_FREE_METHOD`` if no gradient is available. + """ + if method == GRADIENT_FREE_METHOD or has_gradient: + return method + + _logger.warning( + "No gradient is available for this likelihood, so '%s' cannot be used. " + "Falling back to the gradient-free '%s' method.", + method, + GRADIENT_FREE_METHOD, + ) + return GRADIENT_FREE_METHOD + + +def resolve_gradient_method( + pymc_model: "PyMCModel", method: str, objective: Any = None +) -> str: + """Return the method the optimizer will *actually* run, probing the gradient. + + Both ``pymc.find_MAP`` and ``optimize_objective`` silently swap in + the gradient-free method when the likelihood has no gradient, but neither + reports that it did --- so without this the estimate's recorded ``method`` + would name an optimizer that never ran. + + Parameters + ---------- + pymc_model + The model being optimized. + method + The method resolved so far (never ``None``). + objective + The objective that will be maximized. Defaults to the joint + log-density. + + Returns + ------- + str + ``method``, or ``GRADIENT_FREE_METHOD`` if no gradient is available. + """ + # Short-circuited before the probe: the answer cannot change the method, and + # `gradient_available` is the expensive part of this call. + if method == GRADIENT_FREE_METHOD: + return method + return method_for_gradient(method, gradient_available(pymc_model, objective)) + + +def points_equal(left: dict[str, Any], right: dict[str, Any]) -> bool: + """Check whether ``right`` is fully covered by ``left`` and holds equal values. + + ``right`` is the *reference* (the start): every one of its entries must be + present in ``left`` and match. Callers use this to decide whether the + optimizer moved at all, so the reference has to be the *complete* start + point --- which is why ``find_MAP`` expands a user's partial + ``start=`` through ``make_start_point`` before auditing. Comparing only + the keys the two happen to share would otherwise call a run "unmoved" on the + evidence of one unchanged parameter, while the ones the user did not name + are exactly those that may have moved. + + Compared up to a tight relative tolerance rather than bit-for-bit: the + optimizer's point comes back through a transform and its inverse, so an + untouched parameter can differ from the start in the last ulp. This only + feeds a diagnostic message, so a near-match is the useful reading. + + The absolute tolerance carries that reading through zero. A relative + tolerance around ``0.0`` is a tolerance of exactly zero, so the round trip + turning an untouched ``0.0`` into ``1e-17`` would read as movement --- and + zero is the *usual* start for the group-level offsets of a hierarchical + model, which is precisely where the hint below is worth having. It is set + far below any step an optimizer that really moved would take. + """ + shared = set(right) + if not shared or not shared <= set(left): + return False + return all( + np.allclose( + np.asarray(left[key], dtype=float), + np.asarray(right[key], dtype=float), + rtol=1e-10, + atol=1e-12, + ) + for key in shared + ) + + +def audit_result( + opt_result: Any, + point: dict[str, Any], + start: dict[str, Any], +) -> list[str]: + """Return the reasons an optimizer run should be considered a failure. + + An empty list means the run is trustworthy. + + Parameters + ---------- + opt_result + A ``scipy.optimize.OptimizeResult``, or ``None`` --- which is the *only* + signal that the optimizer was interrupted or that ``maxeval`` was + exhausted. + point + The point the optimizer returned. + start + The *complete* point the optimizer started from. A partial point would + make the did-not-move hint below fire on incomplete evidence. + + Notes + ----- + A point identical to the start is not evidence of failure on its own: a + start that already sits at the optimum legitimately returns ``success=True`` + with ``nit=0``. It only sharpens the diagnosis of a run that has *already* + failed on another criterion, where it points at a non-finite gradient at the + start. + + That includes a run that ran out of ``maxeval``, which reports its *last + iterate* rather than its start --- so a point equal to the start there is + not the budget being too small but the very first evaluation having a + non-finite gradient, which is exactly what the hint says. In practice a + budget-limited run has moved and the hint stays silent. + """ + reasons: list[str] = [] + + if opt_result is None: + reasons.append( + "the optimizer was interrupted or exhausted its `maxeval` budget " + "(raise `maxeval=` if the run was still making progress)" + ) + else: + if not bool(getattr(opt_result, "success", True)): + message = str(getattr(opt_result, "message", "")).strip() + reasons.append( + "the optimizer reported failure" + + (f" ({message})" if message else "") + + f" after {getattr(opt_result, 'nit', '?')} iteration(s)" + ) + objective = getattr(opt_result, "fun", None) + if objective is not None and not np.all( + np.isfinite(np.asarray(objective, dtype=float)) + ): + reasons.append("the objective at the returned point is not finite") + + if reasons and points_equal(point, start): + reasons.append( + "the optimizer never moved off its starting point, which usually " + "means the gradient of the log-density is non-finite there. Pass a " + "different `start=` (HSSM's own `model.initvals` is the default) or " + "try `method='Powell'`" + ) + + return reasons + + +def describe_failures(failures: list[list[str]]) -> str: + """Render one message covering every start that failed. + + Reporting only the last start's reasons --- which is what a single + overwritten variable gives you --- hides the case where the starts failed + for *different* reasons, and that difference is the diagnosis. Identical + reasons are collapsed so the common case stays short. + """ + if not failures: + return "no start converged" + + unique: list[list[str]] = [] + for reasons in failures: + if reasons not in unique: + unique.append(reasons) + + if len(unique) == 1: + detail = "; ".join(unique[0]) + if len(failures) > 1: + return f"all {len(failures)} starts failed the same way: {detail}" + return detail + + return " | ".join( + f"start {index + 1}: " + "; ".join(reasons) + for index, reasons in enumerate(failures) + ) + + +def report_failure(kind: str, failures: list[list[str]], strict: bool) -> None: + """Warn --- or raise under ``strict`` --- about failed optimizer runs. + + Parameters + ---------- + kind + ``"MAP"`` or ``"MLE"``. + failures + One list of reasons per start that failed, in the order they were tried. + strict + Raise a ``RuntimeError`` instead of warning. + """ + # Deliberately does not point at `return_raw=True`: a failed run returns + # `None` before that branch is reached, so it is not an inspection channel. + # The two branches differ in what actually happened, so they cannot share a + # closing sentence: under `strict` nothing is returned (this raises) and + # recommending `strict=True` would be advising the mode already in force. + common = ( + f"find_{kind} failed to converge: " + + describe_failures(failures) + + ". No estimate was stored. " + "Try `n_starts>1`, a different `method=`, or a different `start=`" + ) + if strict: + raise RuntimeError(f"{common}.") + warnings.warn( + f"{common}; `strict=True` raises instead of warning. " + f"`find_{kind}` returned None.", + UserWarning, + stacklevel=3, + ) + + +def validate_n_starts(n_starts: int) -> None: + """Reject a ``n_starts`` the multi-start loop cannot honour. + + ``jittered_starts`` always yields the unjittered start first, so anything + below 1 would silently run exactly one start while the estimate reported + ``n_starts=0``. Refusing it keeps the recorded metadata truthful. + """ + if n_starts < 1: + raise ValueError(f"n_starts must be at least 1, got {n_starts}.") + + +def jittered_starts( + initvals: dict[str, Any], + n_starts: int, + rng: np.random.Generator, + jitter: float = 0.1, +): + """Yield ``n_starts`` starting points: the given one, then jittered copies. + + The jitter is *relative* (``value * (1 + U(-jitter, jitter))``), with an + additive fallback for entries that are exactly zero, where a relative + perturbation would be no perturbation at all. + + Both callers pass a point in *transformed* (unconstrained) space, and that + is what keeps the jitter inside the model's support. A relative jitter is + safe at a lower bound of zero but not at an upper one --- on the constrained + scale it would take ``z = 0.95`` to ``1.045``, outside the interval + ``z`` lives in, whose forward transform is then ``NaN`` and whose start PyMC + rejects. In transformed space there are no bounds left to cross: the + transforms map the whole real line onto the supported range, so every + jittered candidate comes back in support. + + Unlike ``HSSMBase._jitter_initvals`` this never mutates the model's stored + initial values and draws from an explicit ``Generator``, so multi-start runs + are reproducible and leave the model untouched. + """ + yield dict(initvals) + for _ in range(max(n_starts - 1, 0)): + candidate = {} + for name, value in initvals.items(): + array = np.asarray(value) + # The perturbation is multiplicative, so it has to land in floating + # point. Casting back to an integer dtype would truncate + # ``2 * (1 + 0.02)`` to ``2`` --- no jitter at all, while the + # estimate still reports ``n_starts`` starts --- or down to ``1``, a + # 50% jump rather than the requested 10%. + dtype = ( + array.dtype if np.issubdtype(array.dtype, np.floating) else np.float64 + ) + noise = rng.uniform(-jitter, jitter, array.shape) + perturbed = np.where(array == 0, noise, array * (1.0 + noise)) + candidate[name] = perturbed.astype(dtype).reshape(array.shape) + yield candidate + + +def value_var_names(pymc_model: "PyMCModel") -> list[str]: + """Return the names of the model's value variables.""" + return [var.name for var in pymc_model.value_vars] + + +def free_rv_names(pymc_model: "PyMCModel") -> list[str]: + """Return the names of the model's free random variables (constrained).""" + return [rv.name for rv in pymc_model.free_RVs] + + +def constrained_params( + pymc_model: "PyMCModel", point: dict[str, Any] +) -> dict[str, Any]: + """Extract the constrained free parameters from a raw optimizer point. + + The raw point mixes constrained values, transformed value-variable names and + trial-wise deterministics; only the first group is a parameter estimate. + """ + return {name: point[name] for name in free_rv_names(pymc_model) if name in point} + + +def drop_transformed(pymc_model: "PyMCModel", point: dict[str, Any]) -> dict[str, Any]: + """Drop the transformed value-variable entries from a raw optimizer point. + + Mirrors ``pm.find_MAP(include_transformed=False)``, which HSSM cannot pass + through directly because scoring and the failure audit need the transformed + names. + """ + # A value variable whose name differs from its RV's is the transformed one + # (``t`` -> ``t_log__``); the rest are already on the constrained scale. + transformed = set(value_var_names(pymc_model)) - set(free_rv_names(pymc_model)) + return {name: value for name, value in point.items() if name not in transformed} + + +def dims_and_coords( + pymc_model: "PyMCModel", names: list[str] +) -> tuple[dict[str, list[str]], dict[str, list[Any]]]: + """Collect the ArviZ ``dims``/``coords`` metadata for ``names``. + + Only the coords the collected ``dims`` actually reference are kept. Copying + ``pymc_model.coords`` wholesale would hang the model's observation index off + every estimate --- an entry per trial, so hundreds or thousands of strings + that nothing ever reads, carried through every ``pickle`` of the result. + """ + dims = { + name: list(pymc_model.named_vars_to_dims[name]) + for name in names + if name in pymc_model.named_vars_to_dims + and pymc_model.named_vars_to_dims[name] is not None + } + used = {dim for entry in dims.values() for dim in entry} + coords = { + key: list(value) + for key, value in pymc_model.coords.items() + if value is not None and key in used + } + return dims, coords + + +def make_start_point( + pymc_model: "PyMCModel", + overrides: dict[str, Any] | None, + seed: int | None = None, + validate: bool = True, +) -> dict[str, Any]: + """Build a transformed-space start point from constrained overrides. + + ``HSSM.initvals`` is keyed by constrained RV names; the optimizer works on + transformed value variables. PyMC's override machinery applies the forward + transform for us, so no manual conversion is needed --- but only + ``pm.find_MAP`` calls it internally, which is why ``find_MLE`` has to do it + here. + """ + point_fn = make_initial_point_fn( + model=pymc_model, + jitter_rvs=set(), + return_transformed=True, + overrides=cast("Any", overrides), + ) + start = point_fn(seed) + if validate: + pymc_model.check_start_vals(start) + return start + + +def flat_gradient(objective: Any, wrt: list) -> Any: + """Return the gradient of ``objective`` w.r.t. ``wrt`` as one flat vector. + + Concatenating in ``wrt`` order matches the layout + ``DictToArrayBijection.map`` produces for the same variables, so the vector + lines up with the optimizer's parameter vector entry for entry. + + ``disconnected_inputs="ignore"`` is deliberate. PyTensor's default is to + raise when a variable does not reach the objective, but that is a legitimate + state here rather than a user error: the observed log-likelihood of a + centrally-parameterized hierarchical model genuinely does not contain the + group-level ``sigma``. Zero --- what ``return_disconnected="zero"`` supplies + --- *is* the gradient in that case, and the flat-ridge problem it signals is + reported up front by ``find_MLE``'s hierarchical check. + + ``rewrite_pregrad`` is applied first, as ``pymc.Model.dlogp`` does before its + own ``grad`` call. Those are PyTensor's ``canonicalize`` and ``stabilize`` + passes: on the exp/log chains that fill a log-density graph they are what + keeps the differentiated form from overflowing where the objective itself is + finite. Skipping them would leave every gradient built here --- the MLE + optimizer's, the availability probe's and the standard-error Hessian's --- + on shakier numerics than the gradient ``pymc.find_MAP`` uses for the very + same model. + """ + gradients = cast( + "list[Any]", + pytensor.grad(rewrite_pregrad(objective), wrt, disconnected_inputs="ignore"), + ) + return pt.concatenate([pt.atleast_1d(gradient).ravel() for gradient in gradients]) + + +def compile_point_fn(pymc_model: "PyMCModel", outputs: Any, inputs: list) -> Any: + """Compile ``outputs`` as a function of ``inputs``, tolerating unused ones. + + Every caller here compiles against the model's *full* set of value + variables, so that one raveled parameter vector can be mapped back onto the + inputs without bookkeeping. PyTensor's default is to reject an input the + graph never reads, which turns a legitimate objective --- the observed + log-likelihood, which excludes prior-only parameters --- into an opaque + ``UnusedInputError``. Ignoring unused inputs is what makes the mapping + uniform. + """ + return pymc_model.compile_fn(outputs, inputs=inputs, on_unused_input="ignore") + + +def make_point_expander( + pymc_model: "PyMCModel", +) -> Callable[[dict[str, Any]], dict[str, Any]]: + """Compile the map from a transformed-space point to a full point. + + The result carries the constrained values, the transformed value-variable + names and the deterministics --- the same shape + ``pm.find_MAP(include_transformed=True)`` returns. + + Compiled once and reused. The graph does not depend on the point, so a + multi-start run that rebuilt it inside the loop would pay the compile on + every candidate to expand points that all but one of them discards. + + Returns + ------- + Callable + A function taking a point over the model's value variables and + returning the expanded point. + """ + unobserved = get_default_varnames(pymc_model.unobserved_value_vars, True) + function = compile_point_fn(pymc_model, unobserved, pymc_model.value_vars) + names = [var.name for var in unobserved] + + def expand(point: dict[str, Any]) -> dict[str, Any]: + return dict(zip(names, function(point))) + + return expand + + +def compile_objective( + pymc_model: "PyMCModel", + objective: Any, +) -> tuple[Callable[[dict[str, Any]], Any], Callable[[dict[str, Any]], Any] | None]: + """Compile ``objective`` and, when available, its gradient. + + Both are returned as functions of a *point dictionary* over the model's + value variables. The gradient is returned as a single flat vector whose + entry order matches ``DictToArrayBijection.map`` over the continuous value + variables. + + Deliberately independent of any start point: the compiled graphs are the + same for every candidate in a multi-start run, so ``find_MLE`` compiles once + and reuses the result. Binding a start (which ``optimize_objective`` does + with ``DictToArrayBijection.mapf``) is the cheap step and stays per-run --- + it also supplies the values of any value variable that is *not* raveled into + the parameter vector, which must come from that run's own start. + """ + inputs = pymc_model.value_vars + logp_fn = compile_point_fn(pymc_model, objective, inputs) + + try: + gradient = flat_gradient(objective, pymc_model.continuous_value_vars) + dlogp_fn = compile_point_fn(pymc_model, gradient, inputs) + except NO_GRADIENT_ERRORS: + dlogp_fn = None + + return logp_fn, dlogp_fn + + +def optimize_objective( + pymc_model: "PyMCModel", + objective: Any, + start: dict[str, Any], + method: str = DEFAULT_METHOD, + maxeval: int = 5000, + progressbar: bool = True, + label: str = "MLE", + compiled: tuple[Callable[[dict[str, Any]], Any], Any] | None = None, + expand: Callable[[dict[str, Any]], dict[str, Any]] | None = None, + **minimize_kwargs, +) -> tuple[dict[str, Any], Any, str]: + """Maximize ``objective`` over the model's continuous value variables. + + A small stand-in for ``pymc.find_MAP`` for objectives PyMC cannot build + itself --- specifically the observed log-likelihood that ``find_MLE`` + maximizes. + + Parameters + ---------- + pymc_model + The model whose value variables are optimized. + objective + A scalar PyTensor variable to *maximize*. + start + A transformed-space start point covering all value variables. + method + A ``scipy.optimize.minimize`` method. Never ``None``: SciPy would then + auto-select a gradient-based method while ``jac`` is still set. Swapped + for ``GRADIENT_FREE_METHOD`` if the objective has no gradient --- + check the third return value rather than assuming. + maxeval + Maximum number of objective evaluations. Exhausting it yields + ``opt_result is None``, mirroring PyMC's behaviour (including its + off-by-one: the budget is checked *before* the counter is bumped, so + ``maxeval + 1`` evaluations are allowed). + progressbar + Whether to display a progress bar, as ``pymc.find_MAP`` does. + label + Title for that progress bar. + compiled + The ``(logp_fn, dlogp_fn)`` pair from ``compile_objective``, to reuse + across the starts of a multi-start run. Compiled here when omitted. + expand + The point expander from ``make_point_expander``, likewise reused across + starts. Compiled here when omitted. + minimize_kwargs + Forwarded to ``scipy.optimize.minimize``. + + Returns + ------- + tuple[dict, OptimizeResult | None, str] + The point --- constrained values, transformed names and deterministics, + exactly like ``pm.find_MAP(include_transformed=True)``; the raw + optimizer result (``None`` if interrupted or out of evaluations); and + the method that actually ran. + + Notes + ----- + When the run is interrupted or runs out of evaluations, the point returned + is the last iterate at which the objective (and gradient, when used) was + finite --- not the start. Falling back to the start instead would make every + exhausted run look like an optimizer that never moved, which is a different + failure with a different remedy. + """ + continuous = pymc_model.continuous_value_vars + if not continuous: + raise ValueError("Model has no unobserved continuous variables.") + + names = [var.name for var in continuous] + x0 = DictToArrayBijection.map({name: start[name] for name in names}) + raw_logp_fn, raw_dlogp_fn = ( + compile_objective(pymc_model, objective) if compiled is None else compiled + ) + # Bound to *this* run's start, so the value variables that stay outside the + # raveled parameter vector take their values from it. + logp_fn = DictToArrayBijection.mapf(raw_logp_fn, start) + dlogp_fn = ( + None if raw_dlogp_fn is None else DictToArrayBijection.mapf(raw_dlogp_fn, start) + ) + + use_gradient = dlogp_fn is not None and method != GRADIENT_FREE_METHOD + if dlogp_fn is None and method != GRADIENT_FREE_METHOD: + _logger.warning( + "Gradient not available for this likelihood; falling back to the " + "gradient-free '%s' method.", + GRADIENT_FREE_METHOD, + ) + method = GRADIENT_FREE_METHOD + + evaluations = 0 + previous_x = x0.data + progress = CustomProgress( + *Progress.get_default_columns(), + TextColumn("{task.fields[loss]}"), + console=Console(theme=default_progress_theme), + disable=not progressbar, + ) + task = progress.add_task(label, total=maxeval, loss="") + + def cost(x: np.ndarray): + nonlocal evaluations, previous_x + raveled = RaveledVars(x, x0.point_map_info) + value = -np.float64(logp_fn(raveled)) + + gradient = None + if use_gradient: + gradient = -np.asarray(dlogp_fn(raveled), dtype=np.float64) # type: ignore[misc] + # Only a point where the gradient is usable is worth reporting back: + # it is the last place the optimizer could actually have moved from. + # Copied, not referenced: SciPy hands over an array it owns, and + # whether it reallocates or writes through it on the next iteration + # is an implementation detail that varies by method and version. + if np.all(np.isfinite(gradient)): + previous_x = np.array(x, copy=True) + elif np.isfinite(value): + previous_x = np.array(x, copy=True) + + if evaluations > maxeval: + raise StopIteration( + f"Maximum number of objective evaluations ({maxeval}) reached." + ) + evaluations += 1 + # Refreshing the rendered loss on every evaluation costs more than the + # objective does on cheap likelihoods, so throttle it the way PyMC's + # own optimizer progress bar throttles. + if evaluations % 10 == 0: + progress.update(task, completed=evaluations, loss=f"logp = {-value:,.5g}") + else: + progress.update(task, completed=evaluations) + + return (value, gradient) if use_gradient else value + + with progress: + try: + opt_result = sp_optimize.minimize( + cost, x0.data, method=method, jac=use_gradient, **minimize_kwargs + ) + best_x = opt_result["x"] + except (KeyboardInterrupt, StopIteration) as error: + best_x, opt_result = previous_x, None + _logger.info(str(error)) + finally: + progress.update(task, completed=evaluations, refresh=True) + + raveled = RaveledVars(np.asarray(best_x), x0.point_map_info) + expand = make_point_expander(pymc_model) if expand is None else expand + point = expand(DictToArrayBijection.rmap(raveled, start)) + + return point, opt_result, method + + +def make_scorer( + pymc_model: "PyMCModel", + observed_only: bool = False, + function: Callable[[dict[str, Any]], Any] | None = None, +) -> Callable[[dict[str, Any]], float]: + """Compile a reusable log-density scorer for ``pymc_model``. + + Compiling once and reusing the result matters for multi-start runs, which + would otherwise pay the compile cost per candidate. + + Parameters + ---------- + pymc_model + The model to score against. + observed_only + If ``True``, score the observed log-likelihood only --- the objective + ``find_MLE`` maximizes. MAP and MLE estimates are only comparable when + both are scored on the *same* objective. + function + An already-compiled evaluator of the same objective over the model's + value variables, as ``compile_objective`` returns. ``find_MLE`` scores + exactly the objective it optimizes, so passing that one in saves + compiling the observed log-likelihood a second time per call. Compiled + here when omitted. + + Returns + ------- + Callable + A function taking a point (every value variable, transformed names + included) and returning its log-density, or ``-inf`` if it could not be + evaluated. It raises ``KeyError`` if the point is missing any value + variable --- unscoreable input, rather than a point of zero density. + """ + if function is None: + objective = ( + pymc_model.observedlogp + if observed_only + else pymc_model.logp(jacobian=False) + ) + function = compile_point_fn(pymc_model, objective, pymc_model.value_vars) + scoring_fn = function + names = value_var_names(pymc_model) + + def score(point: dict[str, Any]) -> float: + missing = [name for name in names if name not in point] + if missing: + # Not a `-inf`: that is the answer for a point of zero density, and + # reusing it for a malformed point would report "infinitely bad" + # where the truth is "unscoreable". The likeliest way to get here + # is scoring a `find_MAP(include_transformed=False)` estimate, + # which drops exactly the transformed entries needed below, so the + # message names that case rather than only the missing keys. + raise KeyError( + f"Cannot score this point: {', '.join(missing)} missing from " + "it. Scoring requires every value variable, transformed names " + "included; an estimate produced with " + "`include_transformed=False` does not carry them. Re-run " + "without that flag, or score the raw optimizer point." + ) + try: + return float(scoring_fn({name: point[name] for name in names})) + except (ValueError, FloatingPointError): # pragma: no cover - defensive + return -np.inf + + return score + + +def score_point( + pymc_model: "PyMCModel", + point: dict[str, Any], + observed_only: bool = False, +) -> float: + """Evaluate the model's log-density at ``point``. + + A one-shot convenience wrapper around ``make_scorer``; prefer the latter + when scoring more than one point against the same model. + + Parameters + ---------- + pymc_model + The model to score against. + point + A point containing every value variable (transformed names included). + observed_only + If ``True``, score the observed log-likelihood only --- the objective + ``find_MLE`` maximizes. MAP and MLE estimates are only comparable when + both are scored on the *same* objective. + + Returns + ------- + float + The log-density, or ``-inf`` if it could not be evaluated. + + Raises + ------ + KeyError + If ``point`` is missing any of the model's value variables, as an + estimate produced with ``include_transformed=False`` is. + """ + return make_scorer(pymc_model, observed_only=observed_only)(point) + + +def objective_value(opt_result: Any) -> float | None: + """Return the objective SciPy reports at its solution, or ``None``. + + ``pm.find_MAP`` minimizes the *negative* log-density, so ``opt_result.fun`` + is that density negated at the returned point --- the same number a freshly + compiled scorer produces, bit-for-bit, under L-BFGS-B, Powell, Nelder-Mead + and BFGS. Reading it costs nothing, where compiling a scorer costs a full + logp graph. + + ``None`` when the value is unusable: a SciPy method that leaves ``fun`` + unset, a vector-valued residual rather than a scalar, or a non-finite + entry. Callers fall back to scoring the point directly. + + Parameters + ---------- + opt_result + The ``OptimizeResult`` returned alongside the point. + + Returns + ------- + float | None + The reported objective, or ``None`` if it cannot be trusted. + """ + value = getattr(opt_result, "fun", None) + if value is None: + return None + array = np.asarray(value, dtype=float).ravel() + if array.size != 1 or not np.isfinite(array[0]): + return None + return float(array[0]) + + +def standard_errors( + pymc_model: "PyMCModel", + point: dict[str, Any], + observed_only: bool = False, +) -> dict[str, Any] | None: + """Compute standard errors from the inverse Hessian at ``point``. + + The Hessian is taken on the *untransformed* model, so the result is on the + scale users report. This matters twice over: ``pm.find_hessian`` works in + transformed value-variable space with ``jacobian=True``, which is both the + wrong scale *and* the curvature of a different function from the one + ``find_MAP`` maximizes (``jacobian=False``). + + The Hessian itself is obtained by central differences of the analytic + gradient. Symbolic second derivatives do not exist for ONNX/JAX-wrapped + (``approx_differentiable``) likelihoods, and this way the analytical and the + LAN path share one implementation. + + Parameters + ---------- + pymc_model + The model to differentiate. + point + The estimate, containing every constrained parameter. + observed_only + Take the curvature of the observed log-likelihood rather than the joint + log-density --- the observed information, which is what an MLE wants. + + Returns + ------- + dict | None + Errors keyed by constrained parameter name, or ``None`` if no gradient + is available (blackbox likelihoods). All entries are ``NaN`` when the + Hessian is unusable --- singular, not negative definite, or built from a + stencil that leaves the model's support, as happens at an optimum + sitting on a bound. It is all of them rather than the offending one + because the errors come from inverting the joint Hessian. + + For ``observed_only=True`` these are standard errors in the usual sense. + For the joint log-density they are posterior standard deviations under a + Laplace approximation around the mode; the arithmetic is identical but + the object is not. + + Warnings + -------- + The cost grows with the number of scalar parameters ``n``: the central + differences need ``2 * n`` full-data gradient evaluations plus ``2 * n`` + cheaper log-density evaluations for the support check, followed by an + ``n x n`` inverse. On a hierarchical model with hundreds of group-level + offsets, ``se=True`` can take considerably longer than the optimization it + follows. A line is logged at ``INFO`` before starting when ``n`` is large. + """ + untransformed = remove_value_transforms(pymc_model) + objective = ( + untransformed.observedlogp + if observed_only + else untransformed.logp(jacobian=False) + ) + continuous = untransformed.continuous_value_vars + names = [var.name for var in continuous] + + missing = [name for name in names if name not in point] + if missing: # pragma: no cover - defensive + _logger.warning( + "Cannot compute standard errors: %s missing from the estimate.", + ", ".join(missing), + ) + return None + + base = { + name: np.asarray(point[name], dtype=untransformed.named_vars[name].dtype) + for name in names + } + x0 = DictToArrayBijection.map(base) + + try: + gradient = flat_gradient(objective, continuous) + gradient_fn = DictToArrayBijection.mapf( + compile_point_fn(untransformed, gradient, untransformed.value_vars), + base, + ) + except NO_GRADIENT_ERRORS: + warnings.warn( + "Standard errors are not available for this likelihood: it exposes " + "no gradient, so the Hessian cannot be approximated. Returning the " + "estimate without standard errors.", + UserWarning, + # standard_errors <- _make_point_estimate <- find_MAP/find_MLE <- user + stacklevel=4, + ) + return None + + def gradient_at(x: np.ndarray) -> np.ndarray: + return np.asarray(gradient_fn(RaveledVars(x, x0.point_map_info)), dtype=float) + + # The gradient alone cannot tell us whether a stencil point is inside the + # model's support: a bound enters the graph through `pt.switch`, so outside + # it the log-density is `-inf` while its gradient stays perfectly finite. + # Only the density itself reveals the boundary, hence this second function. + logp_fn = DictToArrayBijection.mapf( + compile_point_fn(untransformed, objective, untransformed.value_vars), base + ) + + def in_support(x: np.ndarray) -> bool: + return bool(np.isfinite(logp_fn(RaveledVars(x, x0.point_map_info)))) + + size = x0.data.size + if size > _SE_SIZE_HINT: + _logger.info( + "Computing standard errors for %d parameters: %d full gradient " + "evaluations plus a %dx%d matrix inverse. This may take a while.", + size, + 2 * size, + size, + size, + ) + # Scaled to the precision the *gradient* is computed at, not to float64's. + # The central-difference step that balances truncation against rounding goes + # as the cube root of the relative noise in the differentiated quantity, and + # under `set_floatX("float32")` that noise is ~1e-7 rather than ~2e-16 --- so + # a float64-sized step is roughly three orders of magnitude too small there + # and divides the gradient's rounding noise by far too little. On a smooth + # analytical likelihood the resulting Hessian is still usable (the errors + # move by a fraction of a percent); the margin matters on the noisier + # ONNX/JAX gradients, which is also where float32 is actually used. + step = np.cbrt(np.finfo(pytensor.config.floatX).eps) * np.maximum( + 1.0, np.abs(x0.data) + ) + hessian = np.empty((size, size)) + # An estimate pinned to a bound puts one of its two stencil points outside + # the support. The gradient there is finite, so the resulting Hessian is + # well conditioned and every guard below would pass --- reporting the + # curvature of the *smooth continuation* of a density that is `-inf` on one + # side as if it were a standard error. That is exactly the case the + # docstring promises `NaN` for, so the stencil has to be checked directly. + stencil_in_support = True + for index in range(size): + forward = x0.data.astype(float).copy() + backward = x0.data.astype(float).copy() + forward[index] += step[index] + backward[index] -= step[index] + if not (in_support(forward) and in_support(backward)): + stencil_in_support = False + break + hessian[:, index] = (gradient_at(forward) - gradient_at(backward)) / ( + 2 * step[index] + ) + hessian = 0.5 * (hessian + hessian.T) + + errors = np.full(size, np.nan) + # Every entry, not just the offending one: the errors come from inverting + # the *joint* Hessian, so one unusable column corrupts the whole inverse. + if stencil_in_support and np.all(np.isfinite(hessian)): + try: + # `inv` alone is not a test of the documented precondition: it raises + # only on *exact* singularity, so a merely ill-conditioned Hessian + # comes back as large finite numbers that would pass the + # `variances > 0` check below and be reported as standard errors. + # Cholesky is the honest gate --- it fails on anything that is not + # positive definite --- and the condition number catches the + # near-singular remainder, where the inverse is numerical noise. + np.linalg.cholesky(-hessian) + if np.linalg.cond(-hessian) >= 1.0 / np.finfo(float).eps: + raise np.linalg.LinAlgError("Hessian is numerically singular.") + covariance = np.linalg.inv(-hessian) + variances = np.diag(covariance) + with np.errstate(invalid="ignore"): + errors = np.where(variances > 0, np.sqrt(variances), np.nan) + except np.linalg.LinAlgError: + pass + + if not np.all(np.isfinite(errors)): + warnings.warn( + "The Hessian at the estimate is singular or not negative definite, " + "so some standard errors are NaN. This is expected when the optimum " + "sits on a parameter bound or a parameter is weakly identified.", + UserWarning, + stacklevel=4, + ) + + return DictToArrayBijection.rmap(RaveledVars(errors, x0.point_map_info)) + + +def group_specific_terms(hssm_model) -> list[str]: + """Return the names of any group-specific (hierarchical) terms in the model. + + Detected structurally rather than by name-matching ``*_sigma``. Note that + ``bambi.Model.components`` is a ``dict``, so it must be iterated over + ``.values()`` --- bare iteration yields the component *names* and the check + would silently never fire. + """ + terms: list[str] = [] + for component in hssm_model.model.components.values(): + component_terms = getattr(component, "group_specific_terms", None) + if component_terms: + terms.extend(list(component_terms)) + return terms + + +# endregion diff --git a/tests/rl/test_rlssm.py b/tests/rl/test_rlssm.py index ad04eebc3..04688c939 100644 --- a/tests/rl/test_rlssm.py +++ b/tests/rl/test_rlssm.py @@ -432,6 +432,23 @@ def test_rlssm_sample_smoke(self, rldm_data, rlssm_config) -> None: ) assert trace is not None + @pytest.mark.slow + def test_rlssm_find_map_smoke(self, rldm_data, rlssm_config) -> None: + """MAP estimation should work on an RLSSM, extra fields included. + + RL models carry the trial-wise feedback the learning process consumes as + extra fields, and ``find_MAP`` runs the same extra-field refresh + ``sample()`` does before optimizing. This covers that path end to end; + it does not assert the refresh changes anything, because it currently + cannot -- see ``HSSMBase._point_estimate_setup``. + """ + model = RLSSM(data=rldm_data, model_config=rlssm_config) + estimate = model.find_MAP(progressbar=False) + + assert estimate is not None + assert estimate.success + assert "rl_alpha" in estimate.params + class TestRLSSMSimplifiedInterface: """Public model= kwarg API coverage.""" diff --git a/tests/test_initvals.py b/tests/test_initvals.py index f729d0bc7..8c6604f84 100644 --- a/tests/test_initvals.py +++ b/tests/test_initvals.py @@ -47,11 +47,19 @@ def test_sample_map(caplog, loglik_kind, model, sampler, initvals): ) cav_data = hssm.load_data("cavanagh_theta") caplog.set_level(logging.INFO) + # `_jitter_initvals` perturbs every initial value through the *unseeded* + # global numpy RNG, so `model.initvals` -- and therefore the MAP started + # from it -- would differ on every run, leaving the assertions below to be + # evaluated at a random point each time (which is how this test + # intermittently failed at float32). Turning the jitter off removes the + # dependence on global RNG state entirely; seeding it would only pin the + # draw, and nothing here is about jitter. model_on = hssm.HSSM( data=cav_data, model=model, loglik_kind=loglik_kind, process_initvals=True, + initval_jitter=0.0, ) initial_point = model_on.initial_point(transformed=True) @@ -77,6 +85,20 @@ def test_sample_map(caplog, loglik_kind, model, sampler, initvals): progressbar=False, ) + # This test used to be pure smoke, so it could not tell a real MAP from + # the bug it was meant to guard against: an optimizer that never left + # its starting point and had that start stored as the estimate. Note the + # module runs at float32, where `find_MAP` takes the gradient-free path, + # so the tolerance below is deliberately loose -- the claim is that the + # point moved, not where it landed. + estimate = model_on.map + assert estimate.success, "MAP estimation reported failure" + assert any( + not np.allclose(estimate.params[name], model_on.initvals[name], atol=1e-4) + for name in estimate.params + if name in model_on.initvals + ), "MAP estimate never moved off the initial point" + def _check_initval_defaults_correctness(model) -> None: """Check if initial values from default dictionary are correctly applied.""" diff --git a/tests/test_optimize.py b/tests/test_optimize.py new file mode 100644 index 000000000..45cb3baa4 --- /dev/null +++ b/tests/test_optimize.py @@ -0,0 +1,1270 @@ +"""Tests for point estimation --- ``find_MAP``, ``find_MLE`` and ``PointEstimate``. + +See issue #1102. +""" + +import pickle +import warnings + +import arviz as az +import cloudpickle +import jax +import numpy as np +import pandas as pd +import pymc as pm +import pytensor +import pytest +from scipy.optimize import OptimizeResult +from ssms.basic_simulators.simulator import simulator + +import hssm +from hssm.optimize import ( + REPR_MAX_ROWS, + PointEstimate, + audit_result, + describe_failures, + group_specific_terms, + jittered_starts, + points_equal, + score_point, +) + +TRUE_DDM = {"v": 0.5, "a": 1.5, "z": 0.5, "t": 0.5} +DATA_SEED = 10 + + +@pytest.fixture(autouse=True) +def _pin_float64_precision(): + """Pin every test in this module to float64 and restore the caller's state. + + Point estimation is gradient based and therefore precision sensitive: the + recovery tolerances below only hold at double precision, and the float32 + test deliberately flips the setting. 29 test modules call + ``hssm.set_floatX("float32")`` at *import* time and pytest imports every + module during collection, so without this fixture the precision a test runs + under depends on collection order. + + The raw configs are saved and restored directly rather than going through + ``hssm.set_floatX`` / ``set_jax_precision``, which would pull in the aDDM + kernel and its dtype cache. + """ + previous_floatx = pytensor.config.floatX + previous_x64 = jax.config.jax_enable_x64 + pytensor.config.floatX = "float64" + jax.config.update("jax_enable_x64", True) + yield + pytensor.config.floatX = previous_floatx + jax.config.update("jax_enable_x64", previous_x64) + + +@pytest.fixture(scope="module") +def ddm_data(): + """Return 500 DDM trials simulated from ``TRUE_DDM``. + + Seeded on purpose. The recovery tests assert the estimate lands within + ``0.2`` of the truth, but at 500 trials the drift rate `v` is the least + well-identified parameter: across seeds its estimate sits around `0.585` + with an SD near `0.055`, so an unseeded draw puts the assertion only about + two SDs from its bound and the suite fails a small fraction of runs for no + reason. This seed keeps the whole module deterministic; its worst + deviation from truth is `0.058` across both MAP and MLE. + """ + simulated = simulator( + list(TRUE_DDM.values()), model="ddm", n_samples=500, random_state=DATA_SEED + ) + return pd.DataFrame( + np.column_stack([simulated["rts"][:, 0], simulated["choices"][:, 0]]), + columns=["rt", "response"], + ) + + +@pytest.fixture +def ddm_model(ddm_data): + """Return a plain analytical DDM model over ``ddm_data``.""" + return hssm.HSSM(data=ddm_data, model="ddm", loglik_kind="analytical") + + +@pytest.fixture(scope="module") +def hierarchical_data(): + """Return a three-participant slice of ``cavanagh_theta``. + + Small enough for the fast suite, and it still reproduces the non-finite + gradient at PyMC's default start that motivated this work. + """ + data = hssm.load_data("cavanagh_theta") + keep = sorted(data.participant_id.unique())[:3] + return data[data.participant_id.isin(keep)].reset_index(drop=True) + + +@pytest.fixture +def hierarchical_model(hierarchical_data): + """Return a DDM with a group-specific intercept on ``v``.""" + return hssm.HSSM( + data=hierarchical_data, + model="ddm", + loglik_kind="analytical", + include=[{"name": "v", "formula": "v ~ 1 + (1|participant_id)"}], + ) + + +def failed_optimize_result(): + """Return an ``OptimizeResult`` that looks like a failed L-BFGS-B run.""" + return OptimizeResult( + x=np.zeros(4), fun=np.nan, success=False, message="ABNORMAL:", nit=0 + ) + + +# region ===== find_MAP ===== +def test_find_map_recovers_parameters(ddm_model): + """MAP on a well-specified analytical DDM should recover the truth.""" + estimate = ddm_model.find_MAP(progressbar=False) + + assert estimate is not None + assert estimate.success + assert estimate.kind == "MAP" + for name, true_value in TRUE_DDM.items(): + assert estimate.params[name] == pytest.approx(true_value, abs=0.2), name + + +def test_find_map_defaults_start_to_hssm_initvals(ddm_model, monkeypatch): + """`find_MAP` must hand PyMC HSSM's processed initial values, not PyMC's.""" + seen = {} + original = pm.find_MAP + + def spy(*args, **kwargs): + seen["start"] = kwargs["start"] + return original(*args, **kwargs) + + monkeypatch.setattr(pm, "find_MAP", spy) + ddm_model.find_MAP(progressbar=False) + + assert seen["start"] == ddm_model.initvals + # The whole point of the fix: HSSM's start differs from PyMC's own. + assert seen["start"]["t"] != ddm_model.initial_point()["t"] + + +def test_find_map_converges_where_pymc_default_start_fails(hierarchical_model): + """Regression test for the silently-returned start point (issue #1102). + + At PyMC's default start (``t=2.0``) the gradient of this model's log-density + has non-finite entries. The *logp* is finite there, so ``check_start_vals`` + passes and L-BFGS-B aborts at ``nit=0`` --- which the old code stored as the + MAP without complaint. + """ + start = hierarchical_model.initial_point(transformed=True) + gradient = hierarchical_model.pymc_model.compile_dlogp(jacobian=False)(start) + assert not np.all(np.isfinite(gradient)), ( + "precondition: PyMC's default start must have a non-finite gradient" + ) + + estimate = hierarchical_model.find_MAP(progressbar=False) + + assert estimate is not None + assert estimate.success + assert estimate.opt_result.nit > 0 + # Assert the point actually moved, not merely that the optimizer iterated. + assert estimate.params["a"] != pytest.approx( + float(hierarchical_model.initvals["a"]) + ) + assert estimate.params["t"] != pytest.approx( + float(hierarchical_model.initvals["t"]) + ) + + +def test_find_map_reports_failure_instead_of_returning_the_start(hierarchical_model): + """A failed run warns, returns None, and is never cached.""" + bad_start = hierarchical_model.initial_point() + + with pytest.warns(UserWarning, match="find_MAP failed to converge"): + estimate = hierarchical_model.find_MAP( + start=bad_start, method="L-BFGS-B", progressbar=False + ) + + assert estimate is None + assert hierarchical_model._map_dict is None + with pytest.raises(ValueError, match="compute map first"): + _ = hierarchical_model.map + + +def test_find_map_failure_names_the_non_finite_gradient(hierarchical_model): + """The warning should point at the likely cause, not just say 'failed'.""" + with warnings.catch_warnings(record=True) as caught: + warnings.simplefilter("always") + hierarchical_model.find_MAP( + start=hierarchical_model.initial_point(), + method="L-BFGS-B", + progressbar=False, + ) + + messages = [str(warning.message) for warning in caught] + assert any("never moved off its starting point" in message for message in messages) + + +def test_find_map_strict_raises(hierarchical_model): + """`strict=True` turns the failure warning into an error.""" + with pytest.raises(RuntimeError, match="find_MAP failed to converge"): + hierarchical_model.find_MAP( + start=hierarchical_model.initial_point(), + method="L-BFGS-B", + progressbar=False, + strict=True, + ) + + +def test_find_map_failure_clears_an_earlier_estimate(ddm_model): + """A failed run must not leave the *previous* run's estimate reachable. + + The "a failed estimate is never cached" contract is what `model.map` and + `sample(initvals="map")`'s `_map_dict is None` guard both rest on, so it + has to hold on the second call as much as the first. + """ + ddm_model.find_MAP(progressbar=False) + assert ddm_model.map is not None + + with pytest.warns(UserWarning, match="find_MAP failed to converge"): + assert ddm_model.find_MAP(maxeval=3, progressbar=False) is None + + # The attribute, not just the property: it is what `sample` branches on. + assert ddm_model._map_dict is None + with pytest.raises(ValueError, match="Please compute map first"): + _ = ddm_model.map + + +def test_find_map_unprocessed_initvals_still_guarded(hierarchical_data): + """With `process_initvals=False` the fix no-ops, so the guard must catch it. + + Without initval processing (and without jitter) ``model.initvals`` falls + back to PyMC's own start --- the very point that breaks the optimizer. The + default start no longer helps here, so the failure guard is the only thing + standing between the user and a bogus estimate. + """ + model = hssm.HSSM( + data=hierarchical_data, + model="ddm", + loglik_kind="analytical", + include=[{"name": "v", "formula": "v ~ 1 + (1|participant_id)"}], + process_initvals=False, + initval_jitter=0.0, + ) + assert model.initvals["t"] == pytest.approx(float(model.initial_point()["t"])) + + with pytest.warns(UserWarning, match="find_MAP failed to converge"): + estimate = model.find_MAP(method="L-BFGS-B", progressbar=False) + + assert estimate is None + + +def test_find_map_multi_start(ddm_model): + """`n_starts` tries jittered starts and reports how many converged.""" + estimate = ddm_model.find_MAP(n_starts=3, seed=0, progressbar=False) + + assert estimate is not None + assert estimate.n_starts == 3 + assert estimate.n_converged >= 1 + assert estimate.logp >= ddm_model.find_MAP(progressbar=False).logp - 1e-6 + + +def test_find_map_multi_start_leaves_initvals_untouched(ddm_model): + """Multi-start jitter must not mutate the model's stored initial values.""" + before = {name: np.copy(value) for name, value in ddm_model.initvals.items()} + ddm_model.find_MAP(n_starts=3, seed=0, progressbar=False) + + for name, value in before.items(): + np.testing.assert_array_equal(ddm_model.initvals[name], value) + + +def test_find_map_return_raw(ddm_model): + """`return_raw=True` returns the `(point, OptimizeResult)` tuple pymc does.""" + point, raw = ddm_model.find_MAP(return_raw=True, progressbar=False) + + assert isinstance(point, PointEstimate) + assert isinstance(raw, OptimizeResult) + assert raw.success + + +def test_find_map_include_transformed_false(ddm_model): + """`include_transformed=False` drops the transformed value-var entries.""" + estimate = ddm_model.find_MAP(include_transformed=False, progressbar=False) + + assert "t" in estimate + assert "t_log__" not in estimate + assert "z_interval__" not in estimate + + +def _spy_on_extra_fields(model, monkeypatch, calls): + """Make `_check_extra_fields` report work to do and record both calls. + + The check must return something *truthy*, or the update it guards never + runs and the test cannot tell a working refresh from a missing one. + """ + + def check(self, *args): + calls.append("check") + return True + + def update(self, *args): + calls.append("update") + + monkeypatch.setattr(type(model), "_check_extra_fields", check) + monkeypatch.setattr(type(model), "_update_extra_fields", update) + + +def test_find_map_refreshes_extra_fields(ddm_model, monkeypatch): + """RL/aDDM extra fields must be refreshed, as `sample()` does.""" + calls = [] + _spy_on_extra_fields(ddm_model, monkeypatch, calls) + ddm_model.find_MAP(progressbar=False) + + assert calls[:2] == ["check", "update"] + + +def test_find_map_reports_the_method_that_actually_ran(ddm_data): + """A blackbox likelihood forces Powell --- the estimate must say so. + + `pm.find_MAP` swaps the method out internally and never reports it, so + recording the requested method would name an optimizer that never ran. + """ + model = hssm.HSSM(data=ddm_data, model="ddm", loglik_kind="blackbox") + estimate = model.find_MAP(progressbar=False) + + assert estimate is not None + assert estimate.method == "Powell" + + +def test_find_map_audits_against_the_expanded_start(ddm_model, monkeypatch): + """A partial `start=` must be expanded before the did-not-move audit. + + The hint is decided by comparing the returned point against the start, so a + one-key `start=` would let a single unchanged parameter stand in for "the + optimizer never moved" and send the user after the wrong cause. Asserted on + the argument rather than on the resulting message: with a partial start the + optimizer usually moves anyway, so a message-level assertion would pass + whether or not the expansion happened. + """ + seen = {} + original = hssm.optimize.audit_result + + def spy(opt_result, point, start): + seen["start"] = start + return original(opt_result, point, start) + + monkeypatch.setattr(hssm.optimize, "audit_result", spy) + ddm_model.find_MAP(start={"t": np.array(0.1)}, progressbar=False) + + # Every value variable, not just the one the caller named. + expected = {var.name for var in ddm_model.pymc_model.value_vars} + assert set(seen["start"]) == expected + assert len(expected) > 1 + + +def test_find_map_partial_start_is_still_honoured(ddm_model): + """Expanding the start for the audit must not discard the caller's values.""" + estimate = ddm_model.find_MAP(start={"t": np.array(0.1)}, progressbar=False) + + assert estimate is not None + assert estimate.success + + +def test_find_map_propagates_optimizer_usage_errors(ddm_model): + """A bad `method=` is the caller's mistake, not a rejected starting point. + + SciPy raises `ValueError` from inside `pm.find_MAP` for an unknown solver. + Catching it alongside the `SamplingError` that really does mean "PyMC + rejected this start" would relabel a typo as a convergence failure and send + the user off to debug their initial values instead. + """ + with pytest.raises(ValueError, match="Unknown solver"): + ddm_model.find_MAP(method="LBFGSB", progressbar=False) + + +def test_find_map_multi_start_jitter_cannot_leave_the_support(ddm_model): + """Jitter near an upper bound must not burn starts. + + `z` lives in an interval, and a *relative* jitter of a start at 0.95 reaches + 1.045 --- outside it, so the forward transform is NaN and PyMC rejects the + start. Applying the jitter in transformed space instead leaves no bound to + cross. Seeded: with this seed the constrained-space jitter put one of the + four extra starts out of range. + """ + near_bound = { + "v": np.array(0.5), + "a": np.array(1.5), + "z": np.array(0.95), + "t": np.array(0.3), + } + estimate = ddm_model.find_MAP( + start=near_bound, n_starts=5, seed=0, progressbar=False + ) + + assert estimate is not None + assert estimate.n_converged == 5 + + +def test_find_map_reports_every_failed_start(ddm_model): + """With `n_starts>1` the warning must cover all starts, not just the last. + + `maxiter=0` fails every start identically, which is the case the reporting + is expected to collapse rather than repeat three times. + """ + with warnings.catch_warnings(record=True) as caught: + warnings.simplefilter("always") + estimate = ddm_model.find_MAP( + n_starts=3, seed=0, progressbar=False, options={"maxiter": 0} + ) + + assert estimate is None + messages = [str(warning.message) for warning in caught] + assert any("all 3 starts failed the same way" in message for message in messages) + + +# endregion + + +# region ===== float32 guard ===== +def test_float32_warns_and_falls_back_to_powell(ddm_data): + """float32 must warn and, absent an explicit method, switch to Powell. + + At single precision L-BFGS-B stops early *and reports success*, so the + warning cannot be conditioned on the success flag. + """ + hssm.set_floatX("float32", update_jax=True) + model = hssm.HSSM(data=ddm_data, model="ddm", loglik_kind="analytical") + + with pytest.warns(UserWarning, match="float32"): + estimate = model.find_MAP(progressbar=False) + + assert estimate is not None + assert estimate.method == "Powell" + # The gradient-free path still lands near the truth, unlike L-BFGS-B here. + for name, true_value in TRUE_DDM.items(): + assert estimate.params[name] == pytest.approx(true_value, abs=0.25), name + + +def test_float32_honours_an_explicit_method(ddm_data): + """An explicit `method=` is respected, but the warning still fires.""" + hssm.set_floatX("float32", update_jax=True) + model = hssm.HSSM(data=ddm_data, model="ddm", loglik_kind="analytical") + + with pytest.warns(UserWarning, match="float32"): + estimate = model.find_MAP(method="L-BFGS-B", progressbar=False) + + assert estimate.method == "L-BFGS-B" + + +def test_no_float32_warning_at_float64(ddm_model): + """The precision warning must not fire at double precision.""" + with warnings.catch_warnings(record=True) as caught: + warnings.simplefilter("always") + ddm_model.find_MAP(progressbar=False) + + assert not [w for w in caught if "float32" in str(w.message)] + + +# endregion + + +# region ===== find_MLE ===== +def test_find_mle_recovers_parameters(ddm_model): + """MLE on a well-specified analytical DDM should recover the truth.""" + estimate = ddm_model.find_MLE(progressbar=False) + + assert estimate is not None + assert estimate.kind == "MLE" + for name, true_value in TRUE_DDM.items(): + assert estimate.params[name] == pytest.approx(true_value, abs=0.2), name + + +def test_mle_beats_map_on_the_same_objective(ddm_model): + """The MLE must not score below the MAP on the observed log-likelihood. + + Scoring both points on `observedlogp` is the whole comparison: the MAP + maximizes the joint density, so comparing its objective to the MLE's would + say nothing at all. + """ + map_estimate = ddm_model.find_MAP(progressbar=False) + mle_estimate = ddm_model.find_MLE(progressbar=False) + + map_score = score_point(ddm_model.pymc_model, map_estimate, observed_only=True) + mle_score = score_point(ddm_model.pymc_model, mle_estimate, observed_only=True) + + # A strict `>` flakes when the priors are near-flat or the optimum sits on a + # bound; the meaningful claim is that MLE does not do worse. + assert mle_score >= map_score - 1e-3 + + +def test_find_mle_raises_on_hierarchical_models(hierarchical_model): + """Group-level scale is unidentified without the priors --- refuse to fit.""" + assert group_specific_terms(hierarchical_model) == ["1|participant_id"] + + with pytest.raises(ValueError, match="not well posed for hierarchical models"): + hierarchical_model.find_MLE() + + +def test_find_mle_hierarchical_escape_hatch(hierarchical_model): + """`allow_unidentified=True` proceeds anyway, for penalized-likelihood work.""" + estimate = hierarchical_model.find_MLE(allow_unidentified=True, progressbar=False) + + assert estimate is not None + assert estimate.kind == "MLE" + # The whole point of the escape hatch is that the group-level terms are + # still estimated -- unidentifiably, but present. + assert "v_1|participant_id_sigma" in estimate.params + assert "v_1|participant_id_offset" in estimate.params + + +def test_find_mle_escape_hatch_on_a_centered_hierarchical_model(hierarchical_data): + """The escape hatch must survive a parameter the likelihood never reads. + + Under the centered parameterization `sigma` does not enter `observedlogp` at + all. Compiling the objective against every value variable then hits + PyTensor's unused-input check, and differentiating it hits the disconnected + -input check -- both of which have to be told that this is expected here. + """ + model = hssm.HSSM( + data=hierarchical_data, + model="ddm", + loglik_kind="analytical", + include=[{"name": "v", "formula": "v ~ 1 + (1|participant_id)"}], + noncentered=False, + ) + assert "v_1|participant_id_sigma" in [rv.name for rv in model.pymc_model.free_RVs] + + estimate = model.find_MLE(allow_unidentified=True, progressbar=False) + + assert estimate is not None + # `sigma` is absent from the objective, so its gradient is exactly zero and + # the optimizer leaves it where it started. That is the flat ridge the + # default refusal exists to warn about, not a crash. + assert estimate.params["v_1|participant_id_sigma"] == pytest.approx( + float(model.initvals["v_1|participant_id_sigma"]) + ) + + +def test_find_mle_reports_the_method_that_actually_ran(ddm_data): + """A blackbox likelihood forces Powell --- the estimate must say so.""" + model = hssm.HSSM(data=ddm_data, model="ddm", loglik_kind="blackbox") + estimate = model.find_MLE(progressbar=False) + + assert estimate is not None + assert estimate.method == "Powell" + + +def test_find_mle_maxeval_failure_does_not_blame_the_starting_point(ddm_model): + """Running out of budget is a different failure from a stalled gradient. + + The point reported after exhausting `maxeval` is the last iterate, not the + start, so the did-not-move hint cannot fire and mislead the user into + changing `start=` when the fix is a larger `maxeval`. + """ + with warnings.catch_warnings(record=True) as caught: + warnings.simplefilter("always") + estimate = ddm_model.find_MLE(maxeval=3, progressbar=False) + + assert estimate is None + messages = [str(warning.message) for warning in caught] + assert any("maxeval" in message for message in messages) + assert not any("never moved off" in message for message in messages) + + +def test_find_mle_accepts_progressbar(ddm_model): + """`progressbar=` works on both methods; `find_MAP`'s docs invite the idiom.""" + estimate = ddm_model.find_MLE(progressbar=False) + + assert estimate is not None + + +def test_find_mle_raises_on_potentials(ddm_model): + """`observedlogp` silently drops potentials, so refuse the model outright.""" + with ddm_model.pymc_model: + pm.Potential("extra_term", ddm_model.pymc_model.free_RVs[0] * 0.0) + + with pytest.raises(ValueError, match="pm.Potential"): + ddm_model.find_MLE(progressbar=False) + + +def test_mle_property_guard(ddm_model): + """`.mle` raises until `find_MLE` has produced an estimate.""" + with pytest.raises(ValueError, match="compute mle first"): + _ = ddm_model.mle + + ddm_model.find_MLE(progressbar=False) + assert ddm_model.mle is ddm_model._mle_dict + + +def test_find_mle_multi_start(ddm_model): + """`n_starts` tries jittered starts and keeps the best likelihood.""" + estimate = ddm_model.find_MLE(n_starts=3, seed=0, progressbar=False) + + assert estimate is not None + assert estimate.n_starts == 3 + assert estimate.n_converged >= 1 + assert estimate.logp >= ddm_model.find_MLE(progressbar=False).logp - 1e-6 + + +def test_find_mle_multi_start_compiles_the_objective_once(ddm_model, monkeypatch): + """The compile is hoisted out of the start loop --- it is the dominant cost. + + Asserted on the call count rather than on wall-clock, which is far too noisy + to distinguish one PyTensor compile from three. + """ + calls = [] + original = hssm.optimize.compile_objective + + def spy(*args, **kwargs): + calls.append(args[:2]) + return original(*args, **kwargs) + + monkeypatch.setattr(hssm.optimize, "compile_objective", spy) + ddm_model.find_MLE(n_starts=3, seed=0, progressbar=False) + + assert len(calls) == 1 + + +def test_find_mle_does_not_probe_the_gradient_separately(ddm_model, monkeypatch): + """`compile_objective` already answers the gradient question. + + Probing it beforehand would build the gradient graph of the observed + log-likelihood twice per call --- the expensive half of both operations. + """ + calls = [] + monkeypatch.setattr( + hssm.optimize, + "gradient_available", + lambda *args, **kwargs: calls.append(args) or True, + ) + ddm_model.find_MLE(progressbar=False) + + assert calls == [] + + +def test_find_mle_scores_with_the_compiled_objective(ddm_model, monkeypatch): + """The scorer must reuse the compiled objective, not compile a second copy. + + `find_MLE` scores candidates on exactly the objective it maximizes, so the + two compiles would be of the same graph. + """ + seen = {} + original = hssm.optimize.make_scorer + + def spy(pymc_model, observed_only=False, function=None): + seen["function"] = function + return original(pymc_model, observed_only=observed_only, function=function) + + monkeypatch.setattr(hssm.optimize, "make_scorer", spy) + ddm_model.find_MLE(progressbar=False) + + assert seen["function"] is not None + + +def test_find_mle_multi_start_expands_points_with_one_compile(ddm_model, monkeypatch): + """The point expander is hoisted out of the start loop like the objective. + + Every candidate's raveled result has to be mapped back to a full point, but + the graph doing it is the same for all of them --- and all but the winning + candidate's expansion is discarded. + """ + calls = [] + original = hssm.optimize.make_point_expander + + def spy(*args, **kwargs): + calls.append(args) + return original(*args, **kwargs) + + monkeypatch.setattr(hssm.optimize, "make_point_expander", spy) + ddm_model.find_MLE(n_starts=3, seed=0, progressbar=False) + + assert len(calls) == 1 + + +def test_find_mle_return_raw(ddm_model): + """`return_raw=True` returns the `(estimate, OptimizeResult)` tuple.""" + point, raw = ddm_model.find_MLE(return_raw=True, progressbar=False) + + assert isinstance(point, PointEstimate) + assert isinstance(raw, OptimizeResult) + + +def test_find_mle_strict_raises(ddm_model): + """`strict=True` turns the MLE failure warning into an error.""" + with pytest.raises(RuntimeError, match="find_MLE failed to converge"): + ddm_model.find_MLE(maxeval=3, progressbar=False, strict=True) + + +def test_find_mle_failure_clears_an_earlier_estimate(ddm_model): + """`find_MLE` owes `model.mle` the same contract `find_MAP` owes `model.map`.""" + ddm_model.find_MLE(progressbar=False) + assert ddm_model.mle is not None + + with pytest.warns(UserWarning, match="find_MLE failed to converge"): + assert ddm_model.find_MLE(maxeval=3, progressbar=False) is None + + assert ddm_model._mle_dict is None + with pytest.raises(ValueError, match="Please compute mle first"): + _ = ddm_model.mle + + +def test_find_mle_to_datatree_round_trips(ddm_model): + """An MLE feeds the ArviZ toolchain exactly as a MAP does.""" + posterior = ddm_model.find_MLE(progressbar=False).to_datatree().posterior + + assert posterior.sizes["chain"] == 1 + assert posterior.sizes["draw"] == 1 + assert set(TRUE_DDM) <= set(posterior.data_vars) + + +def test_find_mle_float32_warns_and_falls_back(ddm_data): + """The float32 safety net covers `find_MLE`, not just `find_MAP`.""" + hssm.set_floatX("float32", update_jax=True) + model = hssm.HSSM(data=ddm_data, model="ddm", loglik_kind="analytical") + + with pytest.warns(UserWarning, match="float32"): + estimate = model.find_MLE(progressbar=False) + + assert estimate is not None + assert estimate.method == "Powell" + + +@pytest.mark.parametrize("method_name", ["find_MAP", "find_MLE"]) +def test_n_starts_below_one_is_rejected(ddm_model, method_name): + """`n_starts=0` would run one start while reporting `n_starts=0`.""" + with pytest.raises(ValueError, match="n_starts must be at least 1"): + getattr(ddm_model, method_name)(n_starts=0, progressbar=False) + + +def test_find_mle_refreshes_extra_fields(ddm_model, monkeypatch): + """The extra-fields refresh applies to MLE too, not just MAP.""" + calls = [] + _spy_on_extra_fields(ddm_model, monkeypatch, calls) + ddm_model.find_MLE(progressbar=False) + + assert calls[:2] == ["check", "update"] + + +# endregion + + +# region ===== standard errors ===== +def test_standard_errors_are_reported_on_the_constrained_scale(ddm_model): + """`se=True` adds finite, plausibly-sized standard errors.""" + estimate = ddm_model.find_MAP(se=True, progressbar=False) + + assert estimate.se is not None + frame = estimate.to_dataframe() + assert "se" in frame.columns + assert np.all(np.isfinite(frame["se"].to_numpy())) + # With 500 trials the errors should be small but non-zero. + assert np.all(frame["se"].to_numpy() > 0) + assert np.all(frame["se"].to_numpy() < 0.5) + + +def test_standard_errors_for_mle_use_the_observed_information(ddm_model): + """MLE standard errors come from `observedlogp`, not the joint density.""" + estimate = ddm_model.find_MLE(se=True) + + assert estimate.se is not None + assert np.all(np.isfinite(list(estimate.se.values()))) + + +def test_standard_errors_are_nan_when_the_hessian_is_not_negative_definite( + ddm_model, monkeypatch +): + """A Hessian that is not negative definite must yield NaN, not a number. + + Driven through the Cholesky gate rather than by hunting for a model with a + degenerate optimum: the branch under test is "the curvature check failed", + and how it failed makes no difference to what is reported. + """ + + def not_negative_definite(matrix): + raise np.linalg.LinAlgError("Matrix is not positive definite") + + monkeypatch.setattr(np.linalg, "cholesky", not_negative_definite) + + with pytest.warns(UserWarning, match="singular or not negative definite"): + estimate = ddm_model.find_MAP(se=True, progressbar=False) + + assert estimate.se is not None + assert np.all(np.isnan(np.asarray(list(estimate.se.values()), dtype=float))) + + +def test_standard_errors_refuse_an_estimate_sitting_on_a_bound(ddm_model): + """A boundary estimate must give NaN, not the curvature of a continuation. + + The differences are taken in constrained space, so a parameter on its bound + puts one stencil point outside the support. The log-density is `-inf` + there, but the bound enters the graph through `pt.switch`, so the *gradient* + stays finite --- the Hessian then comes back well conditioned and every + downstream guard passes, reporting a confident standard error for an + estimate that is pinned. Only checking the density catches it. + """ + estimate = ddm_model.find_MAP(progressbar=False) + on_bound = { + name: np.asarray(value).copy() for name, value in estimate.params.items() + } + on_bound["z"] = np.array(0.0) # `z` is Uniform(0, 1); this is its lower bound + + with pytest.warns(UserWarning, match="singular or not negative definite"): + errors = hssm.optimize.standard_errors(ddm_model.pymc_model, on_bound) + + assert errors is not None + # Every entry, not only `z`: the errors invert the joint Hessian, so one + # unusable column makes the whole inverse meaningless. + assert np.all(np.isnan(np.asarray(list(errors.values()), dtype=float))) + + +def test_standard_errors_still_finite_at_an_interior_optimum(ddm_model): + """The support check must not fire on an ordinary interior estimate.""" + with warnings.catch_warnings(): + warnings.simplefilter("error", UserWarning) + estimate = ddm_model.find_MAP(se=True, progressbar=False) + + assert estimate.se is not None + assert np.all(np.isfinite(np.asarray(list(estimate.se.values()), dtype=float))) + + +def test_standard_errors_unavailable_without_a_gradient(ddm_data): + """Blackbox likelihoods have no gradient: warn and return the estimate.""" + model = hssm.HSSM(data=ddm_data, model="ddm", loglik_kind="blackbox") + + with pytest.warns(UserWarning, match="Standard errors are not available"): + estimate = model.find_MAP(se=True, progressbar=False) + + assert estimate is not None + assert estimate.se is None + + +# endregion + + +# region ===== PointEstimate ===== +def test_point_estimate_is_a_dict(ddm_model): + """Existing code paths type-check the estimate with `isinstance(..., dict)`.""" + estimate = ddm_model.find_MAP(progressbar=False) + + assert isinstance(estimate, dict) + assert estimate["a"] == estimate.params["a"] + assert set(estimate.params) == {"v", "a", "z", "t"} + # The raw mapping keeps the transformed names and deterministics. + assert "t_log__" in estimate + + +def test_point_estimate_carries_only_the_coords_its_dims_use(ddm_model): + """The model's observation index has no business on a point estimate. + + `coords` used to be copied from the model wholesale, hanging one entry per + trial off every estimate --- never read, but pickled along with it. + """ + estimate = ddm_model.find_MAP(progressbar=False) + + used = {dim for entry in estimate.dims.values() for dim in entry} + assert set(estimate.coords) <= used + assert "__obs__" not in estimate.coords + # Still enough metadata for the ArviZ conversion to succeed. + assert estimate.to_datatree().posterior.sizes["draw"] == 1 + + +def test_point_estimate_to_dataframe(ddm_model): + """`to_dataframe` gives one row per scalar entry, no `se` column by default.""" + frame = ddm_model.find_MAP(progressbar=False).to_dataframe() + + assert list(frame.columns) == ["estimate"] + assert set(frame.index) == {"v", "a", "z", "t"} + + +def test_point_estimate_to_datatree_feeds_arviz(ddm_model): + """The estimate must round-trip into the ArviZ/HSSM toolchain.""" + estimate = ddm_model.find_MAP(progressbar=False) + datatree = estimate.to_datatree() + + assert datatree.posterior.sizes["chain"] == 1 + assert datatree.posterior.sizes["draw"] == 1 + + summary = az.summary(datatree) + # arviz 1.3 column names --- the legacy `hdi_3%`/`hdi_97%` no longer exist. + assert {"mean", "sd", "eti89_lb", "eti89_ub", "ess_bulk", "r_hat"} <= set( + summary.columns + ) + # A single draw has no sampling distribution: sd collapses and the + # convergence diagnostics are undefined. Documented, not a bug. + assert np.allclose(summary["sd"].to_numpy(), 0.0) + assert np.all(np.isnan(summary["r_hat"].to_numpy().astype(float))) + + ddm_model.restore_traces(datatree) + assert ddm_model.traces is datatree + + +def test_point_estimate_to_datatree_keeps_vector_dims(hierarchical_model): + """Vector-valued (group-level) parameters keep their dims and coords.""" + estimate = hierarchical_model.find_MAP(progressbar=False) + posterior = estimate.to_datatree().posterior + + offsets = posterior["v_1|participant_id_offset"] + assert offsets.sizes["chain"] == 1 + assert offsets.sizes["draw"] == 1 + assert offsets.ndim == 3 + + +def test_point_estimate_survives_pickle(ddm_model): + """Extra attributes on a dict subclass are easy to lose in serialization.""" + estimate = ddm_model.find_MAP(se=True, progressbar=False) + restored = pickle.loads(pickle.dumps(estimate)) + + assert isinstance(restored, PointEstimate) + assert restored.kind == estimate.kind + assert restored.success == estimate.success + assert restored.method == estimate.method + assert restored.logp == pytest.approx(estimate.logp) + assert restored.se.keys() == estimate.se.keys() + assert restored.params.keys() == estimate.params.keys() + assert dict(restored).keys() == dict(estimate).keys() + + +def test_point_estimate_copy_keeps_the_metadata(ddm_model): + """`dict.copy` on a dict subclass returns a plain dict --- override it.""" + estimate = ddm_model.find_MAP(se=True, progressbar=False) + duplicate = estimate.copy() + + assert isinstance(duplicate, PointEstimate) + assert duplicate.kind == estimate.kind + assert duplicate.method == estimate.method + assert duplicate.params == estimate.params + assert duplicate.se.keys() == estimate.se.keys() + # A copy, not an alias: mutating one must not reach the other. + duplicate.params["a"] = 999.0 + assert estimate.params["a"] != 999.0 + + +def test_point_estimate_repr_truncates_long_estimates(hierarchical_model): + """A hierarchy has one row per offset; `repr` must not dump them all.""" + estimate = hierarchical_model.find_MAP(progressbar=False) + estimate.params = {f"p{index}": float(index) for index in range(REPR_MAX_ROWS + 5)} + + text = repr(estimate) + + assert "and 5 more row(s)" in text + assert f"p{REPR_MAX_ROWS - 1}" in text + assert f"p{REPR_MAX_ROWS}" not in text + # The full frame is still reachable, which is what the hint promises. + assert len(estimate.to_dataframe()) == REPR_MAX_ROWS + 5 + + +def test_empty_estimate_does_not_read_as_not_computed(ddm_model): + """`.map`'s guard is `is None`; an empty-but-computed estimate is computed.""" + ddm_model._map_dict = PointEstimate({}) + + assert ddm_model.map == {} + + +def test_point_estimate_is_exported_at_the_top_level(ddm_model): + """`isinstance` checks need the class without reaching into a submodule.""" + assert hssm.PointEstimate is PointEstimate + + +# endregion + + +# region ===== sample() integration ===== +def test_sample_initvals_map_raises_when_map_fails(ddm_model, monkeypatch): + """A failed MAP must not silently degrade to PyMC's default start. + + `find_MAP` caches nothing on failure, so falling through would hand bambi + `initvals=None` --- sampling from a point the user never asked for. + """ + monkeypatch.setattr(type(ddm_model), "find_MAP", lambda self, **kwargs: None) + + with pytest.raises(RuntimeError, match="did not converge"): + ddm_model.sample(initvals="map", draws=1, tune=1, chains=1, cores=1) + + +def test_sample_initvals_map_passes_constrained_params_only(ddm_model, monkeypatch): + """The handoff should carry the constrained free parameters and nothing else.""" + seen = {} + + def fake_fit(*args, **kwargs): + seen["initvals"] = kwargs["initvals"] + raise RuntimeError("stop here") + + ddm_model.find_MAP(progressbar=False) + monkeypatch.setattr(ddm_model.model, "fit", fake_fit) + + with pytest.raises(RuntimeError, match="stop here"): + ddm_model.sample(initvals="map", draws=1, tune=1, chains=1, cores=1) + + assert set(seen["initvals"]) == {"v", "a", "z", "t"} + + +def test_laplace_sampler_warns(ddm_model, monkeypatch, caplog): + """Bambi's laplace path cannot receive HSSM's initial values --- warn.""" + + def fake_fit(*args, **kwargs): + raise RuntimeError("stop here") + + monkeypatch.setattr(ddm_model.model, "fit", fake_fit) + + with caplog.at_level("WARNING", logger="hssm"): + with pytest.raises(RuntimeError, match="stop here"): + ddm_model.sample(sampler="laplace") + + assert "laplace" in caplog.text + assert "find_MAP" in caplog.text + + +# endregion + + +# region ===== helper units ===== +def test_audit_result_accepts_a_converged_run_that_did_not_move(): + """A start already at the optimum returns `nit=0` --- that is not a failure.""" + point = {"a": np.array(1.5)} + result = OptimizeResult(x=np.zeros(1), fun=-1.0, success=True, nit=0) + + assert audit_result(result, point, point) == [] + + +def test_audit_result_flags_a_none_result(): + """`opt_result is None` is the only signal for interrupt / maxeval.""" + reasons = audit_result(None, {"a": np.array(1.5)}, {"a": np.array(2.0)}) + + assert any("maxeval" in reason for reason in reasons) + + +def test_points_equal_requires_the_reference_to_be_covered(): + """A partial reference cannot establish that the optimizer did not move.""" + point = {"a": np.array(1.5), "t": np.array(0.3)} + + assert points_equal(point, {"a": np.array(1.5), "t": np.array(0.3)}) + # `t` moved, so the points differ even though `a` matches. + assert not points_equal(point, {"a": np.array(1.5), "t": np.array(0.9)}) + # A key the point does not carry means the comparison is not decidable. + assert not points_equal(point, {"a": np.array(1.5), "z": np.array(0.5)}) + assert not points_equal(point, {}) + + +def test_describe_failures_collapses_identical_reasons(): + """Five starts failing the same way should read as one sentence, not five.""" + text = describe_failures([["gradient is non-finite"]] * 5) + + assert text == "all 5 starts failed the same way: gradient is non-finite" + + +def test_describe_failures_keeps_distinct_reasons_apart(): + """When the starts fail differently, that difference *is* the diagnosis.""" + text = describe_failures([["rejected by PyMC"], ["ran out of maxeval"]]) + + assert "start 1: rejected by PyMC" in text + assert "start 2: ran out of maxeval" in text + + +def test_describe_failures_handles_no_failures(): + """`report_failure` is also reached when every start was skipped.""" + assert describe_failures([]) == "no start converged" + + +def test_audit_result_adds_the_gradient_hint_only_on_failure(): + """The did-not-move hint sharpens a failure; it never creates one.""" + point = {"a": np.array(1.5)} + + reasons = audit_result(failed_optimize_result(), point, point) + assert any("never moved off its starting point" in reason for reason in reasons) + + moved = audit_result(failed_optimize_result(), point, {"a": np.array(2.0)}) + assert not any("never moved off" in reason for reason in moved) + + +def test_jittered_starts_yields_the_original_first(): + """The unjittered start is always tried, so `n_starts=1` is a no-op.""" + initvals = {"t": np.array(0.025), "a": np.array(1.5)} + starts = list(jittered_starts(initvals, 3, np.random.default_rng(0))) + + assert len(starts) == 3 + assert starts[0] == initvals + # Relative jitter cannot push a positive parameter through zero. + assert all(candidate["t"] > 0 for candidate in starts) + assert starts[1] != starts[2] + + +def test_jittered_starts_does_not_truncate_integer_valued_entries(): + """An integer-dtype entry must still be jittered by the requested fraction. + + Casting the perturbed value back to the input dtype would round `2 * 1.02` + to `2` --- no jitter at all, while the estimate still claims `n_starts` + starts --- or down to `1`, a 50% jump rather than the requested 10%. + """ + starts = list( + jittered_starts({"a": np.array(2)}, 6, np.random.default_rng(1), jitter=0.1) + ) + jittered = [float(candidate["a"]) for candidate in starts[1:]] + + assert all(value != 2.0 for value in jittered) + assert all(1.8 <= value <= 2.2 for value in jittered) + + +def test_points_equal_sees_an_untouched_zero_as_unmoved(): + """A zero entry must survive the transform round trip as "did not move". + + Zero is the usual start for a hierarchical model's group-level offsets, and + the round trip through a transform and its inverse can return `1e-17`. A + purely relative tolerance is a tolerance of exactly zero there, which would + read that as movement and suppress the non-finite-gradient hint. + """ + assert points_equal({"offset": np.array(1e-17)}, {"offset": np.array(0.0)}) + assert not points_equal({"offset": np.array(0.5)}, {"offset": np.array(0.0)}) + + +def test_scoring_a_constrained_only_estimate_says_what_is_wrong(ddm_model): + """A point missing its transformed entries is unscoreable, not zero-density. + + `include_transformed=False` drops exactly the `*_log__`/`*_interval__` + entries the scorer indexes, so the failure has to name that cause. Silently + returning `-inf` would report the estimate as infinitely improbable. + """ + estimate = ddm_model.find_MAP(include_transformed=False, progressbar=False) + + with pytest.raises(KeyError, match="include_transformed=False"): + score_point(ddm_model.pymc_model, estimate) + + +def test_strict_failure_message_fits_the_strict_path(): + """Under `strict` the call raises, so the warning's wording does not fit. + + Nothing is returned, and recommending `strict=True` would be advising the + mode already in force. + """ + with pytest.raises(RuntimeError) as excinfo: + hssm.optimize.report_failure("MAP", [["it failed"]], strict=True) + + message = str(excinfo.value) + assert "No estimate was stored" in message + assert "returned None" not in message + assert "strict=True" not in message + + +def test_failure_message_does_not_advertise_return_raw(): + """A failed run returns None, so `return_raw=True` is not an escape hatch.""" + with pytest.warns(UserWarning, match="find_MAP failed to converge") as caught: + hssm.optimize.report_failure("MAP", [["it failed"]], strict=False) + + message = str(caught[0].message) + assert "return_raw" not in message + assert "n_starts" in message + + +# endregion + + +# region ===== slow ===== +@pytest.mark.slow +def test_find_map_full_hierarchical_model(): + """The full `cavanagh_theta` hierarchical fit --- the original repro. + + The fast suite exercises a three-participant slice; this covers the whole + 19-parameter model, which takes far too many L-BFGS-B iterations for a fast + test. + """ + model = hssm.HSSM( + data=hssm.load_data("cavanagh_theta"), + model="ddm", + loglik_kind="analytical", + include=[{"name": "v", "formula": "v ~ 1 + (1|participant_id)"}], + ) + estimate = model.find_MAP(progressbar=False) + + assert estimate is not None + assert estimate.success + assert estimate.params["a"] != pytest.approx(float(model.initvals["a"])) + # The group-level scale should not collapse: hierarchical MAP is well + # behaved once the start is right. + assert estimate.params["v_1|participant_id_sigma"] > 0.1 + + +@pytest.mark.slow +def test_point_estimate_survives_cloudpickle(ddm_model): + """`save_model` uses cloudpickle, which must preserve the extra attributes. + + Extra attributes on a `dict` subclass are exactly the kind of state that + quietly disappears through a serializer that treats the object as a plain + mapping. + """ + estimate = ddm_model.find_MAP(se=True, progressbar=False) + restored = cloudpickle.loads(cloudpickle.dumps(estimate)) + + assert isinstance(restored, PointEstimate) + assert restored.kind == "MAP" + assert restored.success + assert restored.logp == pytest.approx(estimate.logp) + assert restored.se.keys() == estimate.se.keys() + assert restored.dims == estimate.dims + + +@pytest.mark.slow +def test_save_load_does_not_carry_a_point_estimate(ddm_model, tmp_path): + """Model persistence rebuilds from constructor args and drops inference state. + + `HSSMBase.__getstate__` stores only the constructor arguments, so a saved + model carries no traces, no VI approximation and no point estimate --- those + are persisted (or not) separately. This pins that contract down for point + estimates so the behaviour is a decision rather than a surprise. + """ + ddm_model.find_MAP(progressbar=False) + ddm_model.save_model( + model_name="point_estimate_roundtrip", + base_path=tmp_path, + allow_absolute_base_path=True, + ) + + loaded = hssm.HSSM.load_model(path=tmp_path / "point_estimate_roundtrip") + + assert loaded._map_dict is None + with pytest.raises(ValueError, match="compute map first"): + _ = loaded.map + + +@pytest.mark.slow +@pytest.mark.parametrize("model_name", ["ddm", "angle"]) +def test_find_map_approx_differentiable(model_name): + """Point estimation works through the ONNX/JAX likelihoods too.""" + theta = { + "ddm": [0.5, 1.5, 0.5, 0.3], + "angle": [0.5, 1.5, 0.5, 0.3, 0.2], + }[model_name] + simulated = simulator(theta, model=model_name, n_samples=500) + data = pd.DataFrame( + np.column_stack([simulated["rts"][:, 0], simulated["choices"][:, 0]]), + columns=["rt", "response"], + ) + model = hssm.HSSM(data=data, model=model_name, loglik_kind="approx_differentiable") + + estimate = model.find_MAP(se=True, progressbar=False) + + assert estimate is not None + assert estimate.success + assert estimate.params["a"] == pytest.approx(1.5, abs=0.3) + assert estimate.params["v"] == pytest.approx(0.5, abs=0.3) + # Second derivatives do not exist symbolically for these Ops; standard + # errors come from finite differences of the analytic gradient instead. + assert estimate.se is not None + + +@pytest.mark.slow +def test_find_map_and_mle_on_a_regression_model(cavanagh_test): + """A fixed-effects regression is fine for MLE --- only group terms are not.""" + model = hssm.HSSM( + data=cavanagh_test, + model="ddm", + loglik_kind="analytical", + include=[{"name": "v", "formula": "v ~ 1 + theta"}], + ) + + map_estimate = model.find_MAP(progressbar=False) + mle_estimate = model.find_MLE(progressbar=False) + + assert map_estimate is not None + assert mle_estimate is not None + assert "v_theta" in map_estimate.params + assert ( + score_point(model.pymc_model, mle_estimate, observed_only=True) + >= score_point(model.pymc_model, map_estimate, observed_only=True) - 1e-3 + ) + + +# endregion