Parallel Ensemble Simulations Interface

abstract type AbstractEnsembleProblem <: SciMLBase.AbstractSciMLProblem

Base interface for ensemble problems.

An AbstractEnsembleProblem describes many related solves generated from a template problem. Concrete subtypes should expose the template problem, a trajectory-generation hook, an output hook, a batch reduction hook, and any initial reduction state needed by ensemble solvers. The standard concrete implementation is EnsembleProblem.

struct EnsembleProblem{T, T2, T3, T4, T5} <: SciMLBase.AbstractEnsembleProblem

Container for a template problem and the user hooks used to run an ensemble of related SciML solves.

An EnsembleProblem is solved by repeatedly calling prob_func(prob, ctx) to construct the trajectory-specific problem, solving that problem with the requested numerical algorithm, passing the result through output_func(sol, ctx), and combining batches with reduction(u, data, I). The ctx argument is an EnsembleContext that identifies the trajectory and carries per-trajectory RNG state when rng or seed is supplied to solve.

Constructor

EnsembleProblem(
    prob::AbstractSciMLProblem;
    output_func = (sol, ctx) -> (sol, false),
    prob_func = (prob, ctx) -> prob,
    reduction = (u, data, I) -> (append!(u, data), false),
    u_init = [], safetycopy = prob_func !== DEFAULT_PROB_FUNC
)

Positional Arguments

• prob: The canonical problem used as the template for each trajectory.

Keyword Arguments

• prob_func: A function (prob, ctx) that modifies the problem for each trajectory. ctx is an EnsembleContext providing ctx.sim_id (unique id 1:trajectories), ctx.repeat (rerun counter, starts at 1), ctx.rng (per-trajectory RNG or nothing), ctx.sim_seed, and ctx.master_rng. prob_func must preserve the problem type of prob; for example, a JumpProblem must remain a JumpProblem, an ODEProblem must remain an ODEProblem.
• output_func: A function (sol, ctx) that determines what is saved from each trajectory. It returns (out, rerun), where out is stored in the batch output and rerun requests that the same trajectory be rerun with ctx.repeat incremented.
• reduction: A function (u, data, I) that combines the current accumulator u with the outputs data from the trajectory index range I. It returns (new_u, converged), where converged=true stops the ensemble early.
• u_init: The initial accumulator passed to reduction. When nothing, the accumulator is initialized from the first batch output.
• safetycopy: Determines whether a safety deepcopy is called on the prob before the prob_func. By default, this is true for any user-given prob_func, as without this, modifying the arguments of something in the prob_func, such as parameters or caches stored within the user function, are not necessarily thread-safe. If you know that your function is thread-safe, then setting this to false can improve performance when used with threads. For nested problems, e.g., SDE problems with custom noise processes, deepcopy might be insufficient. In such cases, use a custom prob_func.

Example

function prob_func(prob, ctx)
    remake(prob, u0 = randn(ctx.rng, length(prob.u0)))
end

output_func(sol, ctx) = (sol[end, 2], false)
ensemble_prob = EnsembleProblem(prob; prob_func, output_func)

WeightedEnsembleProblem(ensembleprob, weights)
WeightedEnsembleProblem(args...; weights, kwargs...)

Associate an EnsembleProblem with one weight per trajectory.

The direct constructor wraps an existing ensemble problem. The keyword constructor forwards args and kwargs to EnsembleProblem, verifies that weights sum to one, and verifies that their length matches the generated problem collection.

Arguments

• ensembleprob: The ensemble problem whose trajectories are weighted.
• weights: An AbstractVector containing one weight per trajectory.

Keywords

• weights: Required by the forwarding constructor.
• kwargs...: Forwarded to EnsembleProblem.

Fields

• ensembleprob: The wrapped ensemble problem.
• weights: The trajectory weights.

Returns

• WeightedEnsembleProblem: A wrapper that delegates ordinary ensemble-problem properties and exposes weights.

Throws

• AssertionError: The forwarding constructor throws when the weights do not sum to one or their length does not match the generated problem collection.

Example

ensembleprob = EnsembleProblem([:first, :second])
weighted = WeightedEnsembleProblem(ensembleprob, [0.25, 0.75])
weighted.weights

