Integrators (MC Sampling)

Monte Carlo integration modules for likelihood evaluation and parameter inference.

Overview

Integration in RIFT is primarily focused on the computation of the evidence (marginal likelihood) and the characterization of the posterior distribution. Given the high dimensionality and potential complexity of the gravitational-wave parameter space, RIFT employs a variety of Monte Carlo (MC) integration strategies.

The core goal of these integrators is to evaluate integrals of the form:

\[Z = \int \mathcal{L}(\theta) \pi(\theta) d\theta\]

where \(\mathcal{L}\) is the likelihood and \(\pi\) is the prior.

Integration Strategies

RIFT provides several specialized integrators, each optimized for different regimes of the likelihood surface or computational budgets.

Standard MC Sampler (``mcsampler``) The default integrator. It performs basic Monte Carlo integration by sampling from the prior and averaging the likelihood. While robust, it can be inefficient for highly peaked likelihoods.

GPU-Accelerated Sampler (``mcsamplerGPU``) A high-performance implementation designed to leverage CUDA-enabled GPUs. By parallelizing the likelihood evaluations across thousands of GPU cores, it can achieve orders-of-magnitude speedups over the CPU implementation, especially when using vectorized likelihoods. Reference: See Wysocki et al. for details on GPU acceleration.

Ensemble / GMM Sampler (``mcsamplerEnsemble``) Based on Gaussian Mixture Model (GMM) approximations of the likelihood. This approach identifies the primary modes of the posterior and uses an ensemble of GMMs to efficiently sample the most important regions of the parameter space. Reference: See Wofford et al. for details on ensemble-based marginalization.

Adaptive Volume Sampler (``mcsamplerAdaptiveVolume``) An integrator that dynamically adjusts the sampling volume based on the observed likelihood distribution. By focusing samples in regions of high likelihood, it reduces the variance of the evidence estimate. Reference: See Wagner et al. for adaptive sampling techniques.

Normalizing Flows Sampler (``mcsamplerNFlow``) Leverages deep learning-based Normalizing Flows to learn a bijective mapping between a simple base distribution (e.g., a Gaussian) and the complex posterior. This allows for nearly exact sampling from the posterior and highly efficient integration. Reference: See Wagner et al..

Portfolio Integrator (``mcsamplerPortfolio``) A meta-integrator that manages a “portfolio” of different sampling strategies. It can dynamically allocate computational resources among various integrators to optimize the trade-off between accuracy and wall-clock time. Reference: See Wagner et al..

API Reference

The following sections provide the detailed API for the integration modules.

Base Integrators

class RIFT.integrators.mcsampler_generic.MCSamplerGeneric(**kwargs)

Bases: object

Class to define a set of parameter names, limits, and probability densities.

add_parameter(params, pdf, left_limit=None, right_limit=None, prior_pdf=None, adaptive_sampling=False, **kwargs)

Add one (or more) parameters to sample dimensions. params is either a string describing the parameter, or a tuple of strings. The tuple will indicate to the sampler that these parameters must be sampled together. left_limit and right_limit are on the infinite interval by default, but can and probably should be specified. If several params are given, left_limit, and right_limit must be a set of tuples with corresponding length. Sampling PDF is required, and if not provided, the cdf inverse function will be determined numerically from the sampling PDF.

draw_simplified(n_samples, *args, **kwargs)

draw n_samples

Should not be implemented in generic class

integrate(lnF, *args, xpy=<module 'numpy' from '/opt/hostedtoolcache/Python/3.11.18/x64/lib/python3.11/site-packages/numpy/__init__.py'>, **kwargs)

integrate_log(lnF, *args, xpy=<module 'numpy' from '/opt/hostedtoolcache/Python/3.11.18/x64/lib/python3.11/site-packages/numpy/__init__.py'>, **kwargs)

static match_params_from_args(args, params)

Given two unordered sets of parameters, one a set of all “basic” elements (strings) possible, and one a set of elements both “basic” strings and “combined” (basic strings in tuples), determine whether the sets are equivalent if no basic element is repeated. e.g. set A ?= set B (“a”, “b”, “c”) ?= (“a”, “b”, “c”) ==> True ((“a”, “b”, “c”)) ?= (“a”, “b”, “c”) ==> True ((“a”, “b”), “d”)) ?= (“a”, “b”, “c”) ==> False # basic element ‘d’ not in set B ((“a”, “b”), “d”)) ?= (“a”, “b”, “d”, “c”) ==> False # not all elements in set B represented in set A

setup(**kwargs)

update_sampling_prior(ln_weights, n_history, tempering_exp=1, log_scale_weights=True, floor_integrated_probability=0, external_rvs=None, **kwargs)

Default setting is for lnL. The choices adopted are fairly arbitrary. For simplicity, we are using ‘log scale weights’ (derived from lnL, not L) which lets us sample well away from the posterior at the expense of not being as strongly sampling near the peak.

NOTE: Currently deployed for mcsamplerPortfolio, NOT yet part of core code here

Standard MC Sampler

class RIFT.integrators.mcsampler.HealPixSampler(skymap, massp=1.0)

Bases: object

Class to sample the sky using a FITS healpix map. Equivalent to a joint 2-D pdf in RA and dec.

static decra2thph(dec, ra)

theta/phi to RA/dec theta (north to south) (0, pi) phi (east to west) (0, 2*pi) declination: north pole = pi/2, south pole = -pi/2 right ascension: (0, 2*pi)

theta = pi/2 - dec ra = phi

property massp

pseudo_cdf_inverse(dec_in=None, ra_in=None, ndraws=1, stype='vecthist')

Select points from the skymap with a distribution following its corresponding pixel probability. If dec_in, ra_in are suupplied, they are ignored except that their shape is reproduced. If ndraws is supplied, that will set the shape. Will return a 2xN numpy array of the (dec, ra) values. stype controls the type of sampling done to retrieve points. Valid choices are ‘rejsamp’: Rejection sampling: accurate but slow ‘vecthist’: Expands a set of points into a larger vector with the multiplicity of the points in the vector corresponding roughly to the probability of drawing that point. Because this is not an exact representation of the proability, some points may not be represented at all (less than quantum of minimum probability) or inaccurately (a significant fraction of the fundamental quantum).

pseudo_pdf(dec_in, ra_in)

Return pixel probability for a given dec_in and ra_in. Note, uses healpy functions to identify correct pixel.

renormalize()

Identify the points contributing to the overall cumulative probability distribution, and set the proper normalization.

static thph2decra(th, ph)

theta/phi to RA/dec theta (north to south) (0, pi) phi (east to west) (0, 2*pi) declination: north pole = pi/2, south pole = -pi/2 right ascension: (0, 2*pi)

dec = pi/2 - theta ra = phi

class RIFT.integrators.mcsampler.MCSampler

Bases: SamplerOutputMixin, object

Class to define a set of parameter names, limits, and probability densities.

add_parameter(params, pdf, cdf_inv=None, left_limit=None, right_limit=None, prior_pdf=None, adaptive_sampling=False)

Add one (or more) parameters to sample dimensions. params is either a string describing the parameter, or a tuple of strings. The tuple will indicate to the sampler that these parameters must be sampled together. left_limit and right_limit are on the infinite interval by default, but can and probably should be specified. If several params are given, left_limit, and right_limit must be a set of tuples with corresponding length. Sampling PDF is required, and if not provided, the cdf inverse function will be determined numerically from the sampling PDF.

cdf_function(param)

Numerically determine the CDF from a given sampling PDF. If the PDF itself is not normalized, the class will keep an internal record of the normalization and adjust the PDF values as necessary. Returns a function object which is the interpolated CDF.

cdf_inverse(param)

Numerically determine the inverse CDF from a given sampling PDF. If the PDF itself is not normalized, the class will keep an internal record of the normalization and adjust the PDF values as necessary. Returns a function object which is the interpolated CDF inverse.

clear()

Clear out the parameters and their settings, as well as clear the sample cache.

draw(rvs, *args, **kwargs)

Draw a set of random variates for parameter(s) args. Left and right limits are handed to the function. If args is None, then draw all parameters. ‘rdict’ parameter is a boolean. If true, returns a dict matched to param name rather than list. rvs must be either a list of uniform random variates to transform for sampling, or an integer number of samples to draw.

integrate(func, *args, **kwargs)

Integrate func, by using n sample points. Right now, all params defined must be passed to args must be provided, but this will change soon.

Limitations:

func’s signature must contain all parameters currently defined by the sampler, and with the same names. This is required so that the sample values can be passed consistently.

kwargs: nmax – total allowed number of sample points, will throw a warning if this number is reached before neff. neff – Effective samples to collect before terminating. If not given, assume infinity n – Number of samples to integrate in a ‘chunk’ – default is 1000 save_integrand – Save the evaluated value of the integrand at the sample points with the sample point history_mult – Number of chunks (of size n) to use in the adaptive histogramming: only useful if there are parameters with adaptation enabled tempering_exp – Exponent to raise the weights of the 1-D marginalized histograms for adaptive sampling prior generation, by default it is 0 which will turn off adaptive sampling regardless of other settings temper_log – Adapt in min(ln L, 10^(-5))^tempering_exp tempering_adapt – Gradually evolve the tempering_exp based on previous history. floor_level – total probability of a uniform distribution, averaged with the weighted sampled distribution, to generate a new sampled distribution n_adapt – number of chunks over which to allow the pdf to adapt. Default is zero, which will turn off adaptive sampling regardless of other settings convergence_tests - dictionary of function pointers, each accepting self._rvs and self.params as arguments. CURRENTLY ONLY USED FOR REPORTING

tripwire_fraction - fraction of nmax of iterations after which n_eff needs to be greater than 1+epsilon for epsilon a small number tripwire_epsilon - small number used in tripwire test

Pinning a value: By specifying a kwarg with the same of an existing parameter, it is possible to “pin” it. The sample draws will always be that value, and the sampling prior will use a delta function at that value.

integrate_vegas(func, *args, **kwargs)

Uses vegas to do the integral. Does not return sample points Remember: pdf, cdf_inv refer to the sampling prior, so I need to multiply the integrand by a PDF ratio product!

static match_params_from_args(args, params)

Given two unordered sets of parameters, one a set of all “basic” elements (strings) possible, and one a set of elements both “basic” strings and “combined” (basic strings in tuples), determine whether the sets are equivalent if no basic element is repeated.

e.g. set A ?= set B

(“a”, “b”, “c”) ?= (“a”, “b”, “c”) ==> True ((“a”, “b”, “c”)) ?= (“a”, “b”, “c”) ==> True ((“a”, “b”), “d”)) ?= (“a”, “b”, “c”) ==> False # basic element ‘d’ not in set B ((“a”, “b”), “d”)) ?= (“a”, “b”, “d”, “c”) ==> False # not all elements in set B represented in set A

RIFT.integrators.mcsampler.M_samp_vector(Mmin, Mmax, x)

exception RIFT.integrators.mcsampler.NanOrInf(value)

Bases: Exception

RIFT.integrators.mcsampler.clip_angle_limits(lo, hi, kind)

Clip an angular range [lo,hi] (radians) to the physical domain of kind (‘declination’ -> [-pi/2,pi/2], ‘inclination’ -> [0,pi]).

Raises ValueError if the requested range is empty/inverted, or if it does not overlap the physical domain. Returns (lo, hi) as floats with lo < hi.

RIFT.integrators.mcsampler.convergence_test_MostSignificantPoint(pcut, rvs, params)

RIFT.integrators.mcsampler.convergence_test_NormalSubIntegrals(ncopies, pcutNormalTest, sigmaCutRelativeErrorThreshold, rvs, params)

RIFT.integrators.mcsampler.cos_samp(x)

RIFT.integrators.mcsampler.cos_samp_cdf_inv_vector(p)

RIFT.integrators.mcsampler.cos_samp_vector(x)

RIFT.integrators.mcsampler.cosine_sampler_limits(lo, hi, kind)

Map an angular range [lo,hi] (radians) into the coordinate actually sampled by RIFT’s ‘cosine’ sky/orientation samplers.

kind=’declination’: sampled variable is z = sin(dec); sin is increasing, so

the limit order is preserved: [lo,hi] -> [sin(lo), sin(hi)].

kind=’inclination’: sampled variable is z = cos(iota); cos is DECREASING, so

the limit order SWAPS: [lo,hi] -> [cos(hi), cos(lo)].

The range is clipped to the physical angular domain first, and the result is clipped to the sampler domain [-1,1]. Raises ValueError on an empty or inverted request. Returns (z_lo, z_hi) with z_lo < z_hi.

RIFT.integrators.mcsampler.dec_samp(x)

RIFT.integrators.mcsampler.dec_samp_cdf_inv_vector(p)

RIFT.integrators.mcsampler.dec_samp_vector(x)

RIFT.integrators.mcsampler.delta_func_pdf(x_0, x)

RIFT.integrators.mcsampler.delta_func_samp(x_0, x)

RIFT.integrators.mcsampler.distance_limit_range(value, d_min, d_max)

Parse and validate a –limit-distance ‘LO,HI’ request (Mpc).

Returns (lo, hi) as floats with d_min <= lo < hi <= d_max. Raises ValueError on an unparseable, non-finite, empty, inverted, or out-of-prior-range request: a box outside [d_min,d_max] would ask the sampler for distances the prior gives zero mass, which is a configuration error, not something to clip silently.

RIFT.integrators.mcsampler.distance_sampler_kwargs(sampler_module, sampling_range, prior_range, d_prior='Euclidean', xpy=<module 'numpy' from '/opt/hostedtoolcache/Python/3.11.18/x64/lib/python3.11/site-packages/numpy/__init__.py'>, adaptive_sampling=False)

Build the add_parameter() kwargs for the distance coordinate.

sampler_module supplies the backend-appropriate uniform sampler helpers (mcsampler / mcsamplerGPU / mcsamplerAdaptiveVolume all export them, and the GPU ones return device arrays). sampling_range sets the sampler’s proposal and its (llim, rlim); prior_range sets the normalization of prior_pdf and is the ONLY thing that determines the reported evidence scale. Pass them equal for the historical behaviour.

RIFT.integrators.mcsampler.gauss_samp(mu, std, x)

RIFT.integrators.mcsampler.gauss_samp_withfloor(mu, std, myfloor, x)

RIFT.integrators.mcsampler.infer_array_module(x, xpy=None)

Return the array module (numpy, cupy, …) that should be used to operate on x.

The truncated samplers below are handed to whichever backend the run selected: mcsampler and mcsamplerAdaptiveVolume call them with host (numpy) arrays, while mcsamplerGPU.draw_simplified() calls self.pdf[param](samples) / self.cdf_inv[param](p) with a SINGLE positional argument holding a cupy array. Defaulting to numpy would then hit numpy.asarray(cupy_array), which raises, so the backend is inferred from the argument instead of assumed. An explicit xpy= always wins.

