Source code for djura.vulnerability_modeller.eal

# SPDX-License-Identifier: AGPL-3.0-or-later
# Copyright (C) 2025-2026 Djura | Risk - Data - Engineering S.r.l.
from typing import List
import numpy as np




class EAL:
    def __init__(self, elr: List[float], mafe: List[float]) -> None:
        """Object related to expected annual loss (EAL) computations

        Parameters
        ----------
        elr : List[float]
            Expected loss ratio (ELR)
        mafe : List[float]
            Mean annual frequency of exceedance (MAFE) of a limit state
        """
        self.elr = elr
        self.mafe = mafe



    def fit_loss_curve(self) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
        """
        Fitting of a refined loss curve that passes through the performance
        limit states

        Returns
        -------
        tuple[np.ndarray, np.ndarray, np.ndarray]
            Coefficients of the analytical form
            Fitted ELRs
            Fitted MAFEs

        Notes
        -----
        The analytical form is defined as:

        .. math::
            MAFE = \\mathrm{coef}_0 * exp(
                -\\mathrm{coef}_1 * log(ELR)
                - \\mathrm{coef}_2 * log(ELR)^2)

        # noqa: E501, W605
        """
        if len(self.elr) != len(self.mafe) != 3:
            raise ValueError(
                "Length of ELR and MAFE should be equal to 3 for fitting the "
                "refined curve")

        coef = np.zeros(len(self.elr))

        r1 = np.zeros(len(self.elr))
        r2 = np.zeros(len(self.elr))
        r3 = np.zeros(len(self.elr))
        r1[0] = 1
        r2[0] = 1
        r3[0] = 1
        r1[1] = -np.log(self.elr[0])
        r2[1] = -np.log(self.elr[1])
        r3[1] = -np.log(self.elr[2])
        r1[2] = -np.log(self.elr[0]) ** 2
        r2[2] = -np.log(self.elr[1]) ** 2
        r3[2] = -np.log(self.elr[2]) ** 2

        temp1 = np.log(self.mafe)
        temp2 = np.array([r1, r2, r3])
        temp3 = np.linalg.inv(temp2).dot(temp1)
        temp3 = temp3.tolist()

        coef[0] = np.exp(temp3[0])
        coef[1] = temp3[1]
        coef[2] = temp3[2]

        elr_fit = np.linspace(0.01, 1., 100)
        mafe_fit = coef[0] * np.exp(
            -coef[1] * np.log(elr_fit) - coef[2] * np.log(elr_fit) ** 2)

        elr_fit = np.insert(elr_fit, 0, 0.0)
        mafe_fit = np.insert(mafe_fit, 0, mafe_fit[0])

        return coef, elr_fit, mafe_fit




    def get_eal(
        self, im_range, elr: List[float] = None, mafe: List[float] = None
    ) -> float:
        """Computes approximate EAL using the ELRs and MAFEs associated with
        the performance limit states

        Parameters
        -------
        elr : List[float], optional
            Expected loss ratio (ELR), by default EAL.elr
        mafe : List[float], optional
            Mean annual frequency of exceedance (MAFE) of a limit state,
            by default EAL.mafe

        Returns
        -------
        float
            EAL
        """
        if mafe is None:
            mafe = np.array(self.mafe)
        if elr is None:
            elr = np.array(self.elr)

        # Hazard IML tests
        diml = np.diff(im_range)

        # dMAFEdIML, logarithmic gradient divided by IML step
        dmafe_diml = np.log(mafe[1:] / mafe[:-1]) / diml

        # Loss ratio step
        dmdf = np.diff(elr)

        # EAL contributions of each subgroup
        eal_bins = np.zeros(dmafe_diml.shape)
        eal_bins[dmafe_diml != 0] = elr[:-1] * mafe[:-1] * (
            1 - np.exp(dmafe_diml * diml)) - dmdf / diml * mafe[:-1] * (
            np.exp(dmafe_diml * diml) *
            (diml - 1 / dmafe_diml) + 1 / dmafe_diml)

        eal = sum(eal_bins) + mafe[-1]

        # eal = mafe[0] * elr[0]
        # eals = (mafe[:-1] + mafe[1:]) / 2 * np.diff(elr)

        # return eal + sum(eals)
        return eal