stochops

A fluent Python framework for parallel Latin Hypercube Sampling (LHS) and uncertainty-quantification simulations on top of OpenSeesPy.

stochops runs large batches of structural analysis realizations in parallel, each with a different set of sampled input parameters (material strengths, section choices, ground motions...), and streams every result to disk — flat RAM regardless of campaign size.

How it works

A simulation campaign is assembled from four pluggable strategies, wired together by a fluent builder:

Strategy	Protocol	What it does
ParametricModel	callable(dict) -> dict	Runs a single realization of the FE model and returns Engineering Demand Parameters (EDPs) such as {"max_drift": 0.012}.
Sampler	sample(count) -> list[dict]	Generates parameter realizations. Implementations: LHSSampler (stratified LHS), MCSampler (independent draws).
ExecutionEngine	run_batch(model, samples) -> list[dict]	Dispatches batches across worker processes (ProcessPoolEngine, spawn-isolated for OpenSeesPy).
ResultSink	write(batch) / close()	Persists results. ParquetSink streams to disk; ArrowStreamSink buffers in memory for zero-disk handoff.

SimulationEngine.run(total_samples) loops in chunks: sample a batch → dispatch it → stream it to the sink, with an optional Rich progress bar (uv / pip style). A single failing realization is trapped and recorded with status="DIVERGED" / status="FAILED" metadata instead of killing the campaign.

Quick start

from stochops import SimulationBuilder
from stochops.execution.process_pool import ProcessPoolEngine
from stochops.parameters import ContinuousParam, DiscreteParam, Distribution
from stochops.samplers.lhs import LHSSampler
from stochops.samplers.mc import MCSampler
from stochops.sinks.parquet import ParquetSink

def evaluate_realization(params: dict) -> dict:
    """Single OpenSeesPy analysis. Here a cheap stand-in."""
    # ... build the OpenSees model from `params`, run the analysis ...
    return {
        "max_drift": (params["fy"] / params["fc"]) * 0.001,
        "base_shear_kN": (params["fc"] / 1e6) * 12.5,
    }

if __name__ == "__main__":
    # 1. Define the uncertain parameters
    parameters = [
        ContinuousParam("fc", Distribution.NORMAL, mean=30e6, std=3e6, min_value=20e6),
        ContinuousParam("fy", Distribution.LOGNORMAL, mean=400e6, std=20e6),
        DiscreteParam("section", choices=["RECT_300x500", "RECT_400x600"]),
    ]

    # 2. Assemble the simulation
    simulation = (
        SimulationBuilder()
        .model(evaluate_realization)
        # Use MCSampler instead for plain, independent Monte Carlo draws
        .sampler(LHSSampler(parameters=parameters, seed=2026))
        .execution_engine(ProcessPoolEngine(num_workers=8))
        .sink(ParquetSink("results.parquet"))
        .batch_size(100)
        .build()
    )

    # 3. Run 10,000 realizations in batches of 100
    simulation.run(total_samples=10_000)

Important: the model function must be defined at module level (top of the file). Workers run in spawn-isolated processes, which pickle the model by name — a function defined inside if __name__ == "__main__": or another function cannot be sent to the workers. The orchestration itself is wrapped in if __name__ == "__main__": so the workers' re-import of the script does not re-trigger the campaign.

Results can be read with anything that understands Parquet:

import pandas as pd
df = pd.read_parquet("results.parquet")

Columns include the input parameters, sample_id, status (CONVERGED / DIVERGED / FAILED), failure diagnostics, and every returned EDP.

In-memory handoff & progress

For analyses that keep results in memory — e.g. feeding directly into a training or post-processing step — use ArrowStreamSink instead of writing Parquet:

sink = ArrowStreamSink()
simulation = SimulationBuilder().sink(sink).batch_size(100).build()
simulation.run(total_samples=10_000)

import polars as pl          # or: sink.to_pandas()
df = pl.from_arrow(sink.to_arrow())   # zero-disk, no round-trip through a file

tensor = sink.to_torch()     # numeric columns only; pip install "stochops[torch]"

ArrowStreamSink implements the same ResultSink protocol as ParquetSink, so the two are drop-in swaps.

The campaign progress bar is on by default and can be configured through the builder:

SimulationBuilder().progress(show=True, style="uv")   # or "pip", or show=False

(uv is a minimal cyan/magenta theme, pip a classic block bar.)

Sampling notes

• Two strategies. LHSSampler draws stratified, space-filling Latin hypercube points; MCSampler draws independent, identically distributed realizations — the classic Monte Carlo default, with no stratification property but simpler statistical analysis. Both are drop-in swaps for the builder's .sampler(...) slot.
• Deterministic streams. Each sampler's RNG engine is created once and advanced across batches, so consecutive sample() calls draw distinct realizations. Two samplers built with the same seed reproduce the exact same concatenated stream; call sampler.reset() to restart.
• Correlated parameters. Pass a symmetric positive-definite correlation_matrix to impose correlations via a Gaussian copula. For LHSSampler the copula transformation reorders points and forfeits exact LHS stratification — use it only when correlation matters more than space-filling. For MCSampler, independent draws mean the copula has no such trade-off.
• Discrete parameters. DiscreteParam.probabilities are normalized automatically; GroundMotionSetParam draws uniformly from a catalog of ground-motion records.

Analysis robustness

Nonlinear transient analyses diverge. Two guards are built in:

• AdaptiveAnalysisRunner cycles through solution algorithms (Newton → KrylovNewton → NewtonWithLineSearch → BFGS → Broyden) and recursively halves the time step before giving up.
• The process pool isolates each worker with spawn, so OpenSeesPy C++ static state cannot leak across realizations and one segfaulting worker is replaced and retried rather than failing the campaign.

Limitations

• Nataf-correlated sampling (an alternative to the Gaussian copula) is planned but not yet implemented.
• The Parquet schema is frozen the first time the writer opens (after the first success row or max_buffered_rows records). Results introducing brand-new columns after that point emit a warning and drop those columns.
• Solver-divergence classification uses ConvergenceError type checks with a best-effort message fallback; unusual user exceptions may be mislabeled.

Development

uv sync --extra dev
uv run pytest
uv run ruff check src tests
uv run mypy src

Tested against Python 3.10+.