The lookup is by the array type’s top-level module, taken from sys.modules (a cupy array cannot exist unless cupy is already imported), so this file keeps its numpy-only import list and works for any duck-typed backend exposing the numpy API.

RIFT.integrators.mcsampler.inv_uniform_cdf(a, b, x)

RIFT.integrators.mcsampler.linear_down_samp(x, xmin=0, xmax=1)

distribution p(x) propto (xmax-x)

RIFT.integrators.mcsampler.linear_down_samp_cdf(x, xmin=0, xmax=1)

CDF of distribution p(x) propto (xmax-x)

RIFT.integrators.mcsampler.power_down_samp(x, xmin=0, xmax=1, alpha=3)

distribution p(x) propto (xmax-x)

RIFT.integrators.mcsampler.power_down_samp_cdf(x, xmin=0, xmax=1, alpha=3)

CDF of distribution p(x) propto (xmax-x)

RIFT.integrators.mcsampler.pseudo_dist_samp(r0, r)

RIFT.integrators.mcsampler.q_cdf_inv_vector(qmin, qmax, x)

RIFT.integrators.mcsampler.q_samp_vector(qmin, qmax, x)

RIFT.integrators.mcsampler.quadratic_samp(rmax, x)

RIFT.integrators.mcsampler.ret_cos_samp_cdf_inv_vector(incl_lo, incl_hi)

Inverse CDF (p in [0,1] -> inclination) for uniform-in-cos(iota) truncated to [incl_lo, incl_hi]. Monotonically increasing in p: p=0 gives incl_lo (which is arccos of the UPPER cosine limit – note the swap).

RIFT.integrators.mcsampler.ret_cos_samp_vector(incl_lo, incl_hi)

Sampling pdf in INCLINATION for a uniform-in-cos(iota) draw truncated to [incl_lo, incl_hi]. Normalized to unity over that box. Reduces to cos_samp_vector for the full range.

RIFT.integrators.mcsampler.ret_dec_samp_cdf_inv_vector(dec_lo, dec_hi)

Inverse CDF (p in [0,1] -> declination) for uniform-in-sin(dec) truncated to [dec_lo, dec_hi]. Monotonically increasing in p.

RIFT.integrators.mcsampler.ret_dec_samp_vector(dec_lo, dec_hi)

Sampling pdf in DECLINATION for a uniform-in-sin(dec) draw truncated to [dec_lo, dec_hi]. Normalized to unity over that box (the samplers that use an explicit cdf_inv do not renormalize the pdf themselves). Reduces to dec_samp_vector for the full range.

RIFT.integrators.mcsampler.ret_uniform_samp_vector_alt(a, b)

RIFT.integrators.mcsampler.sanityCheckSamplerIntegrateUnity(sampler, *args, **kwargs)

RIFT.integrators.mcsampler.to_host(x)

Return x as a host (numpy) array, copying it off a device backend if it is on one.

numpy.asarray(cupy_array) RAISES – cupy refuses implicit host conversion – so the copy has to go through the device module’s own asnumpy. The module is looked up from the value, the same way infer_array_module does it, so this file keeps its numpy-only import list. A host input is returned UNCHANGED, not re-wrapped: the numpy path must stay bit for bit what it was.

A non-numpy backend it cannot convert RAISES rather than passing the value through. infer_array_module recognizes any module exposing asarray/where/clip, which is a wider set than the ones exposing asnumpy (torch is in the gap), so “return it unchanged” would hand a device array to the host-only arithmetic below the call site – the silent host/device mix this function exists to remove, reintroduced one backend later.

RIFT.integrators.mcsampler.uniform_samp(a, b, x)

RIFT.integrators.mcsampler.uniform_samp_cdf_inv_vector(a, b, p)

RIFT.integrators.mcsampler.uniform_samp_cos_theta(x)

RIFT.integrators.mcsampler.uniform_samp_dec(x, numpy=<module 'numpy' from '/opt/hostedtoolcache/Python/3.11.18/x64/lib/python3.11/site-packages/numpy/__init__.py'>)

RIFT.integrators.mcsampler.uniform_samp_phase(x, numpy=<module 'numpy' from '/opt/hostedtoolcache/Python/3.11.18/x64/lib/python3.11/site-packages/numpy/__init__.py'>)

RIFT.integrators.mcsampler.uniform_samp_psi(x, numpy=<module 'numpy' from '/opt/hostedtoolcache/Python/3.11.18/x64/lib/python3.11/site-packages/numpy/__init__.py'>)

RIFT.integrators.mcsampler.uniform_samp_theta(x, numpy=<module 'numpy' from '/opt/hostedtoolcache/Python/3.11.18/x64/lib/python3.11/site-packages/numpy/__init__.py'>)

RIFT.integrators.mcsampler.uniform_samp_withfloor_vector(rmaxQuad, rmaxFlat, pFlat, x)

MC Sampler with Adaptive Volume

exception RIFT.integrators.mcsamplerAdaptiveVolume.LiveVolumeCollapse

Bases: Exception

The adaptive-volume live set is empty (or carries no usable information).

Raised INSTEAD of the bare numpy/cupy “zero-size array to reduction operation … which has no identity” that the empty-live-volume path used to produce, so callers and logs can tell a degenerate contraction apart from a waveform generation failure. See the collapse discussion in get_likelihood_threshold.

class RIFT.integrators.mcsamplerAdaptiveVolume.MCSampler(n_chunk=400000, **kwargs)

Bases: SamplerOutputMixin, object

add_parameter(params, pdf, cdf_inv=None, left_limit=None, right_limit=None, prior_pdf=None, adaptive_sampling=False)

Add one (or more) parameters to sample dimensions. params is either a string describing the parameter, or a tuple of strings. The tuple will indicate to the sampler that these parameters must be sampled together. left_limit and right_limit are on the infinite interval by default, but can and probably should be specified. If several params are given, left_limit, and right_limit must be a set of tuples with corresponding length. Sampling PDF is required, and if not provided, the cdf inverse function will be determined numerically from the sampling PDF.

bootstrap_from_fisher(mean, fisher, **kwargs)

Warm-start from a Fisher matrix (mean, Gamma): cov = Gamma^{-1}. This is the ‘Fisher-matrix oracle’ – an essentially free substitute for an expensively-trained flow, giving the integrator a correct-to-2nd-order starting proposal.

bootstrap_from_gaussian(mean, cov, n=None, params=None, enc_prob=0.999, seed=None, cover_frac=0.0, dilate=1)

Warm-start from a single Gaussian proposal N(mean, cov) – the Fisher-oracle entry point. Draws n points from the (box-clipped) Gaussian and builds the live-volume grid from them.

IMPORTANT – this is a UNIMODAL seed. A single Gaussian covers only one mode, and because VARAHA’s live volume only ever contracts, any mode the seed misses is lost forever and biases the integral low. Use this only when the target is (locally) unimodal – e.g. a Fisher matrix at the MAP. For known multimodal structure use bootstrap_from_gaussian_mixture(); for an empirical proposal use bootstrap_from_samples() (both cover the full support and stay unbiased).

cover_frac (0..1): optional safety valve – fraction of the seed cloud drawn uniformly from the full box, trading efficiency for coverage on a possibly-misspecified seed. Default 0.

bootstrap_from_gaussian_mixture(means, covs, weights=None, n=None, params=None, enc_prob=0.999, seed=None, dilate=1)

Warm-start from a MIXTURE of Gaussians – the general oracle seed for multimodal targets. This is what a flow oracle, a GMM fit of a previous posterior, or a set of known degenerate modes (e.g. sky reflections) provides. Because the seed cloud covers every component, the resulting live volume is a superset of the support and the integral stays unbiased.

bootstrap_from_samples(samples, params=None, loglkl=None, enc_prob=0.999, cover_frac=0.0, dilate=1, inflate=1.0, seed=None)

Warm-start from an explicit set of reference points populating the high-likelihood region (e.g. a previous run’s posterior draws, a puff of an earlier MAP point, or fair-draw samples from a prior ILE instance). loglkl (optional) is L*prior at those points, used to seed the threshold.

cover_frac (0..1) mixes this fraction of uniform full-box points into the seed cloud, widening the seeded live volume. Leave 0 when reusing a proposal for the SAME problem (e.g. an in-run second pass).

IT IS NOT A COVERAGE GUARANTEE, despite what this docstring claimed until 2026-08. A FINITE set of uniform points occupies only the bins it lands in, so the seeded grid is NOT a superset of a cold (uniform) start – and the shortfall grows fast with dimension. Measured, fraction of the [-5,5]^d prior box covered by the seeded grid (a cold start is 1.0 by construction):

cover_frac: 0.0 0.2 0.5 0.9 d=2 0.027 0.634 0.982 1.000 d=4 0.0015 0.028 0.104 0.620 d=6 6.3e-05 0.00087 0.0033 0.0287

At d=6 even cover_frac=0.9 leaves 97% of the box unsampled. So the claim that “a warm-started integral can never be MORE biased than a cold one” was false; do not rely on it.

WHAT ACTUALLY PROTECTS YOU is having a component with support everywhere. In the default AV+GMM portfolio that is the GMM member: a Gaussian mixture has nonzero density over the whole box, so q_mix never vanishes and a badly-seeded AV member costs efficiency rather than bias (measured with a deliberately displaced seed at d=4 and d=6: |lnZ bias| <= 0.05 in every run). An ALL-AV portfolio has no such member – every component is a hard-edged box – and the same displaced seed gave lnZ bias -1.0 to -6.8 nats.

And note the limit of any coverage fix: in that all-AV test, keeping one member fully cold (V=1) still gave -1.1 to -4.2 nats, because a uniform member at d=6 finds a sharp peak too rarely to carry the integral within the budget (n_eff 3-9). Coverage in principle is necessary, not sufficient. A badly mismatched seed is an efficiency catastrophe no knob repairs – detect it (the L0 rescue) or do not warm-start across dissimilar points.

inflate (>=1) is the HANDOFF SAFETY MARGIN: widen the seed cloud by this factor about its own mean before building the grid. When importing a proposal that came from a neighbouring intrinsic point, the true peak is shifted (and often slightly broader) at this point, so an un-inflated seed may sit just off it; inflate>1 (e.g. 1.5-2) gives margin for that shift while staying far tighter than a cold start. cover_frac is the coarse safety net for gross mismatch; inflate is the fine margin for a modest shift.

clear()

Clear out the parameters and their settings, as well as clear the sample cache.

draw_simple()

draw_simplified(n_to_get, *args, **kwargs)

has_unbounded_support = False

Class to define a set of parameter names, limits, and probability densities.

integrate(func, *args, **kwargs)

Integrate func, by using n sample points. Right now, all params defined must be passed to args must be provided, but this will change soon. Does NOT allow for tuples of arguments, an unused feature in mcsampler

kwargs: nmax – total allowed number of sample points, will throw a warning if this number is reached before neff. neff – Effective samples to collect before terminating. If not given, assume infinity n – Number of samples to integrate in a ‘chunk’ – default is 1000 save_integrand – Save the evaluated value of the integrand at the sample points with the sample point history_mult – Number of chunks (of size n) to use in the adaptive histogramming: only useful if there are parameters with adaptation enabled tempering_exp – Exponent to raise the weights of the 1-D marginalized histograms for adaptive sampling prior generation, by default it is 0 which will turn off adaptive sampling regardless of other settings temper_log – Adapt in min(ln L, 10^(-5))^tempering_exp tempering_adapt – Gradually evolve the tempering_exp based on previous history. floor_level – total probability of a uniform distribution, averaged with the weighted sampled distribution, to generate a new sampled distribution n_adapt – IGNORED as an adaptation schedule by this sampler; accepted only for API

compatibility with mcsampler/mcsamplerGPU, where it does gate update_sampling_prior. AV has no update_sampling_prior: its volume adaptation is intrinsic to the algorithm and runs every cycle regardless of this value. (Verified: output is bit-identical for n_adapt in {10,100,1000}.) The ONLY thing it still does here is participate in the save_intg gate, so n_adapt=0 can still suppress the _rvs cache. Do not reach for this expecting an “adapt then freeze” control – there isn’t one.

convergence_tests - dictionary of function pointers, each accepting self._rvs and self.params as arguments. CURRENTLY ONLY USED FOR REPORTING Pinning a value: By specifying a kwarg with the same of an existing parameter, it is possible to “pin” it. The sample draws will always be that value, and the sampling prior will use a delta function at that value.

integrate_log(lnF, *args, xpy=<module 'numpy' from '/opt/hostedtoolcache/Python/3.11.18/x64/lib/python3.11/site-packages/numpy/__init__.py'>, **kwargs)

Integrate exp(lnF) returning lnI, by using n sample points, assuming integrand is lnF Does NOT allow for tuples of arguments, an unused feature in mcsampler

tempering is done with lnF, suitably modified.

kwargs: nmax – total allowed number of sample points, will throw a warning if this number is reached before neff. neff – Effective samples to collect before terminating. If not given, assume infinity n – Number of samples to integrate in a ‘chunk’ – default is 1000 save_integrand – Save the evaluated value of the integrand at the sample points with the sample point history_mult – Number of chunks (of size n) to use in the adaptive histogramming: only useful if there are parameters with adaptation enabled tempering_exp – Exponent to raise the weights of the 1-D marginalized histograms for adaptive sampling prior generation, by default it is 0 which will turn off adaptive sampling regardless of other settings temper_log – Adapt in min(ln L, 10^(-5))^tempering_exp tempering_adapt – Gradually evolve the tempering_exp based on previous history. floor_level – total probability of a uniform distribution, averaged with the weighted sampled distribution, to generate a new sampled distribution n_adapt – IGNORED as an adaptation schedule by this sampler; accepted only for API

compatibility with mcsampler/mcsamplerGPU, where it does gate update_sampling_prior. AV has no update_sampling_prior: its volume adaptation is intrinsic to the algorithm and runs every cycle regardless of this value. (Verified: output is bit-identical for n_adapt in {10,100,1000}.) The ONLY thing it still does here is participate in the save_intg gate, so n_adapt=0 can still suppress the _rvs cache. Do not reach for this expecting an “adapt then freeze” control – there isn’t one.

convergence_tests - dictionary of function pointers, each accepting self._rvs and self.params as arguments. CURRENTLY ONLY USED FOR REPORTING Pinning a value: By specifying a kwarg with the same of an existing parameter, it is possible to “pin” it. The sample draws will always be that value, and the sampling prior will use a delta function at that value.

