All inputs are dimensional. Returned arrays are read-only.
import numpy as np
from pyimr import NeoHookeanKelvinVoigt, SimulationConfig, simulate
t = np.linspace(0.0, 120e-6, 300)
config = SimulationConfig(
R0=225e-6,
Req=37.5e-6,
material=NeoHookeanKelvinVoigt(shear_modulus_pa=2500.0, viscosity_pa_s=0.1),
)
result = simulate(t, config)result carries time_s, radius_ratio, radius_m, wall_velocity_m_s,
internal_pressure_pa, stress_integral_pa and stats. Thermal
configurations also return gas and liquid temperature fields and, when enabled,
vapor mass fraction. Materials with memory expose their internal stress state.
Inactive fields are None.
Unsupported option values and inconsistent thermal/mass-transfer combinations
fail during configuration. Integration failures and material-domain violations
raise SimulationError.
For repeated solves with one configuration, preparation hoists constant work -- state layout, grids, finite-difference operators, constitutive quadrature, Jacobian sparsity:
problem = prepare(config)
first = problem.solve(t)
second = problem.solve(t) # immutable setup is reused; solve state is freshfrom pyimr import InitialState, PhysicalParameters, SampledForcing
config = SimulationConfig(
R0=225e-6, Req=37.5e-6,
material=NeoHookeanKelvinVoigt(2500.0, 0.1),
physics=PhysicalParameters(polytropic_exponent=1.47),
initial=InitialState(wall_velocity_m_s=2.0, internal_pressure_pa=1.5e5),
)Sampled forcing values are pressure perturbations relative to the far-field baseline. A shape-preserving cubic interpolant is used between samples, and the perturbation is zero outside their time span:
config = SimulationConfig(
R0=225e-6, Req=37.5e-6,
material=NeoHookeanKelvinVoigt(2500.0, 0.1),
sampled_forcing=SampledForcing(
time_s=tuple(measured_time_s),
pressure_pa=tuple(measured_pressure_perturbation_pa),
),
)The tangent-linear solver differentiates the production RHS rather than a reduced surrogate. It covers every operator, every typed material, thermal and mass-transfer states, distributed nonlinear memory, forcing, geometry, initial conditions, and continuous physical parameters.
sensitivity = problem.solve_with_sensitivities(
t,
("R0", "material.shear_modulus_pa", "material.viscosity_pa_s", "physics.polytropic_exponent"),
)
sensitivity.simulation
sensitivity.radius_m # shape: (time, parameter)
sensitivity.state # complete internal-state derivativesParameter paths follow the frozen configuration objects, and derivatives are with respect to dimensional values. All requested directions share one augmented integration; prepared parameter scaling keeps error control dimensionless.
The mechanical path, including distributed memory, uses a cached compiled directional kernel. After its one-time compilation, six simultaneous NHKV gradients take about 1.9 times one prepared forward solve.
Tangent accuracy is not uniform across configurations. The mechanical path is limited by the finite-difference check it is measured against; the coupled thermal paths are limited by how accurately the augmented state/tangent system is integrated, which is a real bound on anything built on those gradients rather than a defect in the tangent equations. Measured error by configuration is in accuracy.md.
Prepared inference uses normalized bounded coordinates, dimensional Gaussian radius likelihoods, analytic sensitivity Jacobians, deterministic Latin-hypercube starts, and optional process-parallel batch evaluation:
from pyimr.inference import InferenceParameter, RadiusObservation, prepare_inference
inference = prepare_inference(
config,
RadiusObservation(measured_time_s, measured_radius_m, standard_deviation_m=2e-6),
(
InferenceParameter("material.shear_modulus_pa", 500.0, 5000.0, "log"),
InferenceParameter("material.viscosity_pa_s", 0.01, 1.0, "log"),
),
)
batch = inference.evaluate_batch(unit_parameter_matrix, workers=4)
fit = inference.fit_multistart(64, seed=7, workers=4)
fit.endpoints # every successful and unsuccessful endpoint is retained
fit.bestUnit parameter vectors always lie in [0, 1]; the configured linear or
logarithmic transform maps them to physical bounds. Multistart results never
discard alternative basins.
Observables do not converge at the same rate, so the right tolerance depends on
which one you fit -- internal pressure is roughly two orders behind radius at
the same setting. rtol=1e-6, atol=1e-8 is ample for likelihood evaluation
against experimental radius data; keep 1e-10, 1e-12 for sensitivities and
1e-9 or tighter for validation. The per-observable table is in
accuracy.md.
Tolerance does not bound how long one solve can take. max_steps (default
1_000_000) turns a trajectory that will not finish into a SimulationError at
a point of your choosing, which is what makes a grid sweep affordable. Three
guards refuse a trajectory the model has left rather than returning it:
max_radius_ratio (default 50) raises when the bubble runs away,
min_radius_ratio (off by default) when it collapses past where the material
is trustworthy, and max_wall_mach (off by default) when a rebound approaches
the sonic line. The last two are different guards for a reason;
boundaries.md has the measurement. When a solve fails and the
message does not say why, pyimr.diagnose does (below).
Nt and tolerance requirements depend on record length, material stiffness and
which observable is fitted, so a setting adequate for one collapse can be badly
wrong over five. pyimr.resolution measures on your own problem:
from pyimr.resolution import choose_resolution
setting = choose_resolution(config, times, target=1e-3, field="radius_ratio")
# Resolution(thermal='fd', Nt=5, rtol=1e-06, atol=1e-08,
# achieved=3.9e-05, seconds=0.0037)
config = setting.apply(config)It builds a reference and checks it is converged, searches both spectral and
fd for the cheapest grid meeting target, then loosens tolerance as far as
that grid allows. Roughly 18 solves -- worth paying before a sampling or
sensitivity campaign, not worth paying for a single run.
target is relative to each field's own peak magnitude. That matters because
observables do not converge together: at identical settings, relative error was
3.4e-07 for radius and 2.8e-05 for internal pressure. It raises rather than
guessing if the reference is not converged or the target is out of reach, since
a number built on either is indistinguishable from a real answer.
Constitutive models for soft matter nest -- NHKV is qKV at zero strain
stiffening and SLS at zero relaxation time -- so comparing best fits always
favours the flexible ones. pyimr.selection scores them by evidence instead:
from pyimr.selection import STANDARD_MODELS, compare, log_evidence, redundancy_over_grid, solve_grid
evidences = {}
for candidate in STANDARD_MODELS.values():
points, normalized, radii, stresses = solve_grid(candidate, solve, count=12)
redundancies = redundancy_over_grid(candidate, STANDARD_MODELS, points, stresses, solve)
evidences[candidate.name], _ = log_evidence(
radii, normalized, redundancies, observed, deviations, dimension=candidate.dimension
)
posterior = compare(evidences)pyimr.noise supplies the strain-rate weighting and the marginalized noise
scale; pyimr.prior the redundancy and Occam penalties. Use ONE grid count
for every model compared -- mixed resolutions let grid luck decide which lands
nearest the truth.
The grid is a screening tool, and quoting its best-fit point as a converged
estimate is a mistake. On the 15C record, refining a qSLS grid from 12 to 20 per
axis moved the log evidence by 17.3 and moved the best-fit point from
g = 658, alpha = 1.52 to g = 144, alpha = 8.86. The "best" cell is simply
whichever node lands nearest the ridge.
Refining further does not rescue it. The posterior standard deviations there are
0.0018, 0.0025, 0.0055, 0.078 in unit coordinates against a grid spacing of
0.053 at 20 per axis, so three of the four directions are narrower than the
spacing -- the posterior occupies about 1e-7 of the prior volume, and resolving
it on a uniform grid would need of order 1e10 points.
The ranking survives this: 17.3 of discretisation error sits well inside a qSLS-to-SLS gap of 101.6 log units. The parameter estimates and the evidence value do not.
For parameter estimates, use pyimr.pymc_op.sample_smc, which concentrates
where the mass is and returns an evidence with an error bar. NUTS
(sample_posterior) also works and is not pathological here -- zero divergences,
tree depth well under its cap -- but costs about 87 ms per gradient and 31
leapfrog steps per iteration, so roughly 1.5 hours per chain against SMC's
half hour (#216).
Always report the best chi-squared per sample alongside the posterior. Model
selection only means something where some candidate actually fits; otherwise the
winner is the least-bad member of an inadequate set, and the posteriors look
just as confident. Worked studies are in examples/.
solve_grid is a Cartesian product at one count on every axis, so it costs
count**dimension and runs out at four or five parameters. At count = 12, one
six-axis candidate is 2,985,984 solves against 139,272 for the whole of
STANDARD_MODELS.
Candidates past that limit live in EXTENDED_MODELS and are scored by expansion
about a fit instead of by quadrature over a grid:
from pyimr.selection import EXTENDED_MODELS, candidate_log_evidence, fit_candidate
candidate = EXTENDED_MODELS["qSLS2"]
fit = fit_candidate(candidate, solve, observed, deviation)
evidence = candidate_log_evidence(candidate, solve, observed, deviation, fit.unit)fit_candidate is multistart least squares in the same unit coordinates, and it is what
locates the point the expansion needs: PreparedInference.fit_multistart works from material
attribute paths, and a candidate's axes need not be material fields at all --- tau_ratio is
a ratio of two of them. It returns the distinct modes it found, not only the winner, because
the expansion is about one mode and the evidence is over all of them; sum with logsumexp
over fit.modes. Check fit.failure_fraction --- near one, the search spent its budget
somewhere the material will not integrate and the result is not a fit --- and note that the
differenced Jacobian costs p + 1 solves an iteration, so an under-converged search reports
a chi_squared worse than the truth's rather than announcing itself.
It costs 1 + 2*dimension solves. The Jacobian is differenced in the unit
coordinates the prior is uniform on, which is both where the Occam factor has to
be measured and what lets it handle candidates whose axes are not material fields
-- qSLS2 has tau_ratio, a ratio of two of them, and oldroydb likewise.
DYNAMICS_MODELS names the fourteen operators, as (dynamics, liquid_eos) pairs. Six
dynamics --- Rayleigh--Plesset, the Keller--Miksis pressure form, and the enthalpy forms
keller-enthalpy, herring, gilmore and lezzi-prosperetti-2 --- of which the last four
take one of three equations of state: tait, mie-gruneisen or nasg. keller-enthalpy and
herring are the lambda = 0 and lambda = 1 members of the Prosperetti--Lezzi (1986)
first-order family, not separate theories; gilmore is Kirkwood--Bethe, with the local wall
sound speed; lezzi-prosperetti-2 is their 1987 second-order equation (8.7) at the authors'
recommended (lambda, theta) = (0.5, 0), and is the only one implicit in the acceleration.
It requires a steady far field, because the second-order far-field terms are dropped. They are not
CandidateModels --- a candidate is a material, and this is the operator it is pushed
through --- so they compare by holding the candidate fixed and varying the solve callback:
from pyimr.selection import DYNAMICS_MODELS
for dynamics, liquid_eos in DYNAMICS_MODELS:
... # build a solve callback at this operator, then fit and score as aboveThe parameter space is identical across the set, so the Occam terms cancel and the difference
in log evidence is a Bayes factor between operators. Every candidate in this package assumes
dynamics="keller-miksis"; on the records analysed in the companion analysis repository, two other operators beat it.
#294 records what that comparison does and does not establish --- in
particular, it must be run in identified coordinates, or the ranking follows the prior box
rather than the data, and the operator is the most absorbed of the model axes, so most of an
operator change can be mimicked by refitting the material.
laplace_log_evidence takes cap_at_prior. Leave it off and the plain expansion
rewards a model for a parameter its design does not probe: an unidentified
direction sends an eigenvalue of J^T J to zero, so -log det / 2 goes to
+infinity. Measured on synthetic one-mode data, the two-mode model beat the
one-mode by 29.7 nats exactly where its second arm did nothing, and lost by 2.6
only where the arm did real work -- backwards, and not by a little.
cap_at_prior=True bounds each eigendirection's contribution at one, which is
what a uniform prior on a bounded cube implies: a posterior cannot be wider than
the prior it came from. It is a strict generalisation, agreeing with the plain
form to round-off wherever every direction is already sharper than the prior, so
it changes an answer only where the plain form was not entitled to one.
candidate_log_evidence turns it on, because it compares across dimensions.
Before a record exists, the question is which experiment to run. Four modules answer it. Every number they return is in nats, so the answers can be compared and added.
pyimr.design scores a design that has not been run. A DesignInference is a
PreparedInference whose observations are placeholders: the same
configuration, time grid, noise level and parameter box, with no data.
expected_information_gain averages the Laplace/Fisher gain over prior draws
and returns its Monte Carlo error bar alongside:
from pyimr.design import design_inference, expected_information_gain
design = design_inference(config, times, standard_deviation_m=2e-6, parameters=parameters)
score = expected_information_gain(design, draws=128, workers=4)
score.expected_information_gain, score.standard_error # natsoptimize_design in pyimr.optimize searches a continuous design space for
the largest gain. You supply build_inference(design) -> DesignInference, so a
design is whatever you say it is -- pulse amplitude, window, radius, or a
mixture -- and the error bar is passed to the surrogate as observation noise so
that a point which scored well by luck is not chased:
from pyimr.optimize import optimize_design
search = optimize_design(build_inference, [(0.0, 8e4)], draws=64, evaluations=24)
search.best_point, search.best_valueA single best design is the answer to the wrong question: an experimenter runs
n_1 bubbles at one setting and n_2 at another, and a search over one point
is not convex, so nothing in its result says whether it stopped at the optimum.
pyimr.measure optimises over a probability measure on a candidate set instead.
The averaged information matrix is linear in the weights, the criterion is
concave, and the Kiefer-Wolfowitz equivalence theorem turns the first-order
condition into a proof of global optimality: gap is the largest directional
derivative toward any candidate, and it is zero at the optimum.
import numpy as np
from pyimr.design import design_information
from pyimr.measure import apportion, optimal_measure
matrices = []
for R0 in (100e-6, 200e-6, 400e-6):
design = design_inference(config_at(R0), times, 2e-6, parameters)
information, requested, failed = design_information(design, draws=64)
matrices.append(information.mean(axis=0)) # one prior-averaged J^T J per candidate
measure = optimal_measure(np.array(matrices))
measure.weights, measure.support, measure.gap # weights [0, 0.54, 0.46], gap ~1e-10
batch = apportion(measure.weights, 12, np.array(matrices))
batch.counts, batch.efficiency # [0, 6, 6], D-efficiency 0.998apportion is the efficient rounding of Pukelsheim and Rieder. Passing
matrices makes it measure the D-efficiency of the integer batch against the
measure rather than assume the rounding was harmless, which for a dozen runs it
need not be.
A plain measure is not always the right object. An information criterion is
happiest concentrating, and a batch on one setting cannot detect that every
model is wrong, so constrained_measure forces weight onto at least settings
distinct candidates. Counting runs asserts that every candidate costs the same,
which fails once a trace of 2N frames competes with a trace of N, so
budgeted_measure normalises by cost. identification_front returns integer
batches that trade parameter precision against separating two models. Each
reports whether its certificate still holds under the constraint.
log det scores the material parameters. A Schur complement scores the model
label. Neither says which question a batch should serve, because they do not
share a scale. pyimr.gain puts every question in nats: the expected
information gain about any subset of coordinates, with the rest treated as
nuisance, is (1/2) log det of a Schur complement of I + M in
prior-standardised coordinates. It is concave, so optimal_measure still
certifies it, and it is finite when M is singular, so a design that cannot
determine every parameter scores low rather than -inf.
from pyimr.gain import Question, expected_gain, gain_criterion
questions = [Question("modulus", (0,)), Question("viscosity", (1,), weight=0.5)]
expected_gain(matrices[0], questions).per_question # nats, unweighted
measure = optimal_measure(np.array(matrices), criterion=gain_criterion(questions))Nats say how much a batch teaches, not whether it settles anything.
runs_to_settle inverts the Gaussian-linear Bayes factor distribution to say
how many runs clear a threshold in favour of the truth at a stated confidence;
runs_to_precision asks the same of one parameter; lack_of_fit_degrees
reports whether the batch leaves any degrees of freedom to discover that every
model in the catalogue is wrong.
pyimr.discriminate scores model discrimination as an integral rather than a
minimum. T-optimality takes min over the rival's parameters, and that inner
problem is multimodal here: a local method landing in the wrong basin returns a
wrong answer with nothing to signal it. expected_log_bayes_factor replaces it
with the expected log Bayes factor between prior-predictive banks, and every
result carries an effective sample size so the silent collapse of a
prior-sample evidence onto one draw is visible. laplace_log_evidence is the
fallback for when it collapses, and is what candidate_log_evidence uses.
screen_models decides which rivals are worth designing for at all. A rival
the records have already decided against by more than decisive nats needs no
experiment, and a design that spends runs separating it is spending them on a
settled matter:
from pyimr.discriminate import screen_models
screen = screen_models(evidences) # one log evidence per model, on the data in hand
screen.live, screen.decided, screen.weightsThe survivors' weights feed gain_criterion as question weights, and the
close pairs are what measure.augmented_information and separability are for:
the model label becomes one more Jacobian column, and what matters is how much
of the difference between two models survives refitting the material.
optimize_design maximises one number. That is the right question only while
the design criteria happen to agree, and on the qSLS study they do not.
Measured over the (R_max, stretch) plane:
| criterion | best design |
|---|---|
identify g against alpha (E-optimality) |
277 µm, stretch 5.0 |
| separate qSLS from SLS — is stiffening real? | 1200 µm, stretch 7.09 |
| separate qSLS from qKV — is relaxation real? | 100 µm, stretch 20 |
The last is the opposite corner from the first, for a physical reason: detecting
a relaxation mode needs a collapse fast enough to excite it, while separating
g from alpha needs a gentle one, because that term of Ze dies as
lambda^-4. No design is best at both, so the answer is a set of designs.
explore_tradeoff traces it. Pass an objective returning several numbers, all
maximised; return a non-finite entry for a design that cannot be run.
from pyimr.pareto import explore_tradeoff
result = explore_tradeoff(criteria, [[50e-6, 1200e-6], [3.0, 20.0]],
evaluations=34, initial=10)
result.front_points # the designs nothing else beats on every criterion
result.front_values # their scores
result.best_for(0) # the design that wins criterion 0 outrightTwo things worth knowing about how it works. The search is ParEGO: one shared archive with a fresh weight vector each iteration, so every expensive evaluation informs every weight — which is what makes it affordable when scoring one design costs hundreds of solves. And the scalarisation is Chebyshev rather than a weighted sum, because a weighted sum only ever returns points on the convex hull of the front. Where two mechanisms genuinely compete the front bulges inward, and there a weighted sum returns the same two endpoints at every weight while the whole interior stays unreachable.
When you score models against each other, let the rival re-fit. T-optimality is a minimum over the rival's parameters, not a difference at its own fitted values; holding them fixed overstates discrimination badly. On this study it inflated qSLS-versus-SLS separation by about an order of magnitude — the honest figure at the present design is 0.668 noise units, meaning SLS imitates qSLS to well within the noise.
SimulationError: maximum number of solver steps was reached is true and
nearly useless: it is the same message whether the trajectory needs a bigger
budget, violates a material's domain, or is so sensitive to its inputs that no
tolerance will be met, and those want opposite responses. pyimr.diagnose runs
the ladder that distinguishes them:
from pyimr.diagnose import diagnose
report = diagnose(config, times)
report.outcome # ok, runaway, ill-conditioned, domain, budget, unresolved
report.summaryIt is several solves, meant for a point that already misbehaved, not for a
sweep. A solve that succeeds is checked too. ill-conditioned means a 1e-9
change in R0 is amplified until the trajectory has shed digits, and the
result is one draw from a sensitive system rather than an answer. budget
means a larger max_steps finishes. domain means the material refused a state
the collapse reached. runaway means no tolerance or solver will help.
unresolved means a loose tolerance completes and the requested one does not,
which is worth a look by hand.
A parameter study is re-run many times while the analysis around it changes.
ResultStore keys on the content of (times, config), so the second run is
free, and it caches failures too, which is most of the value: a point that
exhausts its step budget spends the whole budget before saying so, every time
it is asked.
from pyimr.store import ResultStore
store = ResultStore("~/.cache/pyimr")
result = store.simulate(times, config) # same signature and result as pyimr.simulateThe key cannot see a change to PyIMR itself; pass version= to invalidate a
directory by hand after one.
pyimr.assimilation estimates the initial state of a prepared problem from a
window of observations. four_dvar minimises the strong-constraint 4D-Var
cost with exact gradients from the tangent-linear operator, which is why true
4D-Var is reachable here rather than an ensemble approximation to it.
ensemble_smoother answers the same question from ensemble statistics,
ienks re-linearises about the current estimate each iteration and closes on
the minimum four_dvar finds, and enkf_analysis, ensemble_update and
kalman_analysis are the one-step analyses those are built from. The smoothers
take the PreparedProblem, the times, the observations and a linear
observation operator; noise is a scalar or vector of standard deviations, or
a full covariance. The docstrings give the shapes.
Every workers= argument in the package goes through pyimr.parallel. It
exists because XLA sizes its CPU thread pool from the process's affinity mask
and OpenBLAS does the same at import: sixteen spawn workers on a 128-core host
were measured asking for roughly 6,500 threads between them. worker_pool gives
each worker one core and one thread per library, map_work decides serial or
pooled by timing the first item rather than counting them, and
default_workers reads the affinity mask rather than cpu_count, so a job on a
shared node claims only what the scheduler granted. PYIMR_WORKERS overrides
the default. Set the thread-count environment variables before importing numpy
if you build your own pool; limit_worker_threads does it for you.
pyimr.data covers the step before inference: getting from a measured R(t)
history to the quantities a fit needs.
from pyimr import data
Req = data.equilibrium_radius(R0_m, initial_gas_pressure_pa)
omega_n, beta = data.natural_frequency(R0_m, Req, 2500.0, 0.1)
collapse_times_s, peak_radii_m, peak_times_s = data.collapse_features(
measured_time_s, measured_radius_m
)
data.resolution_convergence(config, times_s, [10, 20, 40])equilibrium_radius inverts the solver's own pressure/radius relation exactly.
natural_frequency linearises Rayleigh-Plesset about Req in a Kelvin-Voigt
medium; it reproduces Minnaert exactly in the gas-only limit and matches the
simulated rebound frequency closely. collapse_features locates interior
extrema with sub-sample parabolic refinement, replacing the manual index windows
of IMR-vanilla calc_3tmins_3Rmaxs. resolution_convergence reports a table for
a ladder you supply, where pyimr.resolution searches for a setting; both scale
deviations by the field's own peak, and both move Mt with Nt only when the
medium is actually solved. Pass (Nt, Mt) pairs to set both.
IMR-vanilla's calc_omega_N is deliberately not ported: it is a scratch script
whose formula treats the gas pressure at Rmax as the equilibrium value,
inflating the stiffness by alpha**(-3*kappa); see upstream.md.
Video processing (calcRofT/) is also out of scope -- that is image analysis,
and scikit-image covers it.