Source code for djura.vulnerability_modeller.fragility

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




def mlefit(param1: float, param2: float, total_count: int, count: int,
           data) -> float:
    """Maximum likelihood method
    Performs a lognormal cumulative distribution function fit to the data
    points based on maximum likelihood method

    Parameters
    ----------
    param1 : float
        Median of the function, parameter of a statistical model to be found
    param2 : float
        Standard deviation of the function, parameter of a statistical model
        to be found
    total_count : int
        Number of data points
    count : int
        Number of failures
    data: Union[List, np.ndarray]
        The function, data points

    Returns
    -------
    float
        Negative Log likelihood to be minimized
    """
    import warnings
    warnings.filterwarnings("ignore", category=RuntimeWarning)

    try:
        p = stats.norm.cdf(np.log(data), loc=np.log(param1), scale=param2)
        likelihood = stats.binom.pmf(count, total_count, p)
        likelihood[likelihood == 0] = 1e-290
        loglik = -sum(np.log10(likelihood))

        warnings.resetwarnings()

        return loglik
    except OverflowError:
        warnings.resetwarnings()
        return 1e+8
    except Exception:
        warnings.resetwarnings()
        return 1e+8





def fit_fragility(xs, counts, total_count, iml_range=None, beta=0.0):

    xs = xs
    counts = counts
    if xs[0] == 0:
        xs, counts = xs[1:], counts[1:]

    theta_hat_mom = np.exp(np.mean(np.log(xs)))
    beta_hat_mom = np.std(np.log(xs))

    x0 = [theta_hat_mom, beta_hat_mom]

    xopt_mle = optimize.fmin(
        func=lambda var: mlefit(
            param1=var[0], param2=var[1],
            total_count=total_count,
            count=np.array(counts), data=xs
        ),
        x0=x0, maxiter=3000, maxfun=3000, disp=False)

    theta_mle = xopt_mle[0]
    beta_mle = (xopt_mle[1] ** 2 + beta ** 2) ** 0.5

    probs = None
    if iml_range is not None:
        probs = stats.norm.cdf(
            np.log(iml_range / theta_mle) / beta_mle, loc=0, scale=1
        )

    return {
        'median': theta_mle,
        'beta': beta_mle,
        'probs': probs,
        'iml_range': iml_range
    }