Bug fix in Hazard.get_mdr and impact calculation in case of a single exposure point and MDR values of 0

[bguillod]

Changes proposed in this PR:

• Fixes a bug in the impact calculation for cases with a single exposure point and a zero MDR value.

This PR fixes #933

PR Author Checklist

• Read the Contribution Guide
• Correct target branch selected (if unsure, select develop)
• Descriptive pull request title added
• Source branch up-to-date with target branch
• Documentation updated
• Tests updated
• Tests passing
• No new linter issues
• Changelog updated

PR Reviewer Checklist

• Read the Contribution Guide
• CLIMADA Reviewer Checklist passed
• Tests passing
• No new linter issues

[bguillod]

@emanuel-schmid The Test Pipelines fail due to no internet connection - can I safely assume this is an issues on the github configuration side and unrelated to changes in this PR?

[emanuel-schmid]

@emanuel-schmid The Test Pipelines fail due to no internet connection - can I safely assume this is an issues on the github configuration side and unrelated to changes in this PR?

@bguillod I do assume it's unrelated to changes in this PR, but I couldn't blame the github configuration either. Atm, I suspect the climada data server had a glitch.

[bguillod]

Ah, it works now and tests pass, so yes it must have been something.
I've assigned @peanutfun to this PR, but actually I guess he's off for a while now, correct? Who should we best assign to this PR @emanuel-schmid (beside @aleeciu who should check it indeed solves all the related issues for his cases)?

[chahank]

Thanks for the fix!

Can you please make a quick check that your modification does not lead to a significant loss of performance for large datasets?

[bguillod]

@chahank According to @peanutfun it should rather lead to a performance improvement, however we can check this if you want. Do you have a suggested case which uses large datasets which I could run? It also depends what you mean by "large" as I don't have an Euler cluster at hand...

[chahank]

Yes, as @peanutfun said, it might lead to an improvement too. I think it is enough if you run it for a case of a few thousand TC events and make a small uncertainty analysis (you can just execute the tutorial with one of your large datasets). There is no need to make extended tests, I would just like to have a quick check to be sure that we are not missing something obvious. Thanks!

[bguillod]

Alright, so based on the uncertainty analysis tutorial I created the following test:

#Define the input variable functions
import numpy as np

from climada.entity import ImpactFunc, ImpactFuncSet, Exposures
from climada.util.constants import EXP_DEMO_H5, HAZ_DEMO_H5
from climada.hazard import Hazard
from functools import partial
import scipy as sp
from climada.engine.unsequa import InputVar
from climada.engine.unsequa import CalcImpact
import time

def impf_func(G=1, v_half=84.7, vmin=25.7, k=3, _id=1):

    def xhi(v, v_half, vmin):
        return max([(v - vmin), 0]) / (v_half - vmin)

    def sigmoid_func(v, G, v_half, vmin, k):
        return G * xhi(v, v_half, vmin)**k / (1 + xhi(v, v_half, vmin)**k)

    #In-function imports needed only for parallel computing on Windows
    import numpy as np
    from climada.entity import ImpactFunc, ImpactFuncSet
    intensity_unit = 'm/s'
    intensity = np.linspace(0, 150, num=100)
    mdd = np.repeat(1, len(intensity))
    paa = np.array([sigmoid_func(v, G, v_half, vmin, k) for v in intensity])
    imp_fun = ImpactFunc("TC", _id, intensity, mdd, paa, intensity_unit)
    imp_fun.check()
    impf_set = ImpactFuncSet([imp_fun])
    return impf_set

haz = Hazard.from_hdf5(HAZ_DEMO_H5)
exp_base = Exposures.from_hdf5(EXP_DEMO_H5)

def prep_full_calc(haz, exp_base, impf_func):
    #It is a good idea to assign the centroids to the base exposures in order to avoid repeating this
    # potentially costly operation for each sample.
    exp_base.assign_centroids(haz)
    def exp_base_func(x_exp, exp_base):
        exp = exp_base.copy()
        exp.gdf['value'] *= x_exp
        return exp
    exp_func = partial(exp_base_func, exp_base=exp_base)
    #Define the InputVars


    exp_distr = {"x_exp": sp.stats.beta(10, 1.1)} #This is not really a reasonable distribution but is used
                                                #here to show that you can use any scipy distribution.

    exp_iv = InputVar(exp_func, exp_distr)

    impf_distr = {
        "G": sp.stats.truncnorm(0.5, 1.5),
        "v_half": sp.stats.uniform(35, 65),
        "vmin": sp.stats.uniform(0, 15),
        "k": sp.stats.uniform(1, 4)
        }
    impf_iv = InputVar(impf_func, impf_distr)
    calc_imp = CalcImpact(exp_iv, impf_iv, haz)
    output_imp = calc_imp.make_sample(N=2**7, sampling_kwargs={'skip_values': 2**8})
    return calc_imp, output_imp

def run_unc_calc(calc_imp, output_imp):
    calc_imp.uncertainty(output_imp, rp = [50, 100, 250])

I then did a few tests on the develop branch and on the branch used in this PR by first running:

calc_imp, output_imp = prep_full_calc(haz, exp_base, impf_func)

and then by profiling using this (everything here was run in a jupyter notebook):

%timeit -r 20 run_unc_calc(calc_imp, output_imp)

The results are not very stable (i.e. they change a little bit each time you run it), so sometimes there is a small difference and overall it tends to be slightly longer in this branch compared to develop. For instance, in the final run I did the outcome was:

1. In the develop branch: 3.19 s ± 13.7 ms per loop (mean ± std. dev. of 20 runs, 1 loop each).
2. In this branch: 3.2 s ± 21.6 ms per loop (mean ± std. dev. of 20 runs, 1 loop each).

Since it already takes several seconds to run this calculation I did not go on with a larger event set. This seems quite acceptable to me as it fixes a bug. Let me know what you think @chahank.

[chahank]

Ok, thanks for checking, this looks good then.

[chahank]

@bguillod did you want to do anything more? Otherwise I would merge. Thanks!

[aleeciu]

I checked the cases that allowed me to find the bug and this fix indeed solves the issue.. thanks @bguillod and @peanutfun

[bguillod]

Great, thanks @aleeciu for checking it resolves your case and @chahank for reviewing! All good from my perspective, so feel free to merge @chahank :-)