joint_p_s_is_normalized_density = False

load_state(path)

Restore a live-volume state saved by save_state(). Verifies the parameter names and box match this sampler before warm-starting.

static match_params_from_args(args, params)

Given two unordered sets of parameters, one a set of all “basic” elements (strings) possible, and one a set of elements both “basic” strings and “combined” (basic strings in tuples), determine whether the sets are equivalent if no basic element is repeated. e.g. set A ?= set B (“a”, “b”, “c”) ?= (“a”, “b”, “c”) ==> True ((“a”, “b”, “c”)) ?= (“a”, “b”, “c”) ==> True ((“a”, “b”), “d”)) ?= (“a”, “b”, “c”) ==> False # basic element ‘d’ not in set B ((“a”, “b”), “d”)) ?= (“a”, “b”, “d”, “c”) ==> False # not all elements in set B represented in set A

prior_prod(x)

Evaluates prior_pdf(x), multiplying together all factors

prior_pdf are host (numpy) functions in general, so evaluate them on a CPU copy of the samples and convert the product back to the active backend. identity_convert / identity_convert_togpu are no-ops when cupy is not in use.

reset_selfish_state()

Start a new integration pass: forget the live set and threshold, keep the grid.

A new pass may use a different integrand (e.g. the calmarg burn-in followed by the full likelihood), so a threshold carried from the previous pass could exclude every new draw. If the previous pass had reached the final threshold, the grid is frozen as adapted for the new pass: the selfish step then does nothing, so the grid is neither contracted again nor re-binned. Otherwise contraction restarts from the current grid.

sampling_density(X)

Pointwise sampling density q(theta) of THIS member, evaluated at ARBITRARY points X (shape (N, ndim), columns in self.params_ordered order). Returns a host (numpy) array of length N, or None if the live-volume state has not been set up yet.

VARAHA draws uniformly over its live volume – the union of the currently-occupied hypercubes (self.binunique), each of width self.dx. The density of the CURRENT grid is therefore the constant

q_live = 1 / (n_occupied_bins * prod(dx))

inside the live volume and 0 outside it. (integrate_log’s own log_joint_s_prior is the mixture of every grid it drew from; see log_retained_density.) We use the geometric form 1/(n_bins*prod(dx)) directly: it is exactly the density of the points draw_simple() produces (equal draws per occupied bin, uniform within a bin), so a multiple-importance-sampling denominator built from it is unbiased regardless of any drift between the tracked scalar V and the actual bin grid.

This method is READ-ONLY – it does not touch any sampler state and does not affect this sampler’s own integrate_log. It exists so the portfolio can form the balance-heuristic mixture density q_mix = sum_m w_m q_m.

save_state(path)

Serialize the compact live-volume state (occupied bins + widths + volume + threshold) to a lightweight .npz. This is RIFT’s cheap alternative to persisting a trained flow: the entire adapted proposal is just an integer bin-index array plus a few scalars.

setup(**kwargs)

update_sampling_prior_selfish(lnF, *args, xpy=<module 'numpy' from '/opt/hostedtoolcache/Python/3.11.18/x64/lib/python3.11/site-packages/numpy/__init__.py'>, no_protect_names=True, **kwargs)

update_sampling_prior

Update VARAHA sampling hypercubes/

Note that external samples are NOT uniform. VARAHA should only be trained on its own samples, not others!

We therefore do a single pure step of VARAHA, including independent draws. We will keep state about ‘V’ etc from previous iterations.

We therefore also have to know about the function we are integrating. However, we do not keep track of the integral result here – the top -level routine does this.

One call is one cycle of integrate_log’s loop, so the cycle state persists between calls in self._selfish_state: the live set (allx/allloglkl/allp), the threshold and the truncated probability trunc_p. Restarted on every call, each call would keep only the top nsel of its own draws, V would shrink by n_chunk/nsel per call, the final threshold would never be reached, and the live volume would contract without limit. The state is cleared by setup(), by installing a warm seed, and by reset_selfish_state(), which mcsamplerPortfolio calls at the start of each pass. The carried live set is capped (_cap_selfish_live_set) to bound memory and cost.

warm_seed_axes()

Column indices a warm seed must span: the ADAPTIVE axes (all of them when nothing is adaptive, since then the grid is one bin per dim and the seed’s only job is to be well-defined). Exposed so a caller building a seed – the ILE’s L0 rescue – can ask the sampler which dimensions its seed will be judged on instead of guessing. A portfolio has no such axes of its own; ask a member.

exception RIFT.integrators.mcsamplerAdaptiveVolume.NanOrInf(value)

Bases: Exception

RIFT.integrators.mcsamplerAdaptiveVolume.build_warm_seed(points, lnL, box_lo, box_hi, axes, deltalnL=15.0, puff_width_frac=0.005, puff_scale='auto', puff_factor=2.0, n_puff=2000, seed=None)

Build the L0 rescue’s warm seed from a pass’s own samples -> (seed, info).

points (n, ndim) and lnL (n,) are the completed pass’s draws. The seed is the points within deltalnL of the peak, PUFFED to full rank if they do not span axes.

RANK, NOT COUNT, is the guard. The rule this replaces was len(seed) < 2, and a count cannot see the failure: measured on zero-noise injections, a 5-point seed at rho_net 102.8 had affine rank 2-4 of 6 and a 2-point seed at rho_net 146.8 had rank 0 of 6. Both passed the count test, and both then warm-started a live volume that had collapsed onto a degenerate subspace (V ~ 3e-06 and ~9e-36 against a healthy ~1e-08), which reports a fine n_eff while lnZ is a lower bound. n points span at most n-1 affine dimensions, so the rank test subsumes the count it replaces.

AUGMENT, DO NOT REPLACE. The handful of real points are the only direct evidence of where the peak is and how wide it is, so they are kept and the puff is added alongside them. This can only help: the grid is built from the union’s extent, so a real point lying outside the puff widens the seeded volume to include it, and VARAHA can only contract afterwards – whereas replacing them throws that information away and pins the support to a guessed width about a single point.

puff_scale:

‘fixed’ – isotropic, puff_width_frac of each parameter’s prior range (the

historical behaviour; puff_width_frac = 1/200 reproduces it exactly).

‘auto’ – the measured posterior scale and correlations from every finite lnL the

pass drew (warm_seed_scale_from_finite_points), falling back to ‘fixed’ when there are too few finite points to estimate one.

puff_factor multiplies the resulting width (variance scales as its square). 2 is the measured optimum and both tails are wrong – see warm_seed_scale_from_finite_points.

RIFT.integrators.mcsamplerAdaptiveVolume.dilate_bins(bins, nbins, axes, layers=1)

Occupied bins plus their axis neighbours (layers steps along each axis in axes), clipped to the grid. The occupied set is built from retained points, so it misses the part of the live region that holds no point yet; the neighbour layer keeps that part drawable.

RIFT.integrators.mcsamplerAdaptiveVolume.get_likelihood_threshold(lkl, lkl_thr, nsel, discard_prob, xpy_here=<module 'numpy' from '/opt/hostedtoolcache/Python/3.11.18/x64/lib/python3.11/site-packages/numpy/__init__.py'>)

Find the likelihood threshold that encolses a probability lkl : array of likelihoods (on bins) lkl_thr: scalar cutoff nsel : integer, has to do with size of array of likelihoods used to evaluate for next array. discard_prob: threshold on CDF to throw away an entire bin. Should be very small

RIFT.integrators.mcsamplerAdaptiveVolume.identity_convert(x)

RIFT.integrators.mcsamplerAdaptiveVolume.identity_convert_togpu(x)

RIFT.integrators.mcsamplerAdaptiveVolume.keep_reserve_from_rvs(sampler, label, integrand_is_log=None)

Attach sampler._warm_seed_reserve from its own _rvs, or leave it None.

Provenance for an export, never a reason to lose a completed integral – the same contract AV and the portfolio already use at their own reserve sites.

RIFT.integrators.mcsamplerAdaptiveVolume.live_volume_collapse_verdict(n_live, ndim, ess=None, khat=None, n_empty_cycles=0, n_live_collapses=0, n_warm_seed=None, n_warm_seed_rank=None, n_warm_seed_dim=None)

Has the adaptive-volume live set degenerated? -> (collapsed, [reasons])

A degenerate contraction must be REPORTED rather than silently exported: the run still returns a lnZ and a sample cloud, but both describe a single mode and the cloud is not a fair posterior draw.

Thresholds, against the separation measured on zero-noise injections at a fixed intrinsic point (rho_net 51 -> 147):

healthy ESS 16.8-20.2, k-hat 1.03-1.50 (rho 51.4, converges cold) collapsed ESS 1.0-1.7, k-hat 21-202 (rho 103-147)

• n_live <= ndim – a live set no larger than the dimension cannot span the space, let alone describe a posterior in it. Geometric, not tuned.

• ESS < 2 – fewer than two effective samples IS one sample.

• ESS < 5 with k-hat > 10 – near-degenerate AND a pathological weight tail.

The gap between the regimes is an order of magnitude wide, so these sit far from both sides of it. k-hat is deliberately NOT a gate on its own: it exceeds its nominal 0.70 “unresolved tail” threshold even in the healthy runs on this problem, and it is not always computable once the live set is tiny.

A WARM start fails the OTHER way, and none of the rules above see it: seeded from too few points the grid contracts onto a sliver of the support, the integrand is then flat across it, and the pass terminates in one cycle looking excellent – large n_live, ESS ~ n, small k-hat – while lnZ is short by the mass outside the sliver. Measured on 12 rho_net=146.8 rescue replicates: the eleven seeded from 2000 puffed points warm-started at V = 7.5e-9 to 1.5e-8 (351-684 live bins) and returned ln(Z/Lmax) = -27.0 to -30.6; the one seeded from 2 points warm-started at V = 9.2e-36 (13 bins) and returned -80.7, i.e. ~50 nats low, with eff_samp 9789 of 10010 samples. So:

• n_warm_seed_rank < n_warm_seed_dim – the seed cloud’s AFFINE RANK (rank of the mean-centred points, per-axis scaled by the box) is below the dimension it must span, so it lies in a lower-dimensional subspace and cannot define a volume there.

RANK, not row count, is the invariant. Rows are neither necessary nor sufficient: thousands of duplicated or collinear points span the same degenerate subspace two points do and fail identically, while d+1 affinely independent points are exactly enough to define a volume in d dimensions and must NOT be flagged. Rank subsumes the count anyway, since n points span at most n-1 affine dimensions.

n_warm_seed* are None on a cold pass -> the rule is skipped. If rank is unavailable (a grid restored by load_state from a run that predates it) we fall back to the count, at the correct simplex boundary: fewer than dim+1 points, i.e. n <= dim.

RIFT.integrators.mcsamplerAdaptiveVolume.lnZ_from_reserve(reserve)

lnZ implied by a warm-seed reserve, on the ORIGINAL PROPOSAL-DRAW normalization.

-> float, or None if the reserve cannot support the estimate.

THE DENOMINATOR IS THE POINT. An importance estimate is mean(w) over the draws that were MADE, and a draw whose likelihood underflowed to -inf is a real draw contributing a real zero. make_warm_seed_reserve drops those rows – they are ballast for every other consumer – and may then cap what survives, so the stored array is neither the draw set nor a uniform sample of it. Averaging over the stored rows silently divides by the wrong n, overestimating by log(n_retained/n_finite). For a PORTFOLIO, whose _rvs holds every draw and whose finite fraction on a collapsed pass is ~1e-5, that is ~11 nats.

It does not cancel in the L0 reject gate, which is what makes it dangerous rather than merely wrong: the cold and warm passes have DIFFERENT finite fractions, so the gate would see the difference of two different-sized errors and could again keep a mass-losing warm pass or discard a good one – the failure the gate change exists to remove, reappearing one layer down.

So rebuild the estimate explicitly. With m rows kept, out of n_finite finite draws, out of n_retained draws made:

sum over the finite population ~= (n_finite / m) * sum over kept rows lnZ = logsumexp(lw_kept) + log(n_finite) - log(m) - log(n_retained)

which collapses to logsumexp(lw) - log(m) exactly when nothing was dropped and nothing was capped – i.e. for AV, whose retained set is already all finite, leaving its behaviour unchanged.

One wrinkle, disclosed rather than corrected: make_warm_seed_reserve force-appends the peak row, so the kept rows are very slightly peak-biased rather than uniform. At one row in n_max (default 20000) that sits far below the MC error on either pass, and correcting it exactly would mean tracking whether the peak was already in the subsample.

RIFT.integrators.mcsamplerAdaptiveVolume.log_retained_density(X, origin, grids, box_lo, dx0, pinned_dims=())

Per-point log density of the retained set, for the deterministic-mixture estimator.

grids[k] = (dx, bins, n_drawn) is the proposal integrate_log drew cycle k from: uniform over the occupied bins (an (M, ndim) array, or a _BinSet over the non-pinned dims), density 1/|B_k|. The pooled draws come from the mixture sum_k n_k 1[x in B_k]/|B_k|; dividing by the retained count gives a density for which logsumexp(lnL + ln p - ln p_s) - ln(n_retained) is the balance-heuristic estimate (1/N) sum L p / q_mix. origin[i] is the cycle that drew point i, which is a member of its own grid by construction (guards rounding at a bin’s upper edge). Pinned dimensions are not drawn from the bins, so they enter with their full width.

RIFT.integrators.mcsamplerAdaptiveVolume.make_reserve_from_rvs(rvs, params_ordered, n_max=20000, integrand_is_log=None, convert=None)

The retained-set reserve for a sampler that still has its rows in _rvs -> dict.

ONE adapter in front of make_warm_seed_reserve, for the backends that reach the fair draw with the retained set still in _rvs: mcsampler, mcsamplerGPU (both entry points), mcsamplerEnsemble and mcsamplerNFlow. AV builds its own from allx/allloglkl because it never keeps them in _rvs; the portfolio aggregates its members’ rows.

WHY THESE BACKENDS NEED ONE AT ALL. Until now the reserve existed only for the L0 rescue, so only AV and the portfolio built it and reserve=None elsewhere was described as “the honest answer, not a gap”. It is a gap for anything that EXPORTS a shape: the fair draw below replaces _rvs with min(–fairdraw-extrinsic-output-n-max (5 by default), 1.5*eff_samp, 1.5*neff) rows taken WITH REPLACEMENT, and the .dgrid exporter then binned those five rows as though they were a likelihood-vs-distance curve. The retained rows are right here; keeping a bounded copy of them costs one subsample and lets the exporter use the sample set instead of its export resample.

