Distribution adapters and analysis cost#

Compare the cost of the same SORM analysis using PySTRA distributions and equivalent SciPy distribution adapters. Timings depend on the machine, dependency versions and solver settings; the measured ratios below are observations for this run.

Before you start: the introductory and SciPy-distribution tutorials. NumPy and SciPy are included in the core installation.

[1]:

import pystra as ra
from scipy.stats import norm, lognorm, uniform
import numpy as np
import timeit

Define a reasonably complex limit-state function

[2]:

def lsf(r, X1, X2, X3):
    g = r - X2 * (1000 * X3) ** (-1) - (X1 * (200 * X3) ** (-1)) ** 2
    return g

Define some objects common to both sets of analyses:

[3]:

limit_state = ra.LimitState(lsf)
options = ra.FORMOptions()

The analysis using the distributions built in to PySTRA:

[4]:

def run_builtin():
    """
    The basic example using built-in distributions
    """
    stochastic_model = ra.StochasticModel()
    stochastic_model.add_variable(ra.Lognormal("X1", 500, 100))
    stochastic_model.add_variable(ra.Normal("X2", 2000, 400))
    stochastic_model.add_variable(ra.Uniform("X3", 5, 0.5))
    stochastic_model.add_variable(ra.Constant("r", 1))

    stochastic_model.set_correlation(
        ra.CorrelationMatrix([[1.0, 0.3, 0.2], [0.3, 1.0, 0.2], [0.2, 0.2, 1.0]])
    )

    sorm = ra.SORM(
        model=stochastic_model,
        limit_state=limit_state,
        options=ra.SORMOptions(form=options),
    )
    sorm_result = sorm.run()
    assert sorm_result.converged
    return sorm_result.approximations["breitung"]

The analysis using the equivalent distributions in SciPy

[5]:

def run_scipy():
    """
    The basic example using Scipy distributions
    """
    stochastic_model = ra.StochasticModel()
    # Lognormal
    zeta = (np.log(1 + (100 / 500) ** 2)) ** 0.5
    lamb = np.log(500) - 0.5 * zeta**2
    stochastic_model.add_variable(
        ra.ScipyDist("X1", lognorm(s=zeta, scale=np.exp(lamb)))
    )
    # Normal
    stochastic_model.add_variable(ra.ScipyDist("X2", norm(loc=2000, scale=400)))
    ## Uniform
    a_b = (0.5**2 * 12) ** (1 / 2)
    a = (2 * 5 - a_b) / 2
    stochastic_model.add_variable(ra.ScipyDist("X3", uniform(loc=a, scale=a_b)))
    # Constant
    stochastic_model.add_variable(ra.Constant("r", 1))

    stochastic_model.set_correlation(
        ra.CorrelationMatrix([[1.0, 0.3, 0.2], [0.3, 1.0, 0.2], [0.2, 0.2, 1.0]])
    )

    sorm = ra.SORM(
        model=stochastic_model,
        limit_state=limit_state,
        options=ra.SORMOptions(form=options),
    )
    sorm_result = sorm.run()
    assert sorm_result.converged
    return sorm_result.approximations["breitung"]

Check agreement and measure execution time#

Check the two implementations before comparing their speed. This isolates distribution-adapter overhead; it does not predict the cost of an external structural model.

[6]:

builtin_probability = run_builtin()
scipy_probability = run_scipy()
np.testing.assert_allclose(builtin_probability, scipy_probability, rtol=1e-4)

number = 100
time_builtin = timeit.timeit(stmt=run_builtin, number=number)
time_scipy = timeit.timeit(stmt=run_scipy, number=number)

print("Total time taken (s)")
print(f"Built-in: {time_builtin}; Scipy: {time_scipy}")
print("Average time per call (s):")
print(f"Built-in: {time_builtin/number}; Scipy: {time_scipy/number}")
print(f"Built-in speed-up: {time_scipy/time_builtin:.2f}")

Total time taken (s)
Built-in: 1.4266015900000184; Scipy: 6.354021983999985
Average time per call (s):
Built-in: 0.014266015900000184; Scipy: 0.06354021983999986
Built-in speed-up: 4.45

Interpretation#

Use the numerical agreement check before drawing performance conclusions. For an expensive structural model, model-evaluation cost can dominate the distribution overhead measured here.

Continue: User guide · API reference · Theory