Slower MonteCarlo simulations for versions > 1.19

[seahawks8]

What happened?

Times should be equal or better than previous versions, not slower.
It does not matter how much I run it, it's consistently slower on the newer OpenTurns versions.

Note that I had to change my code to an example I can post, in production the slowdown was much higher than the time reported in this example, closer to 3x slower

Python 3.8 on Windows 11 64bit

• Version 1.22 = Average runtime of 7.8 seconds
• Version 1.19 = Average runtime of 6.3 seconds.

How to reproduce the issue?

Steps
1. install openturns 1.19 for python 3.8 and run the code below, record the average time
2. install openturns 1.22 for python 3.8 and run the code below, record the average time


import time
import math
import openturns as ot


def sample_func(X):
    x0, x1, x2, x3, x4, x5, x6, x7, x8 = X
    y0 = x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8
    y1 = (x1 - 1.0) * math.exp(x0) * x2 + y0
    return [y1]


if __name__ == "__main__":
    dimension = 9
    cor_mat = ot.CorrelationMatrix(dimension)
    cor_mat[4, 5] = cor_mat[5, 4] = 0.931657
    copula = ot.NormalCopula(cor_mat)

    new_marginals = [
        # wt 0
        ot.Normal(0.413386, 0.00971457),
        # od 1
        ot.Normal(48, 0.0816),
        # d 2
        ot.TruncatedDistribution(ot.Normal(0.0644013, 0.0177285), 0, 0.200035),
        # LG 3
        ot.LogNormal(-3.2626, 1.2359, 0.0004),
        # ys
        ot.TruncatedDistribution(
            ot.Normal(56064, 3137), 52000, ot.TruncatedDistribution.LOWER
        ),
        # TS
        ot.Normal(79528, 5123),
        # cv
        ot.LogNormal(2.99379, 0.571324, 3.2),
        # # model error
        eval(
            "ot.Histogram(0.991865082, [0.10401330133343734, 0.10401330133343734, 0.10401330133343734, 0.10401330133343734, 0.10401330133343734, 0.10401330133343734], [2.0240326293902124, 1.5180244720426592, 4.048065258780425, 0.5060081573475531, 0.5060081573475531, 1.0120163146951062])"
        ),
        # # pressure NOP
        ot.TruncatedDistribution(
            ot.Mixture(
                [
                    ot.Normal(273.914, 28.9517),
                    ot.Normal(167.913, 20.2474),
                    ot.Normal(93.1261, 27.211),
                    ot.Normal(224.136, 22.8987),
                ],
                [0.2169, 0.3903, 0.0748, 0.318],
            ),
            1.15,
            659.4762,
        ),
    ]

    composed_dist = ot.ComposedDistribution(new_marginals, copula)
    rv_vect = ot.RandomVector(composed_dist)
    function = ot.PythonFunction(9, 1, sample_func)

    composed_distribution = ot.ComposedDistribution(new_marginals, copula)
    composed_distribution.setName("Distributions")

    rv_vector = ot.RandomVector(composed_distribution)

    thresh_dir = ot.Less()
    ls_vect = ot.CompositeRandomVector(function, rv_vector)

    # ThresholdEvent
    event = ot.ThresholdEvent(ls_vect, thresh_dir, 0.0)

    # MonteCarloExperiment
    experiment = ot.MonteCarloExperiment()

    # ProbabilitySimulationAlgorithm
    algo = ot.ProbabilitySimulationAlgorithm(event, experiment)

    sims = 50000
    blocks = 20

    algo.setMaximumOuterSampling(sims)
    algo.setBlockSize(blocks)

    algo.setMaximumCoefficientOfVariation(0.000001)
    algo.setName("MC")

    all_times = []
    for i in range(10):
        start = time.time()
        algo.run()
        result = algo.getResult()
        end = time.time()
        print(
            f"{ot.__version__}:{blocks}|{sims}|{result.getProbabilityEstimate()}|{result.getCoefficientOfVariation()}|{end - start}"
        )

        all_times.append(end - start)

    print(f"Average time: {sum(all_times)/len(all_times)} for {ot.__version__}")

### Version

Multiple versions, 1.19 and 1.22

### Operating System

Windows

### Installation media

pip

### Additional Context

_No response_

[jschueller]

I get better results with this example if I set a lower number of cores than the default number (40 HT cores) with ot.TBB.SetThreadsNumber(10) , maybe thats a lead ? what does ot.TBB.GetThreadsNumber() give on your machine ? Maybe its different across OT releases as TBB is upgraded