DEFAULT_PROB_FUNC(prob, ctx)

Default prob_func for EnsembleProblem. It returns the template problem unchanged for every trajectory, so all trajectory-specific variation must come from solver randomness, callbacks, or other solve-time state.

DEFAULT_OUTPUT_FUNC(sol, ctx)

Default output_func for EnsembleProblem. It stores the full solution object from each trajectory and returns false for the rerun flag, meaning the trajectory is accepted on the first completion.

DEFAULT_REDUCTION(u, data, I)

Default ensemble reduction. It appends the current batch data to the accumulator u and returns false for the convergence flag, so the ensemble continues until all requested trajectories have been run.

Base.propertynames(prob::WeightedEnsembleProblem) -> Tuple

Return the delegated property names of prob.ensembleprob, followed by :weights.

This method defines the properties advertised to reflection and tab completion. The returned tuple may include :ensembleprob when the wrapped problem advertises it.

Base.getproperty(prob::WeightedEnsembleProblem, name::Symbol)

Return prob.weights, the wrapped prob.ensembleprob, or a property delegated to the wrapped ensemble problem.

Arguments

• prob: The weighted ensemble wrapper.
• name: :weights, :ensembleprob, or a property supported by the wrapped problem.

Returns

• The stored weights, wrapped problem, or delegated property value.

Throws

Errors from getproperty(prob.ensembleprob, name) propagate for unsupported delegated properties.

sim = solve(enprob, alg, ensemblealg = EnsembleThreads(), kwargs...)

Solve an AbstractEnsembleProblem by running many trajectories of the template problem.

alg is the numerical algorithm used for each generated trajectory, while ensemblealg chooses the scheduling backend. Keywords that are not consumed by the ensemble layer are forwarded to each inner solve call, so common solver keywords such as saveat, tolerances, callbacks, and progress options retain their usual meaning for each trajectory.

The special ensemble keywords are:

• trajectories: Required number of trajectories to run.
• batch_size: Number of trajectories processed before prob.reduction is called. Defaults to trajectories.
• pmap_batch_size: Batch size passed to pmap for distributed ensemble algorithms. Defaults to div(batch_size, 100) > 0 ? div(batch_size, 100) : 1.
• progress_aggregate: Whether per-trajectory progress messages are aggregated into a single total progress message. Defaults to true.
• seed: Master seed for reproducible ensemble solves. Pre-generates per-trajectory seeds.
• rng: Master RNG for reproducible ensemble solves. Takes priority over seed.
• rng_func: Custom per-trajectory RNG factory (ctx::EnsembleContext) -> AbstractRNG. Defaults to default_rng_func which seeds the TaskLocalRNG.

In threaded modes, ctx.master_rng is shared across tasks. A custom rng_func must not mutate it unless that mutation is thread-safe.

EnsembleContext{S, R, M}

Contextual information about the current trajectory within an ensemble solve.

An EnsembleContext is passed to rng_func, prob_func(prob, ctx), and output_func(sol, ctx). It is the stable interface for selecting trajectory-specific data without relying on global counters or backend-specific worker state.

Fields

• sim_id::Int: Unique trajectory index in 1:trajectories.
• repeat::Int: Rerun counter, starting at 1 and incremented when output_func requests a rerun.
• worker_id::Int: 0 for serial and threaded execution, or Distributed.myid() for distributed execution.
• sim_seed::S: Pre-generated seed for this trajectory, or nothing.
• rng::R: Per-trajectory RNG created by rng_func, or nothing while rng_func itself is running.
• master_rng::M: User-provided master RNG, or nothing. Distributed modes set this to nothing to avoid serializing mutable RNG state to workers.

generate_sim_seeds(rng, seed, trajectories)

Pre-generate an array of per-trajectory seeds from a master RNG.

If rng is provided it is used directly; otherwise Random.Xoshiro(seed) is constructed. The returned UInt64 values are stored in each EnsembleContext as ctx.sim_seed so serial, threaded, and distributed ensemble execution can create reproducible per-trajectory RNG state.

default_rng_func(ctx::EnsembleContext)

Default per-trajectory RNG factory.

