API reference
Rendered live from the installed engine docstrings via mkdocstrings — the docstrings in dmipy-sim and dmipy-fit are the source of truth, so this reference can't drift from the code. For the conceptual entry points see the forward and inverse overviews.
dmipy-sim — forward
Simulation
simulate
simulate(
n_walkers: int,
diffusivity=None,
waveform=None,
geometry=None,
seed: int = 123,
T2: float = None,
return_positions: bool = False,
return_compartments=False,
return_walker_signals: bool = False,
r0=None,
walker_batch_size: int = None,
require_gpu=None,
_allow_oom_backoff: bool = True,
)
Run Monte Carlo diffusion simulation.
Parameters
n_walkers : int Number of random walkers. diffusivity : float, optional Diffusion coefficient in m²/s. Required for standard geometries. Omit for MyelinatedCylinder (D values are in the geometry). waveform : Waveform Gradient waveform. G has shape (n_measurements, n_t, 3). geometry : Geometry Boundary geometry. Provides init_positions() and reflect(). seed : int Master PRNG seed (split into per-walker keys internally). T2 : float, optional Transverse relaxation time in seconds. When set, accumulated per-walker inside the scan body as -dt/T2 each step. return_positions : {False, True, 'full'}, optional False (default): no positions. True: final walker positions, (n_walkers, 3). 'full': per-timestep positions, (n_walkers, n_timesteps, 3) — trajectory export for visualisation/analysis (e.g. combine with return_compartments='full' to select walkers that permeated). Supported for standard geometries including Mesh; not the myelin step-fn paths. return_compartments : {False, 'final', 'full'}, optional Controls compartment-ID output. Default False (no change to return value).
- ``False``: no compartment output.
- ``'final'``: return ``(compartment_origin, compartment_current_final)``
as additional outputs. Both are int32 arrays of shape
``(n_walkers,)``.
- ``'full'``: return ``(compartment_origin, compartment_current_full)``
where ``compartment_current_full`` has shape
``(n_walkers, n_timesteps)`` containing the compartment ID at every
timestep.
array-like of shape (n_walkers, 3), optional
Custom initial walker positions in metres (lab frame, float32).
When provided, geometry.init_positions() is skipped and these
positions are used directly. Useful for mixed initial conditions
(e.g., f·N walkers inside cylinders, (1-f)·N walkers outside)
required for Karger-model validation. Default None (use geometry
default positions).
Compartment integer IDs:
Cylinder,Sphere,Ellipsoid,Box1D: 0 = intra, 1 = extra.MyelinatedCylinder: 0 = intra-axonal, 1 = myelin, 2 = extra-axonal.PackedCylinders: 0 = extra-axonal, 1..N = intra cylinder k (1-indexed).
int, optional
If set and smaller than n_walkers, the run is split into walker
chunks of this size, run one at a time, and recombined. Peak device
memory is bounded to one chunk — use this on a small GPU. Each chunk
uses an independent sub-seed, so the ensemble signal is statistically
identical to a single-shot run (not bit-identical). Default None
(all walkers at once).
require_gpu : {None, True, False}, optional
GPU guard against a silent CPU fallback. True raises if no GPU is
visible; False opts out (e.g. a CPU float64 reference check);
None (default) warns when a large run is about to use the CPU.
Returns
signals : np.ndarray of shape (n_measurements,), float32
Normalised signal: Re(
simulate_cpmg
simulate_cpmg(
n_walkers,
diffusivity,
waveform,
geometry,
*,
T2=None,
seed=123,
walker_batch_size=None,
require_gpu=None
)
Multi-echo CPMG signal from a SINGLE diffusion walk.
Walks the spin ensemble once through the full CPMG train (ideal instantaneous
180° refocusing is encoded as the sign flips of waveform.G) and samples the
ensemble signal Re<exp(iφ)·exp(log_w)> at each echo time. This is the
ordinary forward model: one pass through the train, nothing cached or reused.
Build waveform with :func:dmipy_sim.cpmg
(which sets echo_indices).
Parameters
n_walkers : int
diffusivity : float or None
Bulk diffusivity (m²/s); omit for MyelinatedCylinder (D in the geometry).
waveform : Waveform
A multi-echo waveform carrying echo_indices (e.g. from cpmg).
geometry : Geometry
T2 : float, optional
Transverse relaxation time (s), accumulated per-walker in the walk.
seed, walker_batch_size, require_gpu : see :func:simulate.
Returns
signals : np.ndarray, shape (n_echoes, n_measurements), float32 Signal at each echo (echo k = k·TE), one column per gradient direction.
Geometries
Sphere
Bases: Geometry
Reflecting sphere of given radius centred at the origin.
Parameters
radius : float Sphere radius in metres. surface_relaxivity_t2 : float, optional Surface relaxivity ρ₂ in m/s. When set, boundary collisions reduce the walker magnetisation weight by exp(-2·ρ₂·d_perp/D). T2_surface = R / (3·ρ) for a sphere (S/V = 3/R). Default None. permeability : float, optional Membrane permeability κ in m/s. When set, each boundary crossing is probabilistic: the walker transmits with p = min(1, 2κ·d_perp/D) and reflects otherwise. Enables bidirectional exchange — walkers may be inside or outside the sphere at any time. Default None (fully reflecting wall). Exchange time τ = R / (3κ).
volume
volume() -> float
Volume of the sphere: (4/3)·π·R³ (m³).
surface_area
surface_area() -> float
Surface area of the sphere: 4·π·R² (m²).
classify_position
classify_position(r: ndarray) -> jnp.ndarray
Compartment ID: 0=intra (|r| < R), 1=extra (|r| >= R).
init_positions
init_positions(n_walkers, key)
Uniform sampling inside sphere via rejection (CPU numpy).
reflect
reflect(r, step)
Specular reflection off sphere boundary with multiple reflections.
Matches disimpy's while-loop convention: decomposes step into a unit direction + scalar remaining length, then iterates up to MAX_ITER times. At each iteration the distance d to the boundary is computed; if d < remaining the walker is moved to the boundary, the direction is specularly reflected, and remaining is decremented by d + epsilon (epsilon nudges the walker just inside the surface to avoid numerical re-intersection). Uses jax.lax.scan for JAX-compilable fixed iteration.
reflect_with_log_weight
reflect_with_log_weight(r, step, rho_over_D)
Reflect and accumulate per-collision surface-relaxation log-weight.
Uses the same quadratic intersection as reflect(), but also computes the perpendicular penetration depth d_perp = (remaining - d) * cos(α) at each collision, where cos(α) = sqrt(disc) / R.
Returns
r_out : jnp.ndarray, shape (3,) Final walker position. dlog_w : jnp.float32 Log-weight decrement: -2 * rho_over_D * sum(d_perp). For sphere S/V = 3/R → T2_surface = R / (3·ρ).
permeate
permeate(r, step, kappa_over_D, rho_over_D, perm_key)
Probabilistic membrane crossing (Powles 2004) + optional relaxivity.
At each boundary crossing the walker transmits with probability
p = min(1, 2 · κ/D · d_perp)
and reflects otherwise. When rho_over_D > 0 a Brownstein-Tarr weight is applied on reflection only.
The geometry is bidirectional: walkers may be inside (|r| < R) or outside (|r| > R); the appropriate intersection root is selected automatically.
Single-event-per-step approximation. Requires σ/R < 0.1. Exchange time τ = R / (3κ) (V/κS = R/3).
Regime/step-size note: the crossing is correct and step-robust in any fast/slow regime (validated by the closed-system permeability tests: high-κ → free diffusion, monotone-in-κ, etc.). A SINGLE permeable object in an OPEN domain, however, is not a well-mixed reservoir: walkers that exit re-enter, so the apparent residence time is biased ABOVE τ=R/(3κ) (a step-independent effect that grows with time and is larger for lower-dimensional exteriors: 2D cylinder > 3D sphere). For exchange / residence-time work use a periodic packed geometry (a proper reservoir), where τ=R/(3κ) holds; do not compare a single open-domain object's f_inside(t) to the well-mixed τ.
Parameters
r : (3,) float32, current position step : (3,) float32, proposed displacement kappa_over_D : float32, κ/D baked in by make_step_fn rho_over_D : float32, ρ/D (0.0 if no surface relaxivity) perm_key : JAX PRNGKey for the Bernoulli draw
Returns
r_new : (3,) float32, new position dlog_w : float32, log-weight decrement (≤ 0; 0.0 on transmission)
Cylinder
Bases: Geometry
Reflecting infinite cylinder of given radius and orientation.
Restriction acts in the plane perpendicular to orientation.
Walkers move freely along the cylinder axis.
Parameters
radius : float Cylinder inner radius in metres. orientation : array-like of shape (3,) Cylinder axis direction (normalised internally). surface_relaxivity_t2 : float, optional Surface relaxivity ρ₂ in m/s. When set, each boundary collision reduces the walker magnetisation weight by exp(-2·ρ₂·d_out/D), where d_out is the step length that would have exited the cylinder. This implements the surface-T2 model: 1/T2_eff = 1/T2_bulk + ρ₂·S/V with S/V = 2/R for a cylinder. Default None (no surface relaxation). permeability : float, optional Membrane permeability κ in m/s. When set, each boundary crossing is probabilistic: the walker transmits through the wall with probability p = min(1, 2κ·d_perp/D) and reflects otherwise. Enables bidirectional exchange — walkers may be inside or outside the cylinder at any time. Default None (fully reflecting wall).
init_positions
init_positions(n_walkers, key)
Uniform sampling in circular cross-section.
reflect
reflect(r, step)
Specular reflection off cylinder boundary with multiple reflections.
Works in the cylinder frame (orientation → z-axis via _R). The z-axis is free; reflection acts only in the x-y plane (the restricted cross-section). Uses the same unit-direction + scalar-remaining convention as disimpy, with an epsilon nudge inward after each reflection. jax.lax.scan gives a fixed iteration count that is JAX-compilable.
reflect_with_log_weight
reflect_with_log_weight(r, step, rho_over_D)
Specular reflection + surface-relaxation log-weight decrement.
Identical to reflect() but also accumulates the perpendicular outgoing step depths across all boundary collisions within one timestep, and converts them to a magnetisation log-weight decrement:
Δlog_w = -2 · ρ_over_D · Σ d_perp_i
where d_perp_i = (remaining - d) · cos(α) is the perpendicular depth of wall penetration at collision i (cos(α) = dot(d_hat, n_out) at hit point = sqrt(disc)/R), and rho_over_D = ρ₂/D. The factor of 2 comes from the Brownstein-Tarr boundary condition (ρ · ∂M/∂n = D · M) in the partially-absorbing Monte Carlo formulation. Using d_perp (not the full remaining step d_out) is required for the formula to reproduce the correct S/V scaling; numerical verification shows that the coefficient C = π/2 ≈ 1.571 for d_out and C = 2 for d_perp.
Parameters
r : (3,) float32, current walker position (lab frame) step : (3,) float32, proposed displacement (lab frame) rho_over_D : float32, ρ₂/D baked in by make_step_fn
Returns
r_new : (3,) float32, new position (lab frame) dlog_w : float32, log-weight decrement (≤ 0)
permeate
permeate(r, step, kappa_over_D, rho_over_D, perm_key)
Probabilistic membrane crossing (Powles 2004) + optional relaxivity.
At each wall crossing the walker transmits through with probability
p = min(1, 2 · κ/D · d_perp)
and reflects otherwise (specular, same as reflect()). When rho_over_D > 0 a Brownstein-Tarr weight is applied on reflection:
Δlog_w = −2 · ρ/D · d_perp
No weight is applied on transmission. The geometry is bidirectional: walkers may be inside (|r_xy| < R) or outside (|r_xy| > R); the appropriate intersection root is selected automatically.
Single-event-per-step approximation: at most one permeability decision per timestep (consistent with the PackedCylinders exterior method). Requires σ/R < 0.1 for accurate results.
Parameters
r : (3,) float32, current position (lab frame) step : (3,) float32, proposed displacement (lab frame) kappa_over_D : float32, κ/D baked in by make_step_fn rho_over_D : float32, ρ/D (0.0 if no surface relaxivity) perm_key : JAX PRNGKey for the Bernoulli draw
Returns
r_new : (3,) float32, new position dlog_w : float32, log-weight decrement (≤ 0; 0.0 on transmission)
classify_position
classify_position(r: ndarray) -> jnp.ndarray
Compartment ID: 0=intra (|r_xy| < R), 1=extra (|r_xy| >= R).
The check is performed in the cylinder frame (r_xy is the component perpendicular to the cylinder axis).
volume
volume(L: float = 1.0) -> float
Volume of the cylinder: π·R²·L (m³).
Parameters
L : float, optional Cylinder length in metres. Default 1.0 (returns per-unit-length volume, i.e. the cross-sectional area π·R²).
surface_area
surface_area(
L: float = 1.0, include_caps: bool = False
) -> float
Lateral surface area of the cylinder: 2π·R·L (m²).
Caps are excluded by default because the cylinder is modelled as infinite (periodic along its axis) and caps are irrelevant for permeability and relaxivity calculations.
Parameters
L : float, optional Cylinder length in metres. Default 1.0 (returns per-unit-length lateral area, i.e. the circumference 2·π·R). include_caps : bool, optional If True, add the two circular end caps 2·π·R². Default False.
Mesh
Bases: Geometry
Reflecting/permeable triangular-mesh geometry (see module docstring).
Parameters
vertices : (n_vertices, 3) array-like, metres
faces : (n_faces, 3) array-like of int
Triangle vertex-index triples.
periodic : bool or (bool, bool, bool)
Wrap walkers periodically along the given axes (default False = closed mesh).
When any axis is periodic you must pass voxel_min/voxel_max defining
the periodic box (the mesh bbox is usually slightly larger than the true
period, so it is not a safe default).
voxel_min, voxel_max : (3,) array-like, optional
The simulation box. Defaults to the mesh bounding box (closed meshes only).
feature_radius : float, optional
Characteristic feature size (e.g. a cell/pore radius), used to size the
diffusion sub-step (step ~ feature_radius/6, or /25 when permeable)
and the grid. Defaults to half the smallest box side; pass the real cell
radius for packed substrates, otherwise the step may be too coarse.
surface_relaxivity_t2 : float, optional
Surface relaxivity ρ₂ (m/s), symmetric (same on both sides of the wall).
Applies a Brownstein–Tarr weight at the wall. For a side-dependent ρ, use
intra=/extra= instead.
permeability : float or dict, optional
Membrane permeability κ (m/s). A float is symmetric (same both directions,
the default). A dict {"intra_to_extra": κ_out, "extra_to_intra": κ_in}
makes it direction-dependent — note asymmetric κ is a pump (net flux, not
passive equilibrium). None → impermeable.
intra, extra : dict, optional
Per-compartment surface properties for the intra (inside a cell) and extra
sides — currently surface_relaxivity_t2 only, e.g.
intra={"surface_relaxivity_t2": 5e-6}, extra={"surface_relaxivity_t2": 1e-6}:
a spin hitting the wall from inside vs outside then takes a different
relaxivity weight. (Per-compartment diffusivity/T2 is a later layer.)
orientation : (3,) array-like, optional
Direction, in the scanner frame (B0 = +z), along which the mesh's native
+z axis (e.g. a periodic / fibre axis) is placed in the bore. Applied as
an acquisition rotation (the gradient is rotated into the mesh frame in
simulate), so the walk itself is unchanged. The in-plane rotation is
arbitrary — pass R for meshes whose in-plane orientation matters.
R : (3, 3) array-like, optional
Explicit mesh→lab rotation matrix (mutually exclusive with orientation).
cell_size : float, optional
Acceleration-grid cell size. Defaults to 4 * step (safe for the 27-cell
neighbourhood). Larger = fewer/denser cells; must be ≥ the maximum step.
classify_position
classify_position(r)
Compartment tag: 0 = interior (inside a cell), 1 = exterior.
init_positions
init_positions(n_walkers, key, intra=True)
Seed walkers inside (intra=True) or outside the cells by grid rejection.
quality_report
quality_report(verbose=True)
Surface-resolution diagnostics + per-effect accuracy verdict.
Returns a dict; also prints a table when verbose. Uses trimesh (if
installed) to add watertight / component info. The key number is
edge_feature_ratio = median edge / feature_radius: permeability needs
it <~ 0.05 to reach the MC noise floor; diffusion and surface
relaxivity are fine at much coarser ratios.
from_ply
classmethod
from_ply(path, scale=1.0, recenter=False, **kwargs)
Construct a Mesh directly from a mesh file (see :func:load_ply).
Waveforms & encoding
pgse
pgse(
delta,
DELTA,
G_magnitude,
bvecs,
n_t,
slew_rate=DEFAULT_SLEW_RATE,
)
Build a PGSE gradient waveform.
Slew-limited (realizable, trapezoidal lobes) by default -- dmipy-sim is the
forward truth. Pass slew_rate=np.inf for the idealized instantaneous
(square) limit (e.g. A/B against an analytic solution); slew_rate must be
a positive T/m/s value or np.inf (None is rejected). b is set later
via :func:set_b; slew-limiting changes the lobe SHAPE (hence the restricted-
diffusion signal), which is the point.
Parameters
delta : float Gradient pulse duration in seconds. DELTA : float Diffusion time in seconds (centre-to-centre of gradient pulses). G_magnitude : float or array of shape (n_measurements,) Gradient amplitude in T/m. Use set_b() to scale to target b-values. bvecs : array of shape (n_measurements, 3) Unit gradient direction vectors. n_t : int Number of time points. Total duration = DELTA + delta.
Returns
Waveform
set_b
set_b(waveform, b_target)
Return a new Waveform scaled so each measurement has the given b-value.
Parameters
waveform : Waveform
b_target : float or array of shape (n_measurements,)
Target b-values in s/m² (SI units), consistent with calc_b.
Typical clinical values: 1e8–3e9 s/m² (= 100–3000 s/mm²).
.. warning::
A common mistake is passing b-values in **s/mm²** (e.g. 1000) instead
of **s/m²** (e.g. 1e9). ``set_b`` will silently produce gradients
that are 1000× too small, giving essentially b≈0 signals.
Convert: ``b_si = b_mm2 * 1e6``.
Returns
Waveform with scaled G.
dmipy-fit — inverse
Modelling framework
MultiCompartmentModel
Bases: MultiCompartmentModelProperties
The MultiCompartmentModel class allows to combine any number of CompartmentModels and DistributedModels into one combined model that can be used to fit and simulate dMRI data.
Parameters
models : list of N CompartmentModel instances, the models to combine into the MultiCompartmentModel. parameter_links : list of iterables (model, parameter name, link function, argument list), deprecated, for testing only.
fit
fit(
acquisition_scheme,
data,
mask=None,
solver="brute2fine",
Ns=5,
maxiter=300,
N_sphere_samples=30,
use_parallel_processing=False,
number_of_processors=None,
batch_size=None,
loss_fn=None,
sigma_x0=None,
sigma_range=(0.001, 0.5),
)
The main data fitting function of a MultiCompartmentModel.
This function can fit it to an N-dimensional dMRI data set, and returns a FittedMultiCompartmentModel instance that contains the fitted parameters and other useful functions to study the results.
No initial guess needs to be given to fit a model, but a partial or complete initial guess can be given if the user wants to have a solution that is a local minimum close to that guess. The parameter_initial_guess input can be created using parameter_initial_guess_to_parameter_vector().
A mask can also be given to exclude voxels from fitting (e.g. voxels that are outside the brain). If no mask is given then all voxels are included.
An optimization approach can be chosen as either 'brute2fine' or 'mix'. - Choosing brute2fine will first use a brute-force optimization to find an initial guess for parameters without one, and will then refine the result using gradient-descent-based optimization.
Note that given no initial guess will make brute2fine precompute an global parameter grid that will be re-used for all voxels, which in many cases is much faster than giving voxel-varying initial condition that requires a grid to be estimated per voxel.
- Choosing mix will use the recent MIX algorithm based on separation of linear and non-linear parameters. MIX first uses a stochastic algorithm to find the non-linear parameters (non-volume fractions), then estimates the volume fractions while fixing the estimates of the non-linear parameters, and then finally refines the solution using a gradient-descent-based algorithm.
The fitting process can be parallelized across voxels using stdlib concurrent.futures. Pass use_parallel_processing=True to enable it. The algorithm will automatically use all cores in the machine, unless otherwise specified in number_of_processors.
Data with multiple TE are normalized in separate segments using the b0-values according that TE.
Parameters
acquisition_scheme : PGSEAcquisitionScheme instance, An acquisition scheme that has been instantiated using dMipy. data : N-dimensional array of size (N_x, N_y, ..., N_dwis), The measured DWI signal attenuation array of either a single voxel or an N-dimensional dataset. mask : (N-1)-dimensional integer/boolean array of size (N_x, N_y, ...), Optional mask of voxels to be included in the optimization. solver : string, Selection of optimization algorithm. - 'brute2fine' to use brute-force optimization. - 'mix' to use Microstructure Imaging of Crossing (MIX) optimization. Ns : integer, for brute optimization, decised how many steps are sampled for every parameter. maxiter : integer, for MIX optimization, how many iterations are allowed. N_sphere_samples : integer, for brute optimization, how many spherical orientations are sampled for 'mu'. use_parallel_processing : bool, whether or not to use parallel processing (default False). number_of_processors : integer, number of processors to use for parallel processing. Defaults to the number of processors in the computer according to cpu_count(). sigma_x0 : float or None, Initial guess for sigma (noise standard deviation in normalised signal units = 1/SNR0). When provided, sigma is jointly optimized with the diffusion model parameters (fittable sigma mode). Only supported with solver='jax'. Defaults to None (fixed sigma). sigma_range : tuple (float, float), (lower, upper) bounds for sigma during optimization. Only used when sigma_x0 is not None. Default (0.001, 0.5) covers SNR 2–1000.
Returns
FittedCompartmentModel: class instance that contains fitted parameters, Can be used to recover parameters themselves or other useful functions.
simulate_signal
simulate_signal(
acquisition_scheme, parameters_array_or_dict
)
Function to simulate diffusion data for a given acquisition_scheme and model parameters for the MultiCompartmentModel.
Parameters
acquisition_scheme : PGSEAcquisitionScheme instance, An acquisition scheme that has been instantiated using dMipy model_parameters_array : 1D array of size (N_parameters) or N-dimensional array the same size as the data. The model parameters of the MultiCompartmentModel model.
Returns
E_simulated: 1D array of size (N_parameters) or N-dimensional array the same size as x0. The simulated signal of the microstructure model.
Acquisition scheme
acquisition_scheme_from_bvalues
acquisition_scheme_from_bvalues(
bvalues,
gradient_directions,
delta=None,
Delta=None,
TE=None,
min_b_shell_distance=50000000.0,
b0_threshold=10000000.0,
)
Creates an acquisition scheme object from bvalues, gradient directions, pulse duration \(\delta\) and pulse separation time \(\Delta\).
Parameters
bvalues: 1D numpy array of shape (Ndata) bvalues of the acquisition in s/m^2. e.g., a bvalue of 1000 s/mm^2 must be entered as 1000 * 1e6 s/m^2 gradient_directions: 2D numpy array of shape (Ndata, 3) gradient directions array of cartesian unit vectors. delta: float or 1D numpy array of shape (Ndata) if float, pulse duration of every measurements in seconds. if array, potentially varying pulse duration per measurement. Delta: float or 1D numpy array of shape (Ndata) if float, pulse separation time of every measurements in seconds. if array, potentially varying pulse separation time per measurement. min_b_shell_distance : float minimum bvalue distance between different shells. This parameter is used to separate measurements into different shells, which is necessary for any model using spherical convolution or spherical mean. b0_threshold : float bvalue threshold for a measurement to be considered a b0 measurement.
Returns
PGSEAcquisitionScheme: acquisition scheme object contains all information of the acquisition scheme to be used in any microstructure model.
White matter
t2_spectrum_mwf
t2_spectrum_mwf(
signal,
echo_times,
T2_grid=None,
cutoff=0.025,
reg="x2",
x2_factor=1.02,
)
Classic regularised NNLS T2-spectrum MWF (the standard analysis).
Fits a non-negative \(T_2\) spectrum to the multi-echo decay assuming ideal
exp(-t/T2) bases (the instantaneous-pulse assumption), and returns the myelin
water fraction = spectral weight below cutoff (s). Returns
(mwf, T2_grid, spectrum).
reg selects the regularisation: 'x2' (default) chooses the weight
per signal by the chi-square criterion (Whittall--MacKay / Prasloski /
DECAES standard; smoothest fit within x2_factor of the unregularised
misfit), which is the literature-trusted, noise-robust choice. A float
instead applies that fixed zero-order Tikhonov weight (legacy behaviour).