BOTH CONVENTIONS. integrand is lnL on some backends and linear L on others, and the sampler is the only thing that knows which (see RvsRecord.integrand_is_log). The log columns are preferred where they exist; otherwise the linear ones are converted here. log(0) is -inf, which make_warm_seed_reserve drops as a zero-weight row.

RIFT.integrators.mcsamplerAdaptiveVolume.make_warm_seed_reserve(X, lnL, params_ordered, n_max=20000, log_joint_prior=None, log_joint_s_prior=None, rng=None, force_peak=True)

A bounded copy of the points a pass RETAINED, for a later warm start -> dict.

THE ONE BUILDER, because every sampler that can be L0-rescued needs the identical record and they reach this point by different routes: mcsamplerAdaptiveVolume from its own accumulated draws, mcsamplerPortfolio from its aggregated _rvs (it drives members through draw_simplified(), never their integrate_log(), so a member never builds one).

WHY IT HAS TO BE TAKEN EARLY. Both samplers then prune _rvs and fair-draw it down to ~1.5*n_eff rows resampled WITH REPLACEMENT – a resample built for EXPORT. Anything that reads _rvs afterwards and treats it as the sample set sees, on the collapsed pass a rescue exists for, a handful of rows several of which are the same point twice.

STRATIFY BY FINITE-NESS BEFORE BOUNDING, or the bound destroys exactly what it is keeping. AV arrives here with only retained (finite) rows, but the PORTFOLIO’s _rvs holds EVERY draw – and on the collapsed pass this exists for, essentially all of them are -inf. A uniform subsample over all rows then keeps the finite ones in proportion, which is to say almost none: measured, 10 finite rows among 1,000,000 at a 20,000 cap survived as 2 (the forced peak, plus one lucky draw). build_warm_seed would then see a rank-0 core, warm_seed_scale_from_finite_points would decline for want of points, and the puff would fall back to the fixed prior fraction – the very truncation (0.8-8.3 nats, with a healthy-looking ESS) this whole change exists to remove.

So the non-finite rows are dropped outright. They are ballast for every consumer: the seed core is lnL > max - deltalnL, the scale estimator filters to finite, and _lnZ_of_rvs filters to finite before it averages – so dropping them here leaves that lnZ bit-identical while making the cap mean what it says.

NO SECOND STRATUM, deliberately. Once the finite rows are the population, the cap only binds on a run with >n_max FINITE samples – a healthy one, where the peak window is proportionally represented anyway (measured: a 20,001-row reserve out of 4e6 gave a 239-point seed core). Adding a top-lnL stratum on top of that would buy nothing there and would bias both the covariance estimate and any lnZ taken from the reserve, because neither is a uniform sample of the level set any more. The one deliberate exception is the PEAK row, appended unconditionally because the seed is defined relative to it and a subsample can drop it; at one row in n_max its effect on either estimate is negligible.

force_peak IS FOR A SEED, NOT FOR AN EXPORT. Appending the peak row unconditionally costs a seed nothing (“at one row in n_max its effect on either estimate is negligible”) because a seed is a SET of points and the peak defines its centre. A consumer that reads the reserve as a WEIGHTED SAMPLE sees something else. Uniform subsampling without replacement gives every row the same inclusion probability, which cancels in a normalized histogram; that histogram is a RATIO of two unbiased estimators, so it carries an O(1/n) bias rather than none (measured at +1.5 to +3.9 sigma in the small-mass bins over 4000 seeds, i.e. negligible beside the MC noise, but not zero). Forcing one row in at probability 1 while the rest come in at n_max/n_finite is a different thing entirely: a bias that does not shrink with n.

HOW BIG DEPENDS ENTIRELY ON HOW DOMINANT THE PEAK IS, so the numbers below are a measured range and not a constant. On 8000 rows capped at 800, median of 15 subsample seeds, as the top weight is lifted above the rest: no dominant row, ESS 506 forced against 506 free (no effect at all); 3 e-folds up, 223 against 505; 6 e-folds, 3.6 against 505. On a real starved GMM extrinsic pass (180,000 retained rows capped at 20,000, 20 seeds) the forced version reported ESS 33 where the population has 1163, and displaced the exported curve by up to 1.22 nats in a bin against 0.72 free.

force_peak therefore stays True everywhere, and the .dgrid exporter divides that row by its own inclusion probability instead (reserve_distance_and_ln_weights), which is why n_subsample is recorded below. Dropping the row instead is worse than it sounds: the unforced subsample then usually misses the dominant row entirely and reports ESS 506 where the population has 99, against 87 for the forced-and-corrected version.

Neither setting makes the SUBSAMPLE’s own ESS a description of the population – that is what ess_finite, recorded below before the cap, is for.

ISOLATED RNG. This reserve is built unconditionally – including when –sampler-warmstart-retry-neff is unset and nothing will ever read it – so drawing the subsample from the global numpy stream would advance it before the fair draw, before the exported posterior, and before every later event and replica. An opt-in rescue that is switched OFF must not change a seeded run’s output. Default to a private, deterministic generator so the reserve is itself reproducible and costs the caller nothing.

The two prior components ride along when given, so a consumer can rebuild the importance weight – and therefore lnZ – from the retained set rather than from the fair draw.

RIFT.integrators.mcsamplerAdaptiveVolume.profile(fn)

RIFT.integrators.mcsamplerAdaptiveVolume.sample_from_bins(xrange, dx, bu, ninbin, reject_out_of_range=False)

RIFT.integrators.mcsamplerAdaptiveVolume.seed_affine_rank(pts, box_lo, box_hi, axes=None, tol=1e-09)

Affine rank of a warm-seed cloud over axes, i.e. the dimension of the subspace the seed actually spans -> (rank, n_in_box).

THE ONE PLACE this is defined, because two callers must agree exactly: the grid builder records it for the collapse diagnostic, and the ILE’s L0 rescue tests it to decide whether the seed needs puffing. A rescue that puffed on a rank the diagnostic then measured differently would either puff a healthy seed or ship a flagged one.

Measured the way _build_grid_from_points must see it:

• IN-BOX ROWS ONLY. The grid only ever spans the box, so out-of-box rows describe nothing it will build – and left in they inflate the rank, so a seed that is degenerate where it matters could be recorded full-rank.

• mean-centred (AFFINE rank: n points span at most n-1 affine dimensions, so this subsumes the row-count test that used to stand in for it), and

• per-axis scaled by the box, so the tolerance is unit-free – a distance in Mpc and an angle in radians must not get different tolerances.

RIFT.integrators.mcsamplerAdaptiveVolume.set_xpy_to_numpy()

RIFT.integrators.mcsamplerAdaptiveVolume.unseen_fraction(X, lnw, box_lo, dx, axes)

Good-Turing estimate of the posterior mass the grid does not cover: the weight exp(lnw) of live points alone in their bin, as a fraction of the total, binning only axes.

Grids hold only bins with a live point and retained points stay in every later grid, so log_retained_density over-counts where points remain when the live set leaves bins empty. Diagnostic only; see DESIGN_av_coverage_diagnostic.md.

RIFT.integrators.mcsamplerAdaptiveVolume.warm_seed_scale_from_finite_points(points, lnL, box_lo, box_hi, axes, eig_lo=1e-05, eig_hi=0.5)

Estimate the POSTERIOR scale (a box-scaled covariance over axes) from the finite log-likelihoods a collapsed pass already drew -> cov, or None if it cannot be measured.

Why this exists. The L0 rescue’s fallback puff used a hardcoded 1/200 of each parameter’s prior range, which is a property of the PRIOR and knows nothing about the posterior – but the posterior narrows as 1/rho, so one fixed fraction cannot be right across the amplitude range, and a puff narrower than the posterior truncates real mass (VARAHA’s live volume only ever contracts, so the seed is a ceiling on the support).

The information needed is already in hand and was being thrown away. A cold pass at high amplitude draws (near enough) uniformly from the prior box, and returns a finite lnL only inside the region where exp() has not underflowed – i.e. the level set lnL > lnL_max - D. For a locally Gaussian peak that level set is the ellipsoid u^T A u < 2D (u box-scaled about the peak), and points uniform in an ellipsoid have covariance (2D/(d+2)) A^{-1}. So the posterior covariance A^{-1} is recovered as

cov_post = cov(finite points) * (d + 2) / (2 D), D = lnL_max - min(finite lnL)

which is one sample covariance, no fit and nothing to fail to converge. It also delivers the CORRELATIONS – sky position, time and distance are strongly correlated at high amplitude, and an isotropic puff wastes almost all of its points off the ridge.

Deliberately approximate – the draws are only uniform-in-prior until AV starts contracting, and the peak is only locally Gaussian. Measured against a known lnZ on a correlated 6-D peak with the same 745-nat underflow (sigma recovered / true, per axis, 6 replicates): 1.01 - 1.23. So it is good to ~20%, which is what matters, because the error that was being made is a factor of ~5-10.

NEITHER DIRECTION IS FREE, so do not treat “wide is safe” as a licence. Too narrow silently truncates: on that same target a puff at the historical 1/200 of the prior range came in 0.8 - 8.3 nats below the truth, with a healthy-looking ESS. Too wide is not merely inefficient, which is what the surrounding code used to assume. Scanning a multiplier on this estimate, mean (worst) lnZ error over 6 replicates, mean ESS:

x0.5 -8.52 (-19.0) nats, ESS 71 truncated x1 -1.61 (-3.5) nats, ESS 53 x2 +0.08 (-0.2) nats, ESS 52 <– the default x3 +1.13 (+0.6) nats, ESS 26 x6 +3.03 (-1.6) nats, ESS 10 biased HIGH, and efficiency is going x12 -29.99 (-71.9) nats, ESS 1 a cold start in all but name: re-collapses

Both tails are wrong, and the useful range is under a decade wide, so inflate by a small factor and not by an order of magnitude.

Eigenvalues are floored/capped in box-scaled units (eig_lo, eig_hi are standard deviations as a fraction of the box) so no direction can come back degenerate – a zero-width direction would re-create the rank deficiency this is being used to repair.

Ensemble MC Sampler

class RIFT.integrators.mcsamplerEnsemble.HealPixSampler(skymap, massp=1.0)

Bases: object

static decra2thph(dec, ra)

property massp

pseudo_cdf_inverse(dec_in=None, ra_in=None, ndraws=1, stype='vecthist')

pseudo_pdf(dec_in, ra_in)

renormalize()

static thph2decra(th, ph)

class RIFT.integrators.mcsamplerEnsemble.MCSampler

Bases: SamplerOutputMixin, object

add_parameter(params, pdf=None, cdf_inv=None, left_limit=None, right_limit=None, prior_pdf=None, adaptive_sampling=False)

bootstrap_from_samples(samples, params=None, n_comp_warm=2, **kwargs)

Warm-start: fit this GMM’s proposal to a seed cloud so it samples AT the peak from the first draw, instead of having to discover the peak location from scratch. For a needle-in-a-haystack extrinsic posterior (peak ~ 10^-11 of the prior box) this is the difference between converging and never finding the peak by cold draws.

Builds the integrator if setup() has not run yet, then (re)fits one GMM per dim-group from the seed with EQUAL weights – this only shapes the proposal, carries no lnL information, so it cannot bias the estimate (the importance weights still use the true likelihood). AV-only kwargs (cover_frac / inflate) are accepted and ignored, so a portfolio can forward a single seed to every member uniformly.

samples: (N, ndim) array, columns in self.params_ordered order.

calc_pdf(samples)

clear()

draw_simplified(n, *args, **kwargs)

evaluate(samples)

property has_unbounded_support

Does this member’s proposal genuinely have support across the WHOLE prior box?

mcsamplerPortfolio uses this to decide whether it must hold a member cold on a warm start. Getting it wrong costs coverage silently, so it is answered from the INSTALLED MODELS, not from a configuration value.

Two traps this avoids:

• gmm_defensive_frac > 0 is only a REQUEST. add_defensive_component() is called by fit_gmm_adaptive, but the fixed-component fit paths did not call it, and gmm_adaptive defaults to None (off) – so the default configuration asked for a defensive component and never installed one.

• gmm.score() FLOORS its return at 1e-300, so a member always looks like it has nonzero density everywhere. That is a numerical guard against log(0), not coverage: a sample landing there carries weight L*p/q ~ 1e300 and would wreck the estimate rather than support it. Measured, a fixed-component fit to a tight cloud returns exactly that floor at the far corner for every d >= 4.

Reports False whenever it cannot be verified – before the integrator exists, before any group has been trained, or if ANY trained group lacks the component.

integrate(func, *args, **kwargs)

integrate_log(func, *args, **kwargs)

joint_p_s_is_normalized_density = True

static match_params_from_args(args, params)

sampling_density(X)

Pointwise sampling density q(theta) of THIS GMM member, evaluated at ARBITRARY points X (shape (N, ndim), columns in self.params_ordered order). Returns a host (numpy) array of length N, or None if the integrator/GMM has not been built yet.

This is exactly the per-sample product MonteCarloEnsemble._sample stores as sampling_prior_array, but evaluated at supplied points rather than at the member’s own draws: for each grouped set of dimensions it is the fitted mixture density gmm.score(…) (already normalized to integrate to 1 over the box, in ORIGINAL coordinates), or the uniform density 1/vol for a not-yet-fitted (None) group. READ-ONLY; does not affect this sampler’s own integrate(). Used by the portfolio balance heuristic.

setup(n_comp=None, **kwargs)

Build the integrator. REMEMBERS its arguments and re-applies them on later calls.

setup() is not called once: bootstrap_from_samples() re-runs it to rebuild the proposal as a single full-dim group, and mcsamplerPortfolio replays it to reset a member between points. Each of those rebuilt the integrator from ONLY the kwargs of that call, so every option the caller set originally was silently dropped. That has now bitten three separate settings – gmm_dict (the dimension grouping and any seeded models), gmm_defensive_frac, and gmm_defensive_all_paths – each found as its own P1, each patched individually, and a fourth would have followed.

So: merge this call’s kwargs OVER the remembered ones, and remember the result. An explicit argument still wins; an omitted one keeps whatever it was configured to be instead of reverting to a library default. Pass setup_forget=True to start clean.

update_sampling_prior(ln_weights, n_history, tempering_exp=1, log_scale_weights=True, floor_integrated_probability=0, external_rvs=None, **kwargs)

RIFT.integrators.mcsamplerEnsemble.M_samp_vector(Mmin, Mmax, x)

exception RIFT.integrators.mcsamplerEnsemble.NanOrInf(value)

Bases: Exception

RIFT.integrators.mcsamplerEnsemble.complete_dim_group_cover(gmm_dict, ndim, comp_dict=None, gmm_adapt=None, n_comp_default=1, adapt_default=False)