When ctx.sim_seed is available, this seeds Julia's task-local default RNG and returns Random.default_rng(). ctx.rng is nothing when this function is called, since the returned RNG is what will be stored as ctx.rng for prob_func and output_func.

abstract type BasicEnsembleAlgorithm <: SciMLBase.EnsembleAlgorithm

Base interface for built-in ensemble execution algorithms. A BasicEnsembleAlgorithm chooses how trajectories from an EnsembleProblem are scheduled: serially, on Julia threads, across Julia distributed workers, or with a split distributed/threaded strategy.

Concrete subtypes should document their execution backend, worker setup requirements, serialization assumptions, and random number behavior. They are passed as the ensemble algorithm in calls such as solve(ensembleprob, alg, ensemblealg; trajectories, kwargs...), while alg continues to select the numerical solver for each generated trajectory.

struct EnsembleSerial <: SciMLBase.BasicEnsembleAlgorithm

Ensemble execution algorithm that runs trajectories serially in the current Julia process.

EnsembleSerial has the lowest scheduling complexity and is useful for debugging, deterministic single-task execution, or small ensembles where parallel overhead dominates the solve time. Trajectories still receive distinct EnsembleContext values and can use per-trajectory RNG state when seed or rng is supplied to solve.

struct EnsembleThreads <: SciMLBase.BasicEnsembleAlgorithm

Ensemble execution algorithm that schedules trajectories on Julia threads in the current process.

EnsembleThreads is the default basic ensemble backend. It provides shared-memory parallelism with low overhead, but user hooks such as prob_func, output_func, and rng_func must be thread-safe. In particular, mutable objects captured by closures or shared through the template problem should not be mutated unless each trajectory receives independent storage.

struct EnsembleDistributed <: SciMLBase.BasicEnsembleAlgorithm

Ensemble execution algorithm that distributes trajectory batches with Distributed.pmap.

EnsembleDistributed uses the available Julia worker processes; add workers with Distributed.addprocs before solving. The template problem, algorithms, callbacks, and ensemble hooks must be serializable and available on each worker. This backend is usually appropriate when each trajectory is relatively expensive, when trajectories allocate heavily, or when work should be spread across multiple machines. It can also help on a single machine when separate worker processes and garbage collectors outweigh distributed communication overhead.

struct EnsembleSplitThreads <: SciMLBase.BasicEnsembleAlgorithm

Ensemble execution algorithm that combines distributed workers with threaded execution inside each worker process.

EnsembleSplitThreads uses the Julia distributed worker setup provided by the caller and then runs threaded batches on each worker. It is intended for node-based cluster layouts, for example one Julia worker per node with multiple local threads per worker. The same serialization requirements as EnsembleDistributed apply, and threaded hooks must also be thread-safe.

abstract type AbstractEnsembleEstimator <: SciMLBase.AbstractSciMLAlgorithm

Base interface for algorithms that estimate ensemble quantities during an ensemble solve.

Concrete estimator algorithms can be used by ensemble workflows to decide when a Monte Carlo estimate has converged, how many trajectories are required, or how batch reductions should be interpreted. Subtypes should document the statistic they estimate, the stopping criterion, and the reduction state they require.

struct EnsembleSolution{T, N, S} <: SciMLBase.AbstractEnsembleSolution{T, N, S}

Concrete solution container returned by ensemble solves.

The u field stores the accepted output for each trajectory after applying the ensemble problem's output_func. When output_func returns full SciML solutions, u is a collection of solution objects; when it returns reduced values, u stores those reduced values instead. elapsedTime records the wall time spent in the ensemble solve, converged records whether the reduction reported early convergence, and stats stores merged inner-solver statistics when available.

EnsembleSolution supports indexing and iteration through its stored u collection via the AbstractEnsembleSolution interface.

struct EnsembleTestSolution{T, N, S} <: SciMLBase.AbstractEnsembleSolution{T, N, S}

Container for ensemble runs against problems with analytic reference solutions.

EnsembleTestSolution extends the ordinary ensemble solution data with strong and weak error summaries computed by calculate_ensemble_errors. The u field stores the trajectory outputs, while errors, weak_errors, error_means, and error_medians store the per-error-key diagnostics collected across the ensemble. elapsedTime and converged carry the same meaning as in EnsembleSolution.

