Skip to content

PyMCWrapper TrueMeasure - #640

Open
nomadbl wants to merge 11 commits into
developfrom
Composable_True_Measure
Open

nomadbl wants to merge 11 commits into
developfrom
Composable_True_Measure

Conversation

@nomadbl

@nomadbl nomadbl commented Oct 21, 2026

Copy link
Copy Markdown
Collaborator

Summary

  • Issue:
    Evaluate support for composable hierarchical TrueMeasure models
  • Algorithmic or API impact:
    This PR adds a TrueMeasure sub-class that allows sampling from distributions constructed by PyMC.
    This substantially extends the class of distributions which can be sampled by QMCPy "out of the box", and potentially extends application areas (which is the author's intent). The main limitation from the PyMC side is some missing implementations of icdf which can be (and is progressively) patched.
  • Commands run: See the demo notebook, demos/pymc_model.ipynb for usage examples, limitations and workarounds for them.
  • Dependencies: New optional dependencies pymc and pytensor needed for this feature.
  • Generated artifacts: demos/pymc_model.ipynb downloads data from https://github.com/stan-dev/posteriordb to run one of the examples. This data is not commited to the repository.

AI Assistance

  • No substantive AI assistance was used for this PR.
  • AI assistance substantively affected this PR, and I describe that use below.

Checklist

  • I linked the relevant issue or explained why none was needed.
  • I added or updated tests, docs, notebooks, or explained why they were not needed.
  • I verified any equations, citations, benchmarks, numerical claims, or external references changed in this PR.
  • I reviewed AI-assisted content for licensing, provenance, attribution, confidentiality, and security concerns.
  • I described any CI, dependency, notebook runtime, or generated artifact changes if applicable.

Implementation Overview

PyMC was chosen specifically as a widely used, fully featured, python based, framework. This was preferred over the initial proof-of-concept implementation.

PyMC is built on top of PyTensor. PyTensor constructs computational graphs which can be "compiled" (actually, linked) to a program. The distinguishing feature of PyTensor is the ability to transform graphs, so several different computations can be implicitly specified by a single graph that gets transformed in different ways.
PyMC graphs are composed of a combination of two kinds of nodes: random variables like normal_rv, and deterministic operations like pt.exp, pt.log etc.
PyMC Models are (a collection of) these graphs, stored in the model's basic_RVs, unobserved_RVs and a few others. PyMC contains utilities to transform the random variable nodes into the corresponding computations for the log probability, which is then used in the implementation of algorithms.
This PR builds on PyMC's utilities, replacing the random variable nodes with icdf nodes to transform the original graph into a graph performing inverse-transform-sampling. The code in the class itself compiles this transform by calling these utilities. Most of the code in the class itself is glue code transforming arrays to/from QMCPy expected format and the format ingested by the compiled transform function.

Related utilities to add - separate PR or here

Sometimes models get really complicated, with a large number of variables of different shapes etc.
PyMC deals with this by defaulting outputs to a dictionary format {"var_name": ndarray, ...}. QMCPy seems to require/assume flat vectors as the "currency", thus I defaulted PyMCWrapper to output such flat vectors and the user is expected to know how to reconstruct the variables.
A useful utility to accompany here would be a "converter" that reconstructs the dictionary format from the flat format, and the converter would be automatically deduced from the Model passed in. I assumed this was appropriate as a separate PR but this can be added here without much effort.

f"'{_missing_module}'. Install the pymc >= 6.0.1"
)

class PyMCWrapper(AbstractTrueMeasure):
Comment thread test/test_tm_true_measures.py Fixed
Comment thread test/test_tm_true_measures.py Fixed
Comment on lines +417 to +466
# def _detect_issues(icdf_f, dims, split_indices, n_grid=3):
# issues = []

# eps = np.finfo(float).eps
# us = np.meshgrid(*[np.linspace(eps, 1.0 - eps, n_grid) for _ in range(dims)], indexing='ij')
# us = [u.ravel() for u in us]
# us = np.stack(us, axis=-1)
# us = np.split(us, split_indices, axis=-1)
# us = [_x.ravel() for _x in us]
# if np.prod(us[-1].shape) == 0:
# us = us[:-1]
# try:
# x = icdf_f(*us)
# except Exception as e:
# return [f"icdf raised {type(e).__name__}: {e}"]