Give every dimension not already in a dim-group its own uncorrelated 1-d group.

MonteCarloEnsemble._sample() allocates sample_array with xpy.empty and writes only the columns a dim-group names, so a dimension in NO group is uninitialized memory – handed to the integrand and to the prior with no matching factor in sampling_prior_array. That is an unnormalized estimator, and because the value is whatever the allocator returned it is not even reproducibly wrong. Returns the keys added, so a caller can log them.

RIFT.integrators.mcsamplerEnsemble.convergence_test_MostSignificantPoint(pcut, rvs, params)

RIFT.integrators.mcsamplerEnsemble.convergence_test_NormalSubIntegrals(ncopies, pcutNormalTest, sigmaCutRelativeErrorThreshold, rvs, params)

RIFT.integrators.mcsamplerEnsemble.cos_samp(x)

RIFT.integrators.mcsamplerEnsemble.cos_samp_cdf_inv_vector(p)

RIFT.integrators.mcsamplerEnsemble.dec_samp(x)

RIFT.integrators.mcsamplerEnsemble.dec_samp_cdf_inv_vector(p)

RIFT.integrators.mcsamplerEnsemble.delta_func_pdf(x_0, x)

RIFT.integrators.mcsamplerEnsemble.delta_func_samp(x_0, x)

RIFT.integrators.mcsamplerEnsemble.dim_group_frame(sampler, integrate_args, driven_by_integrate)

The parameter-name list that numbers this sampler’s gmm_dict dimension indices.

mcsamplerEnsemble numbers its dimensions two DIFFERENT ways depending on how it is driven, and a caller that builds keys in the wrong one silently groups the wrong parameters:

• standalone: integrate(func, *args) builds raw_bounds from args, so dim i is args[i].

• portfolio member: the portfolio never calls member.integrate(). It calls member.setup(…) – where _setup_impl uses dim = len(self.params_ordered) – and then member.draw_simplified(n, *self.params_ordered). So dim i is params_ordered[i].

Pass driven_by_integrate=False for a portfolio member. Returns a list of names; index into it to build a key, and len() it to know how many dimensions must be covered.

RIFT.integrators.mcsamplerEnsemble.gauss_samp(mu, std, x)

RIFT.integrators.mcsamplerEnsemble.gauss_samp_withfloor(mu, std, myfloor, x)

RIFT.integrators.mcsamplerEnsemble.identity_convert(x)

RIFT.integrators.mcsamplerEnsemble.identity_convert_togpu(x)

RIFT.integrators.mcsamplerEnsemble.inv_uniform_cdf(a, b, x)

RIFT.integrators.mcsamplerEnsemble.pseudo_dist_samp(r0, r)

RIFT.integrators.mcsamplerEnsemble.q_cdf_inv_vector(qmin, qmax, x)

RIFT.integrators.mcsamplerEnsemble.q_samp_vector(qmin, qmax, x)

RIFT.integrators.mcsamplerEnsemble.sanityCheckSamplerIntegrateUnity(sampler, *args, **kwargs)

RIFT.integrators.mcsamplerEnsemble.validate_dim_group_cover(gmm_dict, ndim)

Raise ValueError unless gmm_dict partitions range(ndim) exactly.

Public entry point for a DRIVER that wants to fail before it starts work. The same check runs inside the integrator, but a driver that wraps its per-event call in a broad except can turn that refusal into an empty output file and a zero exit status, which is the one outcome this whole guard exists to prevent.

GPU MC Sampler

class RIFT.integrators.mcsamplerGPU.MCSampler

Bases: SamplerOutputMixin, object

Class to define a set of parameter names, limits, and probability densities.

add_parameter(params, pdf, cdf_inv=None, left_limit=None, right_limit=None, prior_pdf=None, adaptive_sampling=False)

Add one (or more) parameters to sample dimensions. params is either a string describing the parameter, or a tuple of strings. The tuple will indicate to the sampler that these parameters must be sampled together. left_limit and right_limit are on the infinite interval by default, but can and probably should be specified. If several params are given, left_limit, and right_limit must be a set of tuples with corresponding length. Sampling PDF is required, and if not provided, the cdf inverse function will be determined numerically from the sampling PDF.

cdf_function(param)

Numerically determine the CDF from a given sampling PDF. If the PDF itself is not normalized, the class will keep an internal record of the normalization and adjust the PDF values as necessary. Returns a function object which is the interpolated CDF. NOT USED IN THIS ROUTINE, SHOULD NOT BE CALLED

cdf_function_from_histogram(x)

Computes the CDF from a histogram at the points x. Params —— x : array_like, shape = sample_shape Returns ——- P(x) : array_like, shape = sample_shape

cdf_inverse(param)

Numerically determine the inverse CDF from a given sampling PDF. If the PDF itself is not normalized, the class will keep an internal record of the normalization and adjust the PDF values as necessary. Returns a function object which is the interpolated CDF inverse.

cdf_inverse_from_hist(P, param, old_style=False)

• for now, do on the CPU, since this is done rarely and involves fairly small arrays

• this is very wasteful, since we are casting back to the CPU for ALL our sampling points

clear()

Clear out the parameters and their settings, as well as clear the sample cache.

compute_hist(x_samples, param, weights=None, floor_level=0)

draw(rvs, *args, **kwargs)

Draw a set of random variates for parameter(s) args. Left and right limits are handed to the function. If args is None, then draw all parameters. ‘rdict’ parameter is a boolean. If true, returns a dict matched to param name rather than list. rvs must be either a list of uniform random variates to transform for sampling, or an integer number of samples to draw.

draw_simplified(n_samples, *args, **kwargs)

integrate(func, *args, **kwargs)

Integrate func, by using n sample points. Right now, all params defined must be passed to args must be provided, but this will change soon. Does NOT allow for tuples of arguments, an unused feature in mcsampler

kwargs: nmax – total allowed number of sample points, will throw a warning if this number is reached before neff. neff – Effective samples to collect before terminating. If not given, assume infinity n – Number of samples to integrate in a ‘chunk’ – default is 1000 save_integrand – Save the evaluated value of the integrand at the sample points with the sample point history_mult – Number of chunks (of size n) to use in the adaptive histogramming: only useful if there are parameters with adaptation enabled tempering_exp – Exponent to raise the weights of the 1-D marginalized histograms for adaptive sampling prior generation, by default it is 0 which will turn off adaptive sampling regardless of other settings temper_log – Adapt in min(ln L, 10^(-5))^tempering_exp tempering_adapt – Gradually evolve the tempering_exp based on previous history. floor_level – total probability of a uniform distribution, averaged with the weighted sampled distribution, to generate a new sampled distribution n_adapt – number of chunks over which to allow the pdf to adapt. Default is 1000*n, i.e. adaptation runs on ANY run longer than one chunk – it is NOT off by default. The kwarg is in CHUNKS and is scaled by n here, so the default 1000*n is exactly n_adapt=1000; the scaled variable is then compared against self.ntotal, a SAMPLE count. Do NOT pass 1000*n to reproduce the default convergence_tests - dictionary of function pointers, each accepting self._rvs and self.params as arguments. CURRENTLY ONLY USED FOR REPORTING Pinning a value: By specifying a kwarg with the same of an existing parameter, it is possible to “pin” it. The sample draws will always be that value, and the sampling prior will use a delta function at that value.

integrate_log(lnF, *args, xpy=None, **kwargs)

Integrate exp(lnF) returning lnI, by using n sample points, assuming integrand is lnF Does NOT allow for tuples of arguments, an unused feature in mcsampler

tempering is done with lnF, suitably modified.

kwargs: nmax – total allowed number of sample points, will throw a warning if this number is reached before neff. neff – Effective samples to collect before terminating. If not given, assume infinity n – Number of samples to integrate in a ‘chunk’ – default is 1000 save_integrand – Save the evaluated value of the integrand at the sample points with the sample point history_mult – Number of chunks (of size n) to use in the adaptive histogramming: only useful if there are parameters with adaptation enabled tempering_exp – Exponent to raise the weights of the 1-D marginalized histograms for adaptive sampling prior generation, by default it is 0 which will turn off adaptive sampling regardless of other settings temper_log – Adapt in min(ln L, 10^(-5))^tempering_exp tempering_adapt – Gradually evolve the tempering_exp based on previous history. floor_level – total probability of a uniform distribution, averaged with the weighted sampled distribution, to generate a new sampled distribution n_adapt – number of chunks over which to allow the pdf to adapt. Default is 1000*n, i.e. adaptation runs on ANY run longer than one chunk – it is NOT off by default. The kwarg is in CHUNKS and is scaled by n here, so the default 1000*n is exactly n_adapt=1000; the scaled variable is then compared against self.ntotal, a SAMPLE count. Do NOT pass 1000*n to reproduce the default convergence_tests - dictionary of function pointers, each accepting self._rvs and self.params as arguments. CURRENTLY ONLY USED FOR REPORTING Pinning a value: By specifying a kwarg with the same of an existing parameter, it is possible to “pin” it. The sample draws will always be that value, and the sampling prior will use a delta function at that value.

joint_p_s_is_normalized_density = True

static match_params_from_args(args, params)

Given two unordered sets of parameters, one a set of all “basic” elements (strings) possible, and one a set of elements both “basic” strings and “combined” (basic strings in tuples), determine whether the sets are equivalent if no basic element is repeated. e.g. set A ?= set B (“a”, “b”, “c”) ?= (“a”, “b”, “c”) ==> True ((“a”, “b”, “c”)) ?= (“a”, “b”, “c”) ==> True ((“a”, “b”), “d”)) ?= (“a”, “b”, “c”) ==> False # basic element ‘d’ not in set B ((“a”, “b”), “d”)) ?= (“a”, “b”, “d”, “c”) ==> False # not all elements in set B represented in set A

pdf_from_hist(x, param)

reset_sampling(param)

sampling_density(X)

Pointwise sampling density q(theta) of THIS member, evaluated at ARBITRARY points X (shape (N, ndim), columns in self.params_ordered order). Returns a host (numpy) array of length N, or None if this sampler cannot express one.

THE PORTFOLIO MEMBER CONTRACT. mcsamplerPortfolio builds the balance-heuristic mixture denominator q_mix = sum_m frac_m q_m from this method, so what it returns must be a properly NORMALIZED probability density over the sampler’s own box, in original coordinates, rather than a scale factor.

This sampler draws each parameter independently by inverse-CDF sampling, so the density is the product over parameters of pdf[p](x)/_pdf_norm[p] – exactly the per-sample product draw_simplified() reports as joint_p_s, and the two are kept in step. Zero outside the box: cdf_inv cannot produce a point there, and pdf_from_hist would otherwise return the clamped edge bin rather than 0.

PRECONDITION, not enforced here. _pdf_norm[p] is only populated when add_parameter had to build the CDF itself; when the CALLER supplies cdf_inv (what ILE does) it stays at its default 1, and the result is then normalized only if the caller’s pdf already integrates to 1. Every in-tree caller that supplies cdf_inv passes a normalized pdf (ret_uniform_samp_vector_alt is 1/(b-a), cos_samp is sin(x)/2), so this is latent. A caller doing otherwise gets a q_m off by a constant and a biased mixture, with no diagnostic. Pinned by test_portfolio_member_density.py::test_supplied_cdf_inv_with_unnormalized_pdf_is_the_known_hole.

READ-ONLY – touches no sampler state and does not affect this sampler’s own integrate().

setup(n_bins=100, **kwargs)

setup_hist()

Initializes dictionaries for all of the info that needs to be stored for the histograms, across every parameter.

setup_hist_single_param(x_min, x_max, n_bins, param)

update_sampling_prior(ln_weights, n_history, tempering_exp=1, log_scale_weights=True, floor_integrated_probability=0, external_rvs=None, **kwargs)

Default setting is for lnL. The choices adopted are fairly arbitrary. For simplicity, we are using ‘log scale weights’ (derived from lnL, not L) which lets us sample well away from the posterior at the expense of not being as strongly sampling near the peak.

NOTE: Currently deployed for mcsamplerPortfolio, NOT yet part of core code here

RIFT.integrators.mcsamplerGPU.M_samp_vector(Mmin, Mmax, x)

exception RIFT.integrators.mcsamplerGPU.NanOrInf(value)

Bases: Exception

RIFT.integrators.mcsamplerGPU.convergence_test_MostSignificantPoint(pcut, rvs, params)

RIFT.integrators.mcsamplerGPU.convergence_test_NormalSubIntegrals(ncopies, pcutNormalTest, sigmaCutRelativeErrorThreshold, rvs, params)

RIFT.integrators.mcsamplerGPU.cos_samp(x, xpy=None)

RIFT.integrators.mcsamplerGPU.cos_samp_cdf_inv_vector(p, xpy=None)

RIFT.integrators.mcsamplerGPU.cos_samp_vector(x, xpy=None)

RIFT.integrators.mcsamplerGPU.dec_samp(x, xpy=None)

RIFT.integrators.mcsamplerGPU.dec_samp_cdf_inv_vector(p, xpy=None)

RIFT.integrators.mcsamplerGPU.dec_samp_vector(x, xpy=None)

RIFT.integrators.mcsamplerGPU.delta_func_pdf(x_0, x)

RIFT.integrators.mcsamplerGPU.delta_func_samp(x_0, x)

RIFT.integrators.mcsamplerGPU.identity_convert(x)

RIFT.integrators.mcsamplerGPU.identity_convert_togpu(x)

RIFT.integrators.mcsamplerGPU.profile(fn)

RIFT.integrators.mcsamplerGPU.pseudo_dist_samp(r0, r)

RIFT.integrators.mcsamplerGPU.q_cdf_inv_vector(qmin, qmax, x, xpy=None)

RIFT.integrators.mcsamplerGPU.q_samp_vector(qmin, qmax, x)

RIFT.integrators.mcsamplerGPU.ret_uniform_samp_vector_alt(a, b)

RIFT.integrators.mcsamplerGPU.sanityCheckSamplerIntegrateUnity(sampler, *args, **kwargs)

RIFT.integrators.mcsamplerGPU.set_xpy_to_numpy()

RIFT.integrators.mcsamplerGPU.uniform_samp(a, b, x)

RIFT.integrators.mcsamplerGPU.uniform_samp_cdf_inv_vector(a, b, p)

RIFT.integrators.mcsamplerGPU.uniform_samp_dec(x, xpy=None)

Assume range known as

RIFT.integrators.mcsamplerGPU.uniform_samp_phase(x, xpy=None)

Assume range known as 0,2pi

RIFT.integrators.mcsamplerGPU.uniform_samp_psi(x, xpy=None)

Assume range known as 0,pi