struct WeightedEnsembleSolution{T1<:SciMLBase.AbstractEnsembleSolution, T2<:Number}

Solution wrapper that associates an ensemble solution with trajectory weights.

The number of weights must match the number of stored trajectories. Weighted ensemble analysis utilities use the weights to form weighted summary statistics without changing the underlying unweighted solution object.

calculate_ensemble_errors(sim::AbstractEnsembleSolution; kwargs...)
calculate_ensemble_errors(trajectories; elapsedTime = 0.0, converged = false,
    weak_timeseries_errors = false, weak_dense_errors = false)

Collect the strong and weak errors from an ensemble whose trajectories include analytic reference solutions, returning an EnsembleTestSolution.

Strong errors are collected from each trajectory's errors field and summarized by their mean and median. The final-time weak error is always calculated. Set weak_timeseries_errors = true to calculate weak errors at the saved time steps, or weak_dense_errors = true to calculate them on a 100-point interpolation grid.

When an AbstractEnsembleSolution is passed, its elapsedTime and converged fields are preserved. The trajectory-collection form accepts those values as keyword arguments.

EnsembleAnalysis

Namespace for summary statistics over SciMLBase.AbstractEnsembleSolution trajectories. Import it with using SciMLBase.EnsembleAnalysis; its functions provide componentwise, timestep, timepoint, and time-series summaries, including weighted covariance operations where supported.

get_timestep(sim, i)

Return a lazy iterator over sol.u[i] for every trajectory in an ensemble solution.

This is a step-index based accessor: it assumes that the ith saved value is the quantity to compare across trajectories. That is appropriate for fixed-step solutions or ensembles saved with common saveat values. Use get_timepoint when trajectories should be compared at a physical time through interpolation.

get_timepoint(sim, t)

Return a lazy iterator over sol(t) for every trajectory in an ensemble solution.

This is a time-point based accessor: each trajectory is evaluated at the same independent-variable value t, using the solution's callable interpolation interface. Use get_timestep when comparing the same saved index instead of the same physical time.

componentwise_vectors_timestep(sim, i)

Collect the values at saved step index i into componentwise trajectory vectors.

For scalar-valued trajectories, the result is a vector of scalar values. For array-valued trajectories, the result is a vector whose entries contain the values of one state component across all trajectories, preserving the component layout needed by the summary-statistic helpers.

componentwise_vectors_timepoint(sim, t)

Collect interpolated values at time t into componentwise trajectory vectors.

For scalar-valued trajectories, the result is a vector of scalar values. For array-valued trajectories, the result is a vector whose entries contain the values of one state component across all trajectories, using the same component layout as componentwise_vectors_timestep.

componentwise_mean(A)

Compute the arithmetic mean of a collection of scalar or array-valued observations componentwise.

For array-valued observations, every observation must have compatible axes and the returned mean has the same shape as an observation. The collection must be nonempty.

Arguments

• A: An iterable of scalar or array-valued observations to average.

Returns

• mean: The arithmetic mean. This is a scalar for scalar observations and an array with the observation shape for array-valued observations.

Examples

componentwise_mean([[1, 2], [3, 4]]) # [2.0, 3.0]

componentwise_meanvar(A; bessel)

Compute the arithmetic mean and variance of a collection of scalar or array-valued observations componentwise using Welford's algorithm.

For array-valued observations, every observation must have compatible axes and the returned mean and variance have the same shape as an observation. The collection must be nonempty. Fewer than two observations produce NaN instead of a (mean, variance) tuple.

Arguments

• A: An iterable of scalar or array-valued observations from which to compute the statistics.

Keywords

• bessel::Bool = true: If true, use the sample-variance denominator n - 1. If false, use the population-variance denominator n.

Returns

• (mean, variance): A tuple containing the componentwise arithmetic mean and variance. Each entry is a scalar for scalar observations or an array with the observation shape for array-valued observations.

Examples

observations = [[1, 2], [3, 4]]

componentwise_meanvar(observations) # ([2.0, 3.0], [2.0, 2.0])
componentwise_meanvar(observations; bessel = false) # ([2.0, 3.0], [1.0, 1.0])

timestep_mean(sim, i)

Compute the ensemble mean at saved step index i.