# if np.any(~np.isfinite(x)):
# issues.append("ppf() returned non-finite values")

# # TODO: can we generalize this for dims>1?
# # dx = np.diff(x)
# # if np.any(dx < -1e-12):
# # issues.append("ppf() is not nondecreasing (looks non-monotone)")

# # Light-touch pdf/logpdf checks on quantiles (works for bounded/unbounded support)
# # uq = np.linspace(0.01, 0.99, 9)
# # try:
# # xq = np.asarray(icdf_f(uq), dtype=float)
# # except Exception:
# # xq = None

# # if xq is not None and hasattr(dist, "pdf"):
# # try:
# # p = np.asarray(dist.pdf(xq), dtype=float)
# # if np.any(~np.isfinite(p)):
# # issues.append("pdf() returned non-finite values on quantiles")
# # if np.any(p < -1e-12):
# # issues.append("pdf() returned negative values on quantiles")
# # except Exception as e:
# # issues.append(f"pdf() raised {type(e).__name__}: {e}")

# # if xq is not None and hasattr(dist, "logpdf"):
# # try:
# # lp = np.asarray(dist.logpdf(xq), dtype=float)
# # if np.any(~np.isfinite(lp)):
# # issues.append("logpdf() returned non-finite values on quantiles")
# # except Exception as e:
# # issues.append(f"logpdf() raised {type(e).__name__}: {e}")

# return issues
Comment on lines +607 to +626
# try:
# # do some sampling and check no errors
# bs = 10
# inputs = np.ones((bs, self.d))*0.5
# # split to right array lengths and remove trailing empty array if present
# if self.d > 1:
# inputs_list = np.split(inputs, self.split_indices, axis=-1)
# else:
# inputs_list = [inputs.ravel()]
# inputs_list = [_x.ravel() for _x in inputs_list]
# if np.prod(inputs[-1].shape) == 0:
# inputs_list = inputs_list[:-1]
# # reshape to correct shapes
# print(self.shapes)
# print(inputs_list)
# inputs_list = [np.reshape(_x, bs, *shape) for _x, shape in zip(inputs_list, self.shapes)]
# _ = self.icdf(*inputs_list)
# except:
# raise DimensionError(
# f"sampler dimension {sampler.d} must match the model joint samples dimension {self.dimension}."
Comment on lines +630 to +635
# if issues:
# warnings.warn(
# "PyMCWrapper received a PyMC Model "
# "and it failed sanity checks:\n - " + "\n - ".join(issues),
# UserWarning,
# stacklevel=2,

import numpy as np
import warnings
from pytensor.tensor.utils import get_static_shape_from_size_variables
Comment thread test/test_tm_true_measures.py Fixed
@nomadbl nomadbl changed the title Composable true measure PyMCWrapper TrueMeasure Oct 21, 2026
self.models = []
# define models for tests
with pymc.Model() as normal_model:
x = pymc.Normal("x", mu=0.0, sigma=1.0)
nu = 5.0
with pymc.Model() as t_model:
s = pymc.ChiSquared("s", nu=nu)
x = pymc.Normal("x", mu=0.0, sigma=pymc.math.sqrt(nu/s))

with pymc.Model() as observed_t_model:
s = pymc.ChiSquared("s", nu=nu)
x = pymc.Normal("x", mu=0.0, sigma=pymc.math.sqrt(nu/s), observed=np.ones(3))
Comment thread test/test_tm_pymcwrapper.py Fixed
Comment on lines +8 to +11
from qmcpy import (
PyMCWrapper,
DigitalNetB2,
)
@sou-cheng-choi

sou-cheng-choi commented Oct 22, 2026 •

Copy link
Copy Markdown
Member

@nomadbl

Kindly review PR#642 to your branch. It is an attempt to fix the unit test failures. Hope it helps as discussed in Slack with @alegresor .

@sou-cheng-choi sou-cheng-choi left a comment •

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

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

See PR642.

Please requeste for re-review when the tests pass at 100%.

This branch has not been deployed

No deployments
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

3 participants