oq-lite: investigate block size dependency in event based risk

[Michele Simionato]

By changing the number of concurrent tasks one gets different numbers. For instance, try to run oq-risklib/openquake/qa_tests_data/event_based_risk/case_2/ with -c=0 and -c=7.

Result of the investigation: this is not a bug of oq-lite. Instead it is a bug of the vulnerability functions which affects also the engine when the
coefficients of correlation are nonzero. Having fixed that bug, there is
no block size dependence in oq-lite.

[Michele Simionato (michele-simionato)]

https://github.com/gem/oq-risklib/pull/263

[Michele Simionato (michele-simionato)]

https://github.com/gem/oq-engine/pull/1763

[Michele Simionato (michele-simionato)]

The bug potentially affects any calculation with vulnerability
functions with nonzero coefficients of variation. By fixing the bug
and the tests, the numbers produced by the event based calculator and
the scenario calculator change. Usually the change is insignificant
(< 1%) but in some test the numbers are substantially different, see
for instance the change in
openquake/qa_tests_data/event_based_risk/case_3/expected/mean-structural-agg_loss_curve.csv

Notice that while the engine does not exhibit block size dependency (the
reason is that the parallelization algorithm splits by asset and not by
ruptures and the problem is visible only when splitting by rupture, as
oq-lite does) it is still affected by the bug. This means that fixing
the bug will change slightly the numbers obtained by our users, even if they
are using just the engine and not oq-lite.

Having said that, let me explain what the bug is.

The `LogNormalDistribution` class takes in input a matrix N x R of epsilons,
where N is the number of assets and R the number of ruptures. It has a
`sample(means, covs, stddevs)` method which takes three arrays of size R'
and returns an array of probabilities of size R'. What is R'? It is
the number of ground motion values greater than the minimum value
accepted by the vulnerability function. Just to be concrete I will use
the numbers from the test where I discovered the issue. It is a test
where the vulnerability function has the following intensity measure levels:

 imls = [ 0.02, 0.3 , 0.5 , 0.9 , 1.2 ]

there is a single asset (i.e. N=1) and there are 20 ruptures (i.e. R=20)
giving the following 20 ground motion values:

gmvs = [ 0.3307648 , 0.77900947, 0. , 1.2 , 0. ,
        0. , 0.42448847, 0. , 0. , 0. ,
        0.15023323, 0. , 0. , 0. , 0. ,
        0. , 0. , 0. , 0.51451394, 0. ])

The minimum intensity measure level is 0.2, and there are only R' = 6 ground
motions greater or equal than 0.2:

[0.3307648, 0.77900947, 1.2, 0.42448847, 0.15023323, 0.51451394]

The positions of such values with respect to the original array,
counting from zero, are

idxs = [0, 1, 3, 6, 10, 18]

Such positions are extremely important, because EACH GROUND MOTION VALUE
IS ASSOCIATED TO AN EPSILON. In order to be correct, if we discard the
ground motion values smaller than the minimum, we also need to discard
the corresponding epsilons; conversely, if we keep the ground motion
values at positions [0, 1, 3, 6, 10, 18] we must keep the epsilons
at the SAME POSITIONS. Instead the code was just taking the first 6
epsilons, without regard to the positions. This was at the origin
of the block size issue: by splitting the ground motion fields in different
blocks, the epsilons were extracted differently depending on the block splitting.
The fix is to pass the array of the indices `idxs` to the `sample` method
of the `LogNormalDistribution`, so that the right epsilons are taken,
independently from the way the blocks are split.