For array-valued states, the returned value has the same shape as a single state and contains the componentwise mean across trajectories. For scalar states, the returned value is a scalar mean. Passing : computes the full step-indexed mean timeseries via timeseries_steps_mean.

timestep_median(sim, i)

Compute the ensemble median at saved step index i.

For array-valued states, the result is reshaped to match the state at sim.u[1].u[i]; for scalar states, it is the scalar median across trajectories. Passing : computes medians for every saved step.

timestep_quantile(sim, q, i)

Compute the componentwise quantile q at saved step index i.

q is passed to Statistics.quantile for each state component across trajectories. Array-valued states are reshaped to match the state at sim.u[1].u[i]; scalar states return a scalar quantile. Passing : computes the quantile for every saved step.

timestep_meanvar(sim, i)

Compute the ensemble mean and variance at saved step index i.

The result is (mean, variance), computed componentwise across trajectories with Bessel correction by the shared componentwise statistics helper. Passing : computes the full step-indexed mean and variance timeseries.

timestep_meancov(sim, i, j)

Compute componentwise means and covariance between saved step indices i and j.

The result is (mean_i, mean_j, covariance), where each entry is scalar-valued for scalar states or shaped componentwise for array-valued states. Passing (:, :) computes the full step-indexed covariance matrix.

timestep_meancor(sim, i, j)

Compute componentwise means and correlation between saved step indices i and j.

The result is (mean_i, mean_j, correlation), using the covariance and variance computed across trajectories. Passing (:, :) computes the full step-indexed correlation matrix.

timestep_weighted_meancov(sim, W, i, j)

Compute componentwise weighted means and covariance between saved step indices i and j.

W supplies the trajectory weights used by the weighted covariance calculation. The result is (mean_i, mean_j, weighted_covariance). Passing (:, :) computes the full step-indexed weighted covariance matrix.

timepoint_mean(sim, t)

Compute the ensemble mean at physical time t.

Each trajectory is evaluated with sol(t), so this requires a callable solution at t. For array-valued states, the result has the same shape as a single state; for scalar states, it is a scalar mean.

timepoint_median(sim, t)

Compute the componentwise ensemble median at physical time t.

Each trajectory is evaluated with sol(t). Array-valued states are reshaped to match a single saved state layout; scalar states return a scalar median.

timepoint_quantile(sim, q, t)

Compute the componentwise quantile q at physical time t.

Each trajectory is evaluated with sol(t), then Statistics.quantile is applied componentwise across trajectories.

timepoint_meanvar(sim, t)

Compute the ensemble mean and variance at physical time t.

The result is (mean, variance), computed componentwise across interpolated trajectory values at t.

timepoint_meancov(sim, t1, t2)

Compute componentwise means and covariance between physical times t1 and t2.

Each trajectory is evaluated at both times. The result is (mean_t1, mean_t2, covariance).

timepoint_meancor(sim, t1, t2)

Compute componentwise means and correlation between physical times t1 and t2.

Each trajectory is evaluated at both times. The result is (mean_t1, mean_t2, correlation).

timepoint_weighted_meancov(sim, W, t1, t2)

Compute componentwise weighted means and covariance between physical times t1 and t2.

W supplies the trajectory weights used by the weighted covariance calculation. The result is (mean_t1, mean_t2, weighted_covariance).

timeseries_steps_mean(sim)

Compute the ensemble mean at every saved step.

The result is a DiffEqArray with the same time vector as the first trajectory, where each saved value is the componentwise mean across trajectories at the same saved step index.

timeseries_steps_median(sim)

Compute the ensemble median at every saved step.

The result is a DiffEqArray with the first trajectory's time vector and componentwise median values at each saved step index.

timeseries_steps_quantile(sim, q)

Compute the componentwise quantile q at every saved step.

The result is a DiffEqArray with the first trajectory's time vector and componentwise quantile values at each saved step index.

timeseries_steps_meanvar(sim)

Compute the ensemble mean and variance at every saved step.

The result is (means, variances), where both entries are DiffEqArrays sharing the first trajectory's time vector.

timeseries_steps_meancov(sim)

Compute the step-indexed matrix of componentwise mean/covariance summaries.