RIFT.integrators.mcsamplerGPU.uniform_samp_theta(x, xpy=None)

Assume range known as

RIFT.integrators.mcsamplerGPU.uniform_samp_vector(a, b, x, xpy=None)

uniform_samp_vector_lazy:

Implement uniform sampling as multiplication by a constant. Much faster and lighter weight. We never use the cutoffs anyways, because the limits are hardcoded elsewhere.

RIFT.integrators.mcsamplerGPU.uniform_samp_withfloor_vector(rmaxQuad, rmaxFlat, pFlat, x, xpy=None)

Portfolio MC Sampler

class RIFT.integrators.mcsamplerPortfolio.MCSampler(portfolio=None, portfolio_weights=None, oracle_realizations=None, n_chunk=400000, portfolio_freeze_wt=0.05, **kwargs)

Bases: SamplerOutputMixin, object

Class to define a set of parameter names, limits, and probability densities.

add_parameter(params, pdf, **kwargs)

Add one (or more) parameters to sample dimensions. params is either a string describing the parameter, or a tuple of strings. The tuple will indicate to the sampler that these parameters must be sampled together. left_limit and right_limit are on the infinite interval by default, but can and probably should be specified. If several params are given, left_limit, and right_limit must be a set of tuples with corresponding length. Sampling PDF is required, and if not provided, the cdf inverse function will be determined numerically from the sampling PDF.

bootstrap_from_samples(samples, params=None, keep_backstop_cold=None, **kwargs)

Warm-start: forward a seed cloud to every member that supports it (e.g. the AV/VARAHA member’s live volume), EXCEPT the full-support backstop (see below). Members without bootstrap_from_samples are left cold. A warm start only shapes a member’s proposal, and the portfolio combines members with the balance-heuristic mixture density (q_mix), so a mis-seeded member costs efficiency rather than bias – BUT ONLY WHILE SOME MEMBER STILL COVERS THE SUPPORT. Column order matches self.params_ordered, which every member shares (add_parameter forwards to all members in the same order), so no per-member remapping is needed.

Only VARAHA/AV-style members (those exposing bootstrap_from_samples) are seeded directly here. GMM / adaptive-Gaussian members cannot be seeded pre-integration (their internal integrator, hence gmm_dict, does not exist until the first integrate() call), but they still warm up STRUCTURALLY during the run: the portfolio lets the cold GMM member adapt its proposal from the warm AV member’s high-likelihood draws, which is what gives the mixture a faster early n_eff than warm-AV alone. An explicit GMM pre-seed (build an initial gmm_dict from the sample cloud) is a future enhancement. q_mix keeps any cold/mis-seeded member from biasing the estimate.

THE BACKSTOP MUST STAY COLD. Seeding EVERY member destroys the coverage invariant the whole design rests on. cover_frac does not save it: it mixes a FINITE number of uniform points into the seed, and a finite point set occupies only the bins it lands in, so the seeded live volume is NOT a superset of a cold start. Measured on a tight seed in the [-5,5]^d box – fraction of the prior box covered, against 1.0 cold:

d=2: cover_frac 0.0 / 0.2 / 0.5 / 0.9 -> 0.027 / 0.634 / 0.982 / 1.000 d=4: -> 0.0015 / 0.028 / 0.104 / 0.620 d=6: -> 6.3e-05 / 0.00087 / 0.0033 / 0.0287

so at d=6 even cover_frac=0.9 leaves 97% of the box unsampled. Before this change a portfolio warm start narrowed member 0 to V=0.0033 along with everyone else. Keeping member 0 (the designated backstop, which restrict_member_range also refuses to narrow) cold makes the mixture-level guarantee TRUE instead of merely asserted.

BE PRECISE ABOUT WHAT THIS BUYS, because the measurements are not what one expects:

• With a CORRECT seed it is nearly free – n_eff 3630/3466/5695 cold-backstop vs 3876/3461/5621 seeded at d=4, and 5078/5689/5163 vs 5857/5387/5266 at d=6. Within run-to-run scatter, so it is cheap insurance. That is the case for the default.

• It is NOT what protects the DEFAULT AV+GMM portfolio. The GMM member is a Gaussian mixture with nonzero density over the whole box, so q_mix never vanishes there regardless of this setting. With a deliberately displaced seed, |lnZ bias| stayed <= 0.05 in every d=4 and d=6 run in BOTH arms. That unbounded support is the real (previously undocumented) reason production has not been biased by warm starts.

• It does NOT rescue a badly mismatched seed. In an ALL-AV portfolio (every component a hard-edged box) with a displaced seed at d=6, the cold backstop still gave lnZ bias -1.1 to -4.2 nats with n_eff 3-9, versus -1.0 to -6.8 seeded. A uniform member finds a sharp 6-D peak too rarely to carry the integral within budget. Coverage in principle is necessary, not sufficient – for a mismatched seed the mitigations are detection (L0 rescue) and not warm-starting across dissimilar points.

keep_backstop_cold: None (default) uses self.portfolio_warmstart_backstop_cold, itself True by default; pass False to restore the old seed-everything behaviour.

Returns the number of members warm-started (0 is fine; the portfolio still runs).

clear_warm_state()

Clear any warm-start seed AND the installed active grid on every member.

Setting portfolio._warm = None does NOT do this: _warm and the contracted AV grid live on the MEMBERS, not on the portfolio object. Now that a seed is actually installed on the draw path (AV._apply_warm_state), failing to clear it between points would let the next point reuse the PREVIOUS point’s contracted live volume – which can exclude the new point’s support and bias it low with no diagnostic. Called by the driver wherever it used to do sampler._warm = None, including the seed-capture failure and exception paths.

Failures PROPAGATE. A reset that quietly did not happen leaves the next point drawing from the previous point’s grid, which is the exact silent-wrong-answer this method exists to prevent – so it must not be reducible to a log line.

draw(n_samples, *args, **kwargs)

draw n_samples

Draw from portfolio. Uses portfolio weights to calculate desired outcomes. Restricts to ACTIVE members.

integrate(lnF, *args, xpy=<module 'numpy' from '/opt/hostedtoolcache/Python/3.11.18/x64/lib/python3.11/site-packages/numpy/__init__.py'>, **kwargs)

integrate_log(lnF, *args, xpy=<module 'numpy' from '/opt/hostedtoolcache/Python/3.11.18/x64/lib/python3.11/site-packages/numpy/__init__.py'>, **kwargs)

reset_adaptation()

FULL reset: member proposals AND the portfolio’s own adaptive bookkeeping.

clear_warm_state() rebuilds the MEMBERS, but the portfolio itself also learns during a run – draw allocation, per-member quality EMAs and their observation counts, the round-robin probe pointer, the iteration counter, breakpoint progression and the per-member n_ess histories. MC-error replicas that inherit those start with scheduling learned from the earlier replicas, so they are not adaptation-independent and the between-replica scatter still understates the true error – which is the entire quantity the replicas exist to measure. Restores every field to its post-setup value.

restrict_member_range(member_index, param, lo, hi)

Narrow ONE portfolio member’s sampling range for param to [lo, hi] (interval narrowing).

This is a PROPOSAL-only change: the member’s prior callables are untouched, so it keeps reporting the true global prior, and the balance-heuristic mixture density q_mix keeps the estimate unbiased with no renormalization – PROVIDED at least one member retains full support. Member 0 is that backstop by convention and may not be narrowed.

Why the backstop is not optional: measured on a truth-known ladder, a WRONG sub-box costs a STANDALONE sampler up to -1949 nats (while still reporting a healthy n_eff of 220-840, i.e. confidently wrong), but only ~1 nat inside a portfolio whose full-box member keeps q_mix covering the complement. Restriction without a backstop converts a rare pathology into a systematic one.

Call before add_parameter(); narrowing is applied there, and setup() then builds every derived quantity from the narrowed range.

setup(**kwargs)

support_diagnostics()

Finalize the support-mismatch statistics. Returns a dict; see _update_support_diagnostics for definitions and limits.

exception RIFT.integrators.mcsamplerPortfolio.NanOrInf(value)

Bases: Exception

RIFT.integrators.mcsamplerPortfolio.identity_convert(x)

RIFT.integrators.mcsamplerPortfolio.identity_convert_togpu(x)

RIFT.integrators.mcsamplerPortfolio.portfolio_default_weights(n_ess_list, wt_previous, portfolio_probability_floor=0.01, history_factor=0.5, xpy=<module 'numpy' from '/opt/hostedtoolcache/Python/3.11.18/x64/lib/python3.11/site-packages/numpy/__init__.py'>, identity_convert=<function <lambda>>, **kwargs)

RIFT.integrators.mcsamplerPortfolio.profile(fn)

RIFT.integrators.mcsamplerPortfolio.set_xpy_to_numpy()

Normalizing Flows MC Sampler

Monte Carlo Ensemble

Monte Carlo Integrator

Perform an adaptive monte carlo integral.

RIFT.integrators.MonteCarloEnsemble.identity_convert(x)

RIFT.integrators.MonteCarloEnsemble.identity_convert_togpu(x)

class RIFT.integrators.MonteCarloEnsemble.integrator(d, bounds, gmm_dict, n_comp, n=None, prior=None, user_func=None, proc_count=None, L_cutoff=None, use_lnL=False, return_lnI=False, gmm_adapt=None, gmm_epsilon=None, tempering_exp=1, temper_log=False, lnw_failure_cut=None, tempering_adapt=False, ess_target=None, ess_floor=None, gmm_adaptive=None, gmm_defensive_frac=0.05, gmm_inflate=1.0, param_names=None)

Bases: object

Class to iteratively perform an adaptive Monte Carlo integral where the integrand is a combination of one or more Gaussian curves, in one or more dimensions.

Parameters

dint

Total number of dimensions.

boundsdictionary with array bounds, with keys matching gmm_dict

Limits of integration, where each row represents [left_lim, right_lim] for its corresponding dimension.

gmm_dictdict

Dictionary where each key is a tuple of one or more dimensions that are to be modeled together. If the integrand has strong correlations between two or more dimensions, they should be grouped. Each value is by default initialized to None, and is replaced with the GMM object for its dimension(s).

n_compint or {tuple:int}

The number of Gaussian components per group of dimensions. If its type is int, this number of components is used for all dimensions. If it is a dict, it maps each key in gmm_dict to an integer number of mixture model components.

nint

Number of samples per iteration

priorfunction

Function to evaluate prior for samples

user_funcfunction

Function to run each iteration

L_cutofffloat

Likelihood cutoff for samples to store

use_lnLbool

Whether or not lnL or L will be returned by the integrand

integrate(func, min_iter=10, max_iter=20, var_thresh=0.0, max_err=10, neff=inf, nmax=None, progress=False, epoch=None, verbose=True, force_no_adapt=False, use_lnL=False, return_lnI=False, **kwargs)

RIFT.integrators.MonteCarloEnsemble.validate_gmm_dict(bounds, gmm_dict, where='gmm_dict', param_names=None)

Reject a seeded proposal that does not describe the dim-group it is installed against.

_sample() draws from model.sample() and divides the weights by model.score(), the model’s own normalized density over its OWN box. Two things therefore have to agree with the group, and neither is checked anywhere else:

• DIMENSION. A model with more axes than its group contributes a density carrying the extra axes, so lnZ comes out high by their log-volume – +1.84 nats for a 2-D (psi, phi_orb) model against the 1-element (psi,) group the ILE drivers use when phase is marginalized. Fewer axes raises IndexError on the first _sample().

• BOUNDS. A model normalized over a SMALLER box than the group’s is never wrong about any sample it draws, and never draws outside its own box – so the run silently explores a fraction of the group and lnZ is low by ln(area ratio). Replaying an –extrinsic-proposal-breadcrumb recorded without –internal-rotate-phase into a run with it does exactly this: the option doubles psi and phi_orb to (0, 4pi), and the stored seed keeps (0, 2pi). Measured: -1.39 nats, a quarter of the box reachable.

Neither failure disturbs the recovered marginals, which keep their shape while only the evidence moves, so it has to be caught here or not at all.

This is an equality check on the box, NOT an ordering check: two axes sharing a box (psi and phi_orb do) can be permuted without tripping it. Axis order is gmm_dict_from_breadcrumb._permute_group’s job.

Weighted GMM

Gaussian Mixture Models.

This implementation corresponds to frequentist (non-Bayesian) formulation of Gaussian Mixture Models.

Adapted from sklearn, added ability to specify weights.

class RIFT.integrators.weighted_gmm.WeightedGMM(n_components=1, covariance_type='diag', random_state=None, thresh=None, tol=0.001, min_covar=0.001, n_iter=100, n_init=1, params='wmc', init_params='wmc', verbose=0)

Bases: BaseEstimator

Gaussian Mixture Model with weighted samples

Representation of a Gaussian mixture model probability distribution. This class allows for easy evaluation of, sampling from, and maximum-likelihood estimation of the parameters of a GMM distribution.

Initializes parameters such that every mixture component has zero mean and identity covariance.

Read more in the User Guide.

Parameters

n_componentsint, optional

Number of mixture components. Defaults to 1.

covariance_typestring, optional

String describing the type of covariance parameters to use. Must be one of ‘spherical’, ‘tied’, ‘diag’, ‘full’. Defaults to ‘diag’.

random_state: RandomState or an int seed (None by default)

A random number generator instance

min_covarfloat, optional

Floor on the diagonal of the covariance matrix to prevent overfitting. Defaults to 1e-3.

tolfloat, optional

Convergence threshold. EM iterations will stop when average gain in log-likelihood is below this threshold. Defaults to 1e-3.

n_iterint, optional

Number of EM iterations to perform.

n_initint, optional

Number of initializations to perform. the best results is kept

paramsstring, optional

Controls which parameters are updated in the training process. Can contain any combination of ‘w’ for weights, ‘m’ for means, and ‘c’ for covars. Defaults to ‘wmc’.

init_paramsstring, optional

Controls which parameters are updated in the initialization process. Can contain any combination of ‘w’ for weights, ‘m’ for means, and ‘c’ for covars. Defaults to ‘wmc’.

verboseint, default: 0

Enable verbose output. If 1 then it always prints the current initialization and iteration step. If greater than 1 then it prints additionally the change and time needed for each step.

Attributes

weights_array, shape (n_components,)

This attribute stores the mixing weights for each mixture component.

means_array, shape (n_components, n_features)

Mean parameters for each mixture component.

covars_array

Covariance parameters for each mixture component. The shape depends on covariance_type:

(n_components, n_features)             if 'spherical',
(n_features, n_features)               if 'tied',
(n_components, n_features)             if 'diag',
(n_components, n_features, n_features) if 'full'

converged_bool

