Pystra Timing

This notebook examines the speed up possible when using the built-in distribution function implementations, rather than the generic (but more comprehensive) ones provided by SciPy. Where performance is critical, implementation of a user-defined distribution object that avoids error checks and user input validation on each call could be beneficial.

A typical output is (the results will differ by machine):

```
Total time taken (s)
Built-in: 1.8919757150579244; Scipy: 5.2964311360847205
Average time per call (s):
Built-in: 0.018919757150579242; Scipy: 0.052964311360847206
Built-in speed-up: 2.80
```

Import the relevant packages (some may need to be installed as they are not part of the Pystra installation)

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

Define a reasonably complex limit state function

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 analyes:

limit_state = pr.LimitState(lsf)
options = pr.AnalysisOptions()
options.setPrintOutput(False)

The analysis using the distributions built in to Pystra:

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

    stochastic_model.setCorrelation(
        pr.CorrelationMatrix([[1.0, 0.3, 0.2], [0.3, 1.0, 0.2], [0.2, 0.2, 1.0]])
    )

    sorm = pr.Sorm(
        stochastic_model=stochastic_model,
        limit_state=limit_state,
        analysis_options=options
    )
    sorm.run()

The analysis using the equivalent distributions in SciPy

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

    stochastic_model.setCorrelation(
        pr.CorrelationMatrix([[1.0, 0.3, 0.2], [0.3, 1.0, 0.2], [0.2, 0.2, 1.0]])
    )

    sorm = pr.Sorm(
        stochastic_model=stochastic_model,
        limit_state=limit_state,
        analysis_options=options
    )
    sorm.run()

Now call and compare the timing results

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.7366494107991457; Scipy: 6.116523395990953
Average time per call (s):
Built-in: 0.017366494107991456; Scipy: 0.061165233959909526
Built-in speed-up: 3.52