Entry (i, j) contains the result of timestep_meancov. This assumes saved step indices are comparable across trajectories.

timeseries_steps_meancor(sim)

Compute the step-indexed matrix of componentwise mean/correlation summaries.

Entry (i, j) contains the result of timestep_meancor. This assumes saved step indices are comparable across trajectories.

timeseries_steps_weighted_meancov(sim, W)

Compute the step-indexed matrix of componentwise weighted covariance summaries.

Entry (i, j) contains the weighted mean/covariance summary for saved step indices i and j using trajectory weights W.

timeseries_point_mean(sim, ts)

Compute the ensemble mean at each physical time in ts.

The result is a DiffEqArray whose time axis is ts and whose values are the componentwise means of sol(t) across trajectories.

timeseries_point_median(sim, ts)

Compute the componentwise ensemble median at each physical time in ts.

The result is a DiffEqArray over ts; each value is computed from the interpolated trajectory values at that time.

timeseries_point_quantile(sim, q, ts)

Compute the componentwise quantile q at each physical time in ts.

The result is a DiffEqArray over ts; each value is computed from the interpolated trajectory values at that time.

timeseries_point_meanvar(sim, ts)

Compute the ensemble mean and variance at each physical time in ts.

The result is (means, variances), where both entries are DiffEqArrays over ts.

timeseries_point_meancov(sim, ts)

Compute the time-point covariance summary matrix for adjacent entries of ts.

This method pairs ts[1:end-1] with ts[2:end] and returns the same matrix form as timeseries_point_meancov(sim, ts1, ts2).

timeseries_point_meancov(sim, ts1, ts2)

Compute the time-point covariance summary matrix between two time collections.

Entry (i, j) contains the result of timepoint_meancov(sim, ts1[i], ts2[j]).

timeseries_point_meancor(sim, ts)

Compute the time-point correlation summary matrix for adjacent entries of ts.

This method pairs ts[1:end-1] with ts[2:end] and returns the same matrix form as timeseries_point_meancor(sim, ts1, ts2).

timeseries_point_meancor(sim, ts1, ts2)

Compute the time-point correlation summary matrix between two time collections.

Entry (i, j) contains the result of timepoint_meancor(sim, ts1[i], ts2[j]).

timeseries_point_weighted_meancov(sim, W, ts)

Compute the weighted covariance summary matrix for adjacent entries of ts.

This method pairs ts[1:end-1] with ts[2:end] and uses weights W for each trajectory.

timeseries_point_weighted_meancov(sim, W, ts1, ts2)

Compute the weighted covariance summary matrix between two time collections.

Entry (i, j) contains the weighted mean/covariance summary for ts1[i] and ts2[j] using trajectory weights W.

struct EnsembleSummary{T, N, Tt, S, S2, S3, S4, S5} <: SciMLBase.AbstractEnsembleSolution{T, N, S}

The EnsembleSummary type is included to help with analyzing the general summary statistics. Two constructors are provided:

EnsembleSummary(sim; quantiles = [0.05, 0.95])
EnsembleSummary(sim, ts; quantiles = [0.05, 0.95])

The first produces a (mean,var) summary at each time step. As with the summary statistics, this assumes that the time steps are all the same. The second produces a (mean,var) summary at each time point t in ts. This requires the ability to interpolate the solution. Quantile is used to determine the qlow and qhigh quantiles at each timepoint. It defaults to the 5% and 95% quantiles.

Complex-valued trajectories are rejected with a ComplexEnsembleSummaryError: medians and quantiles require a total order, which complex numbers do not have. Summarize real/imaginary parts or magnitudes separately, or use mean/variance helpers that do not need ordering.

Plot Recipe

The EnsembleSummary comes with a plot recipe for visualizing the summary statistics. The extra keyword arguments are:

• idxs: the solution components to plot. Defaults to plotting all components.
• error_style: The style for plotting the error. Defaults to ribbon. Other choices are :bars for error bars and :none for no error bars.
• ci_type : Defaults to :quantile which has (qlow,qhigh) quantiles whose limits were determined when constructing the EnsembleSummary. Gaussian CI 1.96*(standard error of the mean) can be set using ci_type=:SEM.

One useful argument is fillalpha which controls the transparency of the ribbon around the mean.