True when convergence was reached in fit(), False otherwise.

See Also

DPGMMInfinite gaussian mixture model, using the dirichlet

process, fit with a variational algorithm

VBGMMFinite gaussian mixture model fit with a variational

algorithm, better for situations where there might be too little data to get a good estimate of the covariance matrix.

Examples

>>> import numpy as np
>>> from sklearn import mixture
>>> np.random.seed(1)
>>> g = mixture.GMM(n_components=2)
>>> # Generate random observations with two modes centered on 0
>>> # and 10 to use for training.
>>> obs = np.concatenate((np.random.randn(100, 1),
...                       10 + np.random.randn(300, 1)))
>>> g.fit(obs)
GMM(covariance_type='diag', init_params='wmc', min_covar=0.001,
        n_components=2, n_init=1, n_iter=100, params='wmc',
        random_state=None, thresh=None, tol=0.001, verbose=0)
>>> np.round(g.weights_, 2)
array([ 0.75,  0.25])
>>> np.round(g.means_, 2)
array([[ 10.05],
       [  0.06]])
>>> np.round(g.covars_, 2)
array([[[ 1.02]],
       [[ 0.96]]])
>>> g.predict([[0], [2], [9], [10]])
array([1, 1, 0, 0]...)
>>> np.round(g.score([[0], [2], [9], [10]]), 2)
array([-2.19, -4.58, -1.75, -1.21])
>>> # Refit the model on new data (initial parameters remain the
>>> # same), this time with an even split between the two modes.
>>> g.fit(20 * [[0]] +  20 * [[10]])
GMM(covariance_type='diag', init_params='wmc', min_covar=0.001,
        n_components=2, n_init=1, n_iter=100, params='wmc',
        random_state=None, thresh=None, tol=0.001, verbose=0)
>>> np.round(g.weights_, 2)
array([ 0.5,  0.5])

aic(X, w=None)

Akaike information criterion for the current model fit and the proposed data

Parameters

X : array of shape(n_samples, n_dimensions)

Returns

aic: float (the lower the better)

bic(X, w=None)

Bayesian information criterion for the current model fit and the proposed data

Parameters

X : array of shape(n_samples, n_dimensions)

Returns

bic: float (the lower the better)

fit(X, w=None, y=None)

Estimate model parameters with the EM algorithm.

A initialization step is performed before entering the expectation-maximization (EM) algorithm. If you want to avoid this step, set the keyword argument init_params to the empty string ‘’ when creating the GMM object. Likewise, if you would like just to do an initialization, set n_iter=0.

Parameters

Xarray_like, shape (n, n_features)

List of n_features-dimensional data points. Each row corresponds to a single data point.

warray-like, shape = [n_samples] (optional)

Sample weights

Returns

self

fit_predict(X, w=None, y=None)

Fit and then predict labels for data.

Warning: due to the final maximization step in the EM algorithm, with low iterations the prediction may not be 100% accurate

Added in version 0.17: fit_predict method in Gaussian Mixture Model.

Parameters

X : array-like, shape = [n_samples, n_features] w : array-like, shape = [n_samples] (optional)

Sample weights

Returns

C : array, shape = (n_samples,) component memberships

predict(X)

Predict label for data.

Parameters

X : array-like, shape = [n_samples, n_features]

Returns

C : array, shape = (n_samples,) component memberships

predict_proba(X, w=None)

Predict posterior probability of data under each Gaussian in the model.

Parameters

X : array-like, shape = [n_samples, n_features]

Returns

responsibilitiesarray-like, shape = (n_samples, n_components)

Returns the probability of the sample for each Gaussian (state) in the model.

sample(n_samples=1, random_state=None)

Generate random samples from the model.

Parameters

n_samplesint, optional

Number of samples to generate. Defaults to 1.

Returns

Xarray_like, shape (n_samples, n_features)

List of samples

score(X, w=None, y=None)

Compute the log probability under the model.

Parameters

Xarray_like, shape (n_samples, n_features)

List of n_features-dimensional data points. Each row corresponds to a single data point.

Returns

logprobarray_like, shape (n_samples,)

Log probabilities of each data point in X

score_samples(X, w=None)

Return the per-sample likelihood of the data under the model.

Compute the log probability of X under the model and return the posterior distribution (responsibilities) of each mixture component for each element of X.

Parameters

X: array_like, shape (n_samples, n_features)

List of n_features-dimensional data points. Each row corresponds to a single data point.

warray-like, shape = [n_samples] (optional)

Sample weights

Returns

logprobarray_like, shape (n_samples,)

Log probabilities of each data point in X.

responsibilitiesarray_like, shape (n_samples, n_components)

Posterior probabilities of each mixture component for each observation

set_fit_request(*, w: bool | None | str = '$UNCHANGED$') → WeightedGMM

Configure whether metadata should be requested to be passed to the fit method.

Note that this method is only relevant when this estimator is used as a sub-estimator within a meta-estimator and metadata routing is enabled with enable_metadata_routing=True (see sklearn.set_config()). Please check the User Guide on how the routing mechanism works.

The options for each parameter are:

• True: metadata is requested, and passed to fit if provided. The request is ignored if metadata is not provided.

• False: metadata is not requested and the meta-estimator will not pass it to fit.

• None: metadata is not requested, and the meta-estimator will raise an error if the user provides it.

• str: metadata should be passed to the meta-estimator with this given alias instead of the original name.

The default (sklearn.utils.metadata_routing.UNCHANGED) retains the existing request. This allows you to change the request for some parameters and not others.

Added in version 1.3.

Parameters

wstr, True, False, or None, default=sklearn.utils.metadata_routing.UNCHANGED

Metadata routing for w parameter in fit.

Returns

selfobject

The updated object.

set_predict_proba_request(*, w: bool | None | str = '$UNCHANGED$') → WeightedGMM

Configure whether metadata should be requested to be passed to the predict_proba method.

Note that this method is only relevant when this estimator is used as a sub-estimator within a meta-estimator and metadata routing is enabled with enable_metadata_routing=True (see sklearn.set_config()). Please check the User Guide on how the routing mechanism works.

The options for each parameter are:

• True: metadata is requested, and passed to predict_proba if provided. The request is ignored if metadata is not provided.

• False: metadata is not requested and the meta-estimator will not pass it to predict_proba.

• None: metadata is not requested, and the meta-estimator will raise an error if the user provides it.

• str: metadata should be passed to the meta-estimator with this given alias instead of the original name.

The default (sklearn.utils.metadata_routing.UNCHANGED) retains the existing request. This allows you to change the request for some parameters and not others.

Added in version 1.3.

Parameters

wstr, True, False, or None, default=sklearn.utils.metadata_routing.UNCHANGED

Metadata routing for w parameter in predict_proba.

Returns

selfobject

The updated object.

set_score_request(*, w: bool | None | str = '$UNCHANGED$') → WeightedGMM

Configure whether metadata should be requested to be passed to the score method.

Note that this method is only relevant when this estimator is used as a sub-estimator within a meta-estimator and metadata routing is enabled with enable_metadata_routing=True (see sklearn.set_config()). Please check the User Guide on how the routing mechanism works.

The options for each parameter are:

• True: metadata is requested, and passed to score if provided. The request is ignored if metadata is not provided.

• False: metadata is not requested and the meta-estimator will not pass it to score.

• None: metadata is not requested, and the meta-estimator will raise an error if the user provides it.

• str: metadata should be passed to the meta-estimator with this given alias instead of the original name.

The default (sklearn.utils.metadata_routing.UNCHANGED) retains the existing request. This allows you to change the request for some parameters and not others.

Added in version 1.3.

Parameters

wstr, True, False, or None, default=sklearn.utils.metadata_routing.UNCHANGED

Metadata routing for w parameter in score.

Returns

selfobject

The updated object.

RIFT.integrators.weighted_gmm.distribute_covar_matrix_to_match_covariance_type(tied_cv, covariance_type, n_components)

Create all the covariance matrices from a given template

RIFT.integrators.weighted_gmm.load_model(data, random_state=None)

RIFT.integrators.weighted_gmm.log_multivariate_normal_density(X, means, covars, covariance_type='diag')

Compute the log probability under a multivariate Gaussian distribution.

Parameters

Xarray_like, shape (n_samples, n_features)

List of n_features-dimensional data points. Each row corresponds to a single data point.

meansarray_like, shape (n_components, n_features)

List of n_features-dimensional mean vectors for n_components Gaussians. Each row corresponds to a single mean vector.

covarsarray_like

List of n_components covariance parameters for each Gaussian. The shape depends on covariance_type:

(n_components, n_features) if ‘spherical’, (n_features, n_features) if ‘tied’, (n_components, n_features) if ‘diag’, (n_components, n_features, n_features) if ‘full’

covariance_typestring

Type of the covariance parameters. Must be one of ‘spherical’, ‘tied’, ‘diag’, ‘full’. Defaults to ‘diag’.

Returns

lprarray_like, shape (n_samples, n_components)

Array containing the log probabilities of each data point in X under each of the n_components multivariate Gaussian distributions.

RIFT.integrators.weighted_gmm.sample_gaussian(mean, covar, covariance_type='diag', n_samples=1, random_state=None)

Generate random samples from a Gaussian distribution.

Parameters

meanarray_like, shape (n_features,)

Mean of the distribution.

covararray_like, optional

Covariance of the distribution. The shape depends on covariance_type:

scalar if ‘spherical’, (n_features) if ‘diag’, (n_features, n_features) if ‘tied’, or ‘full’

covariance_typestring, optional

Type of the covariance parameters. Must be one of ‘spherical’, ‘tied’, ‘diag’, ‘full’. Defaults to ‘diag’.

n_samplesint, optional

Number of samples to generate. Defaults to 1.

Returns

Xarray, shape (n_features, n_samples)

Randomly generated sample

Statistical Utilities

RIFT.integrators.statutils.block_scatter_sigma(lnZ_blocks, n_blocks)

sigma(lnZ) from the between-chunk scatter of per-chunk mean estimates, via a delete-one jackknife of the n-weighted pooled mean. Each chunk of an adaptive run used a different proposal, so this sees the nonstationarity the pooled within-run variance averages away. (It still cannot see modes that EVERY chunk missed – only independent replicas can.) Returns float sigma or None if fewer than 2 usable chunks.

RIFT.integrators.statutils.bootstrap_lnZ_quantiles(log_wt, n_total=None, n_boot=200, quantiles=(0.05, 0.5, 0.95), rng_seed=None, n_drawn=None)

Bootstrap quantiles of lnZ_hat = ln( sum_i w_i / n_total ) by resampling the stored LOG weights with replacement. When the relative error is O(1) the delta-method +-sigma interval on lnZ is meaningless (the distribution is strongly skewed); these quantiles are an honest same-sample interval. They remain blind to tail mass never sampled – pair with pareto_khat_from_log. n_total: divisor if the stored weights are a pruned subset of a larger run (the pruned-away weights contribute negligibly to the sum). n_drawn: if the stored weights are the nonzero terms of an average over n_drawn draws (the rest rejected as zeros), resample all n_drawn draws: each replicate keeps k ~ Binomial(n_drawn, n/n_drawn) of the stored weights, so the survivor count varies too. The divisor stays n_total, which is the normalization of the stored weights. Returns a numpy array of lnZ quantiles, or None if too few weights.

RIFT.integrators.statutils.cumvar(arr, mean=None, var=None, n=0)

Numerically stable running sample variance measure. If mean and var are supplied, they will be used as the history values. See

http://www.johndcook.com/standard_deviation.html

for algorithm details.

RIFT.integrators.statutils.ess_from_log_weights(log_wt)

Kish effective sample size (sum w)^2 / sum w^2 from LOG weights.

RIFT.integrators.statutils.finalize(existingAggregate)

RIFT.integrators.statutils.finalize_log(existingAggregate, xpy=<module 'numpy' from '/opt/hostedtoolcache/Python/3.11.18/x64/lib/python3.11/site-packages/numpy/__init__.py'>)

fallback mode: if special is fallback, then the aggregate is internally done with the normal numpy calculation!

RIFT.integrators.statutils.init_log(newLogValues_orig, special=<module 'scipy.special' from '/opt/hostedtoolcache/Python/3.11.18/x64/lib/python3.11/site-packages/scipy/special/__init__.py'>, xpy=<module 'numpy' from '/opt/hostedtoolcache/Python/3.11.18/x64/lib/python3.11/site-packages/numpy/__init__.py'>)

fallback mode: if special is fallback, then the aggregate is internally done with the normal numpy calculation!

RIFT.integrators.statutils.int_var(samples)

RIFT.integrators.statutils.pareto_khat_from_log(log_wt, tail_frac=0.2, min_tail=20)

Generalized-Pareto tail index k of the importance-weight distribution, fit to the largest weights (Zhang & Stephens 2009 posterior-mean estimator, as used by Pareto-smoothed importance sampling). Input is LOG weights on any scale (k is invariant under overall rescaling). Interpretation: k < 0.5 : weight variance finite, the naive error estimate is meaningful; 0.5-0.7 : variance marginal, treat the naive sigma as optimistic; k > 0.7 : the weight tail is unresolved – the naive sigma is a lower bound

and the integral itself may be dominated by unseen tail mass.

Returns float k, or None if there are too few finite weights to fit.

RIFT.integrators.statutils.update(existingAggregate, newValues, xpy=<module 'numpy' from '/opt/hostedtoolcache/Python/3.11.18/x64/lib/python3.11/site-packages/numpy/__init__.py'>)

RIFT.integrators.statutils.update_log(existingLogAggregate, newLogValues_orig, special=<module 'scipy.special' from '/opt/hostedtoolcache/Python/3.11.18/x64/lib/python3.11/site-packages/scipy/special/__init__.py'>, xpy=<module 'numpy' from '/opt/hostedtoolcache/Python/3.11.18/x64/lib/python3.11/site-packages/numpy/__init__.py'>)

logsumexp : warning it is implemented but has a different function name, need to wrap it carefully and detect which is used

RIFT.integrators.statutils.welford(x_array, mean=None, var=None, n=0)

https://www.embeddedrelated.com/showarticle/785.php see also https://brenocon.com/blog/2008/11/calculating-running-variance-in-python-and-c/

No reallocations, unlike ‘cumvar’ below!

Gaussian Mixture Model

Gaussian Mixture Model

Fit a Gaussian Mixture Model (GMM) to data and draw samples from it. Uses the Expectation-Maximization algorithm.

Weighted data GMM formulae: different from framework in eg. https://arxiv.org/pdf/1509.01509.pdf

RIFT.integrators.gaussian_mixture_model.add_defensive_component(model, defensive_frac=0.05, width_norm=1.0)