[seahawks8]

Here are my results, which are rather interesting but consistent with my findings on OT 1.19 is much faster (even in this case)

OS	Python	Threads	OT Version 1.22 Avg Time (Seconds)	OT Version 1.19 Avg Time (Seconds)
Windows 11	3.8	2	6.63	5.01
Windows 11	3.8	4	5.44	3.78
Windows 11	3.8	8	5.46	3.82
Windows 11	3.8	9	5.53	4.02
Windows 11	3.8	10	5.57	4.38
Windows 11	3.8	11	5.42	4.53
Windows 11	3.8	12	5.95	4.84
Windows 11	3.8	13	5.79	5.0
Windows 11	3.8	15	5.87	5.35
Windows 11	3.8	20	6.33	5.52

Updated script to loop thread numbers

import time
import math
import openturns as ot


def sample_func(X):
    x0, x1, x2, x3, x4, x5, x6, x7, x8 = X
    y0 = x1 + x2 + x3 + x4 + x5 + x6 + x7 + x8
    y1 = (x1 - 1.0) * math.exp(x0) * x2 + y0
    return [y1]


if __name__ == "__main__":
    dimension = 9
    cor_mat = ot.CorrelationMatrix(dimension)
    cor_mat[4, 5] = cor_mat[5, 4] = 0.931657
    copula = ot.NormalCopula(cor_mat)

    new_marginals = [
        # wt 0
        ot.Normal(0.413386, 0.00971457),
        # od 1
        ot.Normal(48, 0.0816),
        # d 2
        ot.TruncatedDistribution(ot.Normal(0.0644013, 0.0177285), 0, 0.200035),
        # LG 3
        ot.LogNormal(-3.2626, 1.2359, 0.0004),
        # ys
        ot.TruncatedDistribution(
            ot.Normal(56064, 3137), 52000, ot.TruncatedDistribution.LOWER
        ),
        # TS
        ot.Normal(79528, 5123),
        # cv
        ot.LogNormal(2.99379, 0.571324, 3.2),
        # # model error
        eval(
            "ot.Histogram(0.991865082, [0.10401330133343734, 0.10401330133343734, 0.10401330133343734, 0.10401330133343734, 0.10401330133343734, 0.10401330133343734], [2.0240326293902124, 1.5180244720426592, 4.048065258780425, 0.5060081573475531, 0.5060081573475531, 1.0120163146951062])"
        ),
        # # pressure NOP
        ot.TruncatedDistribution(
            ot.Mixture(
                [
                    ot.Normal(273.914, 28.9517),
                    ot.Normal(167.913, 20.2474),
                    ot.Normal(93.1261, 27.211),
                    ot.Normal(224.136, 22.8987),
                ],
                [0.2169, 0.3903, 0.0748, 0.318],
            ),
            1.15,
            659.4762,
        ),
    ]

    composed_dist = ot.ComposedDistribution(new_marginals, copula)
    rv_vect = ot.RandomVector(composed_dist)
    function = ot.PythonFunction(9, 1, sample_func)

    composed_distribution = ot.ComposedDistribution(new_marginals, copula)
    composed_distribution.setName("Distributions")

    rv_vector = ot.RandomVector(composed_distribution)

    thresh_dir = ot.Less()
    ls_vect = ot.CompositeRandomVector(function, rv_vector)

    # ThresholdEvent
    event = ot.ThresholdEvent(ls_vect, thresh_dir, 0.0)

    # MonteCarloExperiment
    experiment = ot.MonteCarloExperiment()

    # ProbabilitySimulationAlgorithm
    algo = ot.ProbabilitySimulationAlgorithm(event, experiment)

    sims = 50000
    blocks = 20

    algo.setMaximumOuterSampling(sims)
    algo.setBlockSize(blocks)

    algo.setMaximumCoefficientOfVariation(0.000001)
    algo.setName("MC")

    print(ot.__version__)
    for b in [2, 4, 8, 9, 10, 11, 12, 13, 15, 20]:
        ot.TBB.SetThreadsNumber(b)
        threads = ot.TBB.GetThreadsNumber()
        all_times = []
        for i in range(5):
            start = time.time()
            algo.run()
            result = algo.getResult()
            end = time.time()
            all_times.append(end - start)

        print(
            f"| Windows 11 | 3.8 | {threads} | {round(sum(all_times)/len(all_times),2)} | |"
        )

[jschueller]

614b5d3 broke eval by blocks (>=1.22)

[jschueller]

thanks for the detailed report, will be fixed for 1.23