Append a broad, box-covering “defensive” component to a fitted mixture so the proposal has heavy enough tails for importance sampling.

This is the single most important fix for the SNR~82 extrinsic posterior: a mixture fit to the (tight) elite cloud UNDER-COVERS the broad, degenerate directions (distance-inclination), so the importance weight L*p/q blows up on the rare draw that lands in a poorly-covered high-likelihood pocket and the effective sample size collapses to ~1. A defensive component (Hesterberg 1995) with weight defensive_frac, wide in the model’s normalized [-1,1] frame, bounds the weights: q >= defensive_frac * q_broad everywhere, so no single sample can dominate. The AV sampler gets the same guarantee from its cover-fraction floor; the fitted GMM had none.

width_norm is the std of the defensive Gaussian in normalized coords (1.0 ~ covers the whole [-1,1] box; truncated to the box it is near-uniform).

class RIFT.integrators.gaussian_mixture_model.estimator(k, max_iters=100, tempering_coeff=1e-08, adapt=None, xpy=None)

Bases: object

Base estimator class for GMM

Parameters

kint

Number of Gaussian components

max_itersint

Maximum number of Expectation-Maximization iterations

fit(sample_array, log_sample_weights)

Fit the model to data

print_params()

Prints the model’s parameters in an easily-readable format

RIFT.integrators.gaussian_mixture_model.fit_gmm_adaptive(sample_array, bounds, log_sample_weights=None, k_max=8, k_min=1, k_candidates=None, epsilon=None, tempering_coeff=1e-08, prune_weight_floor=0.001, defensive_frac=0.05, inflate=1.0)

Fit a GMM whose COMPONENT COUNT is chosen from the data by BIC, then prune near-zero-weight components. Data-driven replacement for a hard-coded per-group component count.

Rationale (measured on the S250114ax extrinsic posterior, SNR~82):

• A fixed SMALL k (e.g. the correlate-all default of 2) cannot wrap a curved distance-inclination degeneracy arc: the elite fit is one broad Gaussian over the ridge, the proposal never locks onto the peak, and the honest effective sample size stays ~1.

• A fixed LARGE k is both statistically fragile (a spurious sharp component collapses onto ~1 elite sample and dominates the importance weights) and computationally costly (score() does an O(k) per-component box normalization on the CPU).

BIC threads between the two: fit k over a ladder, penalize free parameters by ln(N_eff), keep the best, and drop dead components. It allocates more components only where the (importance-weighted) cloud is genuinely non-Gaussian and stays at k=1 for a single blob.

SAFETY FLOOR: k is chosen in [k_min, k_max], and pruning never drops below k_min. Pass k_min = the stress-tested hard-coded per-group count so opting into adaptive can only ADD components where the data earns them, never fewer than the layout that was validated for the primary ILE use case (e.g. a broad multi-modal sky keeps its default components even if the INITIAL elite cloud – fit before the proposal has explored every mode – looks single-peaked).

Parameters

sample_array : (N, d) array in ORIGINAL coordinates. bounds : (d, 2) array of [llim, rlim] per dimension (as gmm expects). log_sample_weights : (N,) importance/elite log-weights (default: equal). k_max : cap on the number of components. k_min : floor on the number of components (default 1); the stress-

tested hard-coded count when used as a refinement layer.

k_candidates : explicit ladder (overrides k_max/k_min-derived ladder). prune_weight_floor : components below this mixture weight are removed (but

never below k_min).

Returns a fitted gmm.

class RIFT.integrators.gaussian_mixture_model.gmm(k, bounds, max_iters=1000, epsilon=None, tempering_coeff=1e-08, memory_factor=3.0, xpy=None)

Bases: object

More sophisticated implementation built on top of estimator class

fit(sample_array, log_sample_weights=None)

Fit the model to data

num_free_params()

Number of free parameters of this k-component d-dim mixture: k means (k*d) + k covariances (k*d*(d+1)/2) + (k-1) mixture weights.

print_params()

Prints the model’s parameters in an easily-readable format

prune_components(weight_floor=0.001, min_keep=1)

Drop mixture components whose weight falls below weight_floor and renormalize. Over-allocated components collapse to ~zero weight under EM; removing them (a) prevents a spurious sharp component from dominating the importance weights and (b) cuts score() cost, which is O(k) in the per-component box normalization (_box_mass, a few ms per component). Keeps at least max(1, min_keep) components (the highest-weight ones) – pass min_keep to preserve a safety floor. No-op if nothing is below the floor.

sample(n, use_bounds=True)

Draw samples from the current model.

Note the model’s means/covariances are stored in normalized coordinates (the [-1, 1] image of self.bounds under self._normalize). Samples are therefore drawn in normalized coordinates and then unnormalized back to the original coordinate frame before being returned, matching the pre-port behavior expected by MonteCarloEnsemble._sample().

score(sample_array, assume_normalized=True)

Score samples under the current model.

update(sample_array, log_sample_weights=None)

Updates the model with new data without doing a full retraining.

RIFT.integrators.gaussian_mixture_model.gpu_logpdf(x, mean, cov, xpy)

GPU-compatible multivariate normal log-pdf. x: (n, d) array mean: (d,) array cov: (d, d) array

Uses Cholesky + a generic linear solve (xpy.linalg.solve) so the same code path works for both numpy and cupy. Note: solve_triangular is NOT in numpy.linalg or cupy.linalg (only in scipy.linalg / cupyx.scipy.linalg), so we deliberately use the generic solver here.

RIFT.integrators.gaussian_mixture_model.identity_convert(x)

RIFT.integrators.gaussian_mixture_model.identity_convert_togpu(x)

RIFT.integrators.gaussian_mixture_model.mvnun(lower, upper, mean, cov, maxpts=None, abseps=1e-05, releps=1e-05)

Modern fallback for scipy.stats.mvn.mvnun using multivariate_normal.cdf. Requires SciPy 1.10.0+ for the ‘lower_limit’ parameter.

Direct Quadrature

RIFT.integrators.direct_quadrature.marginal_function(f, x_rest, args_and_ranges, g=None, epsrel=1e-05, method='default')

f: function, assumed vectorized x_rest: variable with correct range for f (i.e., f(x_rest) is valid). x_rest assumed array of options args_and_ranges: list of [ indx, [low,high],type] showing indexes of parameters and ranges. Usually do one, two, or three dimensions. ‘type’ =[‘uniform’,npts] with npts=100 is default g : optional additional function used to specify priors, if not uniform

Returns: int f(x1,x2,….xn) g(x1) dx1 if one variable provided, etc

Multivariate Truncated Normal

Multivariate Truncnorm

Draw samples from a truncated multivariate Gaussian distribution.

RIFT.integrators.multivariate_truncnorm.sample(mean, cov, bounds, n)

Generate samples

Parameters

meannp.ndarray

Mean of the distribution

covnp.ndarray

Covariance matrix of the distribution

boundsnp.ndarray

Bounds for samples

nint

Number of samples to draw

Returns

np.ndarray

Array of samples

Unreliable Oracle

unreliable_oracle – proposal ‘oracles’ for the RIFT portfolio integrator.

An oracle only PROPOSES candidate points; the portfolio evaluates the true likelihood there and folds them into the training data used to adapt the other integrators. Because oracles never contribute to the integral estimate directly, an inaccurate oracle cannot bias the result – it can only waste a few likelihood evaluations. That makes them a safe channel for injecting cheap, approximate posterior knowledge (a Fisher matrix, a hill-climbed hotspot, or a previous run’s posterior samples).

class RIFT.integrators.unreliable_oracle.ClimbingOracle(**kwargs)

Bases: MCSamplerGeneric

add_parameter(params, pdf, **kwargs)

Add one (or more) parameters to sample dimensions. params is either a string describing the parameter, or a tuple of strings. The tuple will indicate to the sampler that these parameters must be sampled together. left_limit and right_limit are on the infinite interval by default, but can and probably should be specified. If several params are given, left_limit, and right_limit must be a set of tuples with corresponding length. Sampling PDF is required, and if not provided, the cdf inverse function will be determined numerically from the sampling PDF.

draw_simplified(n_samples, *args, **kwargs)

draw n_samples

Should not be implemented in generic class

setup(lnL=None, n_climbers=10, **kwargs)

update_sampling_prior(ln_weights, n_history, lnw_cut=-10, external_rvs=None, verbose=False, **kwargs)

Default setting is for lnL. The choices adopted are fairly arbitrary. For simplicity, we are using ‘log scale weights’ (derived from lnL, not L) which lets us sample well away from the posterior at the expense of not being as strongly sampling near the peak.

NOTE: Currently deployed for mcsamplerPortfolio, NOT yet part of core code here

class RIFT.integrators.unreliable_oracle.FisherGaussianOracle(**kwargs)

Bases: MCSamplerGeneric

add_parameter(params, pdf, **kwargs)

Add one (or more) parameters to sample dimensions. params is either a string describing the parameter, or a tuple of strings. The tuple will indicate to the sampler that these parameters must be sampled together. left_limit and right_limit are on the infinite interval by default, but can and probably should be specified. If several params are given, left_limit, and right_limit must be a set of tuples with corresponding length. Sampling PDF is required, and if not provided, the cdf inverse function will be determined numerically from the sampling PDF.

draw_simplified(n_samples, *args, **kwargs)

Draw n_samples from the box-truncated Gaussian. Returns (p_s, p_prior, rv) with rv shape (n_samples, ndim); p_s is the Gaussian proposal density (so the draws can be used with correct importance weights if desired), p_prior is left None (the portfolio supplies it).

setup(mean=None, cov=None, fisher=None, adapt=False, cov_inflate=1.0, params=None, **kwargs)

mean/cov describe the proposal; alternatively pass a Fisher matrix (cov = fisher^{-1}). params names the columns of mean/cov if they are not already in this oracle’s params_ordered order. adapt=True refreshes (mean,cov) from weighted history on each update_sampling_prior call.

update_sampling_prior(ln_weights, n_history, lnw_cut=-10, external_rvs=None, verbose=False, **kwargs)

Optionally refit (mean, cov) from the recent weighted history. If adapt is False (default), the oracle keeps the supplied Fisher/Gaussian and this is a no-op – the robust behaviour when a trustworthy Fisher is available. With adapt=True it becomes a self-refining Gaussian proposal.

class RIFT.integrators.unreliable_oracle.LadderClimbingOracle(**kwargs)

Bases: MCSamplerGeneric

add_parameter(params, pdf=None, **kwargs)

Add one (or more) parameters to sample dimensions. params is either a string describing the parameter, or a tuple of strings. The tuple will indicate to the sampler that these parameters must be sampled together. left_limit and right_limit are on the infinite interval by default, but can and probably should be specified. If several params are given, left_limit, and right_limit must be a set of tuples with corresponding length. Sampling PDF is required, and if not provided, the cdf inverse function will be determined numerically from the sampling PDF.

climb(x0, scale, label='free', pinned=None, n_generations=None)

Climb from x0 with per-axis step scale; pinned = {index: value}.

draw_simplified(n_samples, *args, delta_lnL=10.0, **kwargs)

Draw from the final Gaussians of climbs whose best is within delta_lnL.

near_peak(delta_lnL=10.0)

History points within delta_lnL of the overall best: (X, lnL, labels).

place(n_samples, delta_lnL=10.0, inflate=1.5, pinned_sd=None)

Covariance placement: equal draws per good climb from N(best_x, inflate^2 cov).

A rung’s pinned coordinate has zero variance in its CMA covariance; give it sd pinned_sd (e.g. half the rung spacing) so seeds fill between rungs.

run(x0, scale, between_climbs=None)

Free climb(s), then the ladder outward from the best free optimum.

x0 is one start (ndim,) or several (n_start, ndim); each gets a free climb. between_climbs(x_start) is called before each climb (e.g. to re-adapt an auxiliary proposal at the start point).

setup(lnL=None, periodic=(), ladder_param=None, ladder_values=(), popsize=10, n_generations=25, n_generations_ladder=None, sigma0=0.3, seed=None, verbose=False, callback=None, **kwargs)

lnL : callable(*cols) -> array. ladder_param : name pinned on rungs. sigma0 : initial step, in units of the scale passed to update (or the box width when no scale is known). callback(label, gen, best_lnL) is called once per generation (for logging / pre-emption).

update_sampling_prior(ln_weights, n_history, external_rvs=None, lnw_cut=-10, **kwargs)

Start from the weighted history: best point, weighted spread as scale.

class RIFT.integrators.unreliable_oracle.PuffballOracle(**kwargs)

Bases: MCSamplerGeneric

add_parameter(params, pdf, **kwargs)

Add one (or more) parameters to sample dimensions. params is either a string describing the parameter, or a tuple of strings. The tuple will indicate to the sampler that these parameters must be sampled together. left_limit and right_limit are on the infinite interval by default, but can and probably should be specified. If several params are given, left_limit, and right_limit must be a set of tuples with corresponding length. Sampling PDF is required, and if not provided, the cdf inverse function will be determined numerically from the sampling PDF.

draw_simplified(n_samples, *args, **kwargs)

draw n_samples

Should not be implemented in generic class

setup(reference_mean=None, reference_cov=None, **kwargs)

update_sampling_prior(ln_weights, n_history, lnw_cut=-10, external_rvs=None, verbose=False, **kwargs)

Default setting is for lnL. The choices adopted are fairly arbitrary. For simplicity, we are using ‘log scale weights’ (derived from lnL, not L) which lets us sample well away from the posterior at the expense of not being as strongly sampling near the peak.

NOTE: Currently deployed for mcsamplerPortfolio, NOT yet part of core code here

class RIFT.integrators.unreliable_oracle.ResamplingOracle(**kwargs)

Bases: MCSamplerGeneric

add_parameter(params, pdf, **kwargs)

Add one (or more) parameters to sample dimensions. params is either a string describing the parameter, or a tuple of strings. The tuple will indicate to the sampler that these parameters must be sampled together. left_limit and right_limit are on the infinite interval by default, but can and probably should be specified. If several params are given, left_limit, and right_limit must be a set of tuples with corresponding length. Sampling PDF is required, and if not provided, the cdf inverse function will be determined numerically from the sampling PDF.

draw_simplified(n_samples, *args, **kwargs)

draw n_samples

Should not be implemented in generic class

setup(reference_samples=None, reference_params=None, **kwargs)

update_sampling_prior(*args, **kwargs)

Default setting is for lnL. The choices adopted are fairly arbitrary. For simplicity, we are using ‘log scale weights’ (derived from lnL, not L) which lets us sample well away from the posterior at the expense of not being as strongly sampling near the peak.

NOTE: Currently deployed for mcsamplerPortfolio, NOT yet part of core code here