Fire Spread

Rothermel Fire Spread Model — Demange-Chryst et al. (2022), Example 4.3¶

This notebook replicates the fire spread reliability-oriented sensitivity analysis from Demange-Chryst et al. [1]. The model is the semi-physical Rothermel model [2] for wildland surface fire spread, with modifications by Albini [3] and Catchpole & Catchpole [4], as adapted by Song et al. [5].

Model. Rate of fire spread $R$ (cm $\cdot$ s$^{-1}$) as a function of 10 fuel, weather, and terrain parameters. The failure threshold is $R > t = 60$ cm $\cdot$ s$^{-1}$.

Input variables (10) with mixed Normal / LogNormal marginals, several truncated, and one correlation: $\rho(m_d, U) = -0.8$.

Reference failure probability: $p_f \approx 1.4 \times 10^{-4}$.

The reference target Shapley effects from the paper [1, Table 4] are used for validation.

[1] Demange-Chryst, J., Bachoc, F., & Morio, J. (2022). IJUQ, arXiv:2202.12679. [2] Rothermel, R.C. (1972). A mathematical model for predicting fire spread in wildland fuels. USDA Forest Service. [3] Albini, F.A. (1976). Estimating wildfire behavior and effects. USDA Forest Service. [4] Catchpole, E.A. & Catchpole, W.R. (1991). Int. J. Wildland Fire, 1:101–106. [5] Song, E., Nelson, B.L., & Staum, J. (2016). SIAM/ASA JUQ, 4(1):1060–1083.

In [1]:

Copied!





# Import dependencies
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.stats import norm, lognorm

from shapleyx.utilities.mc_shapley import (
    MultivariateNormal,
    shapley_effects,
)

from importlib.metadata import version
print(f"Running on ShapleyX v{version('shapleyx')}")
# Import dependencies
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.stats import norm, lognorm

from shapleyx.utilities.mc_shapley import (
    MultivariateNormal,
    shapley_effects,
)

from importlib.metadata import version
print(f"Running on ShapleyX v{version('shapleyx')}")

Running on ShapleyX v0.5.1

Input Distribution¶

Ten variables — four LogNormal, six Normal (several truncated), with one correlation pair: $\rho(m_d, U) = -0.8$.

#	Symbol	Description	Distribution
1	$\delta$	Fuel depth (cm)	LogN(2.19, 0.517)
2	$\sigma$	Fuel particle area-to-volume ratio (cm$^{-1}$)	LogN(3.31, 0.294), truncated $\ge 3/0.6$
3	$h$	Fuel particle low heat content (Kcal$\cdot$kg$^{-1}$)	LogN(8.48, 0.063)
4	$\rho_p$	Oven-dry particle density (D.W. g$\cdot$cm$^{-3}$)	LogN($-$0.592, 0.219)
5	$m_l$	Moisture content of live fuel (H$_2$O g D.W. g$^{-1}$)	N(1.18, 0.377), $\ge 0$
6	$m_d$	Moisture content of dead fuel (H$_2$O g D.W. g$^{-1}$)	N(0.19, 0.047)
7	$S_T$	Fuel particle total mineral content (MIN. g D.W. g$^{-1}$)	N(0.049, 0.011), $\ge 0$
8	$U$	Wind speed at midflame height (km$\cdot$h$^{-1}$)	6.9 $\times$ LogN(1.0174, 0.5569)
9	$\tan\phi$	Slope	N(0.38, 0.186), $\ge 0$
10	$P$	Dead fuel loading / total fuel loading	LogN($-$2.19, 0.64), $\le 1$

Correlation: $\rho(m_d, U) = -0.8$ (windier $\rightarrow$ less moisture).

Rejection rules: All negative values rejected; $S_T \le 1$; $P \le 1$; $\sigma \ge 3/0.6$.

In [2]:

Copied!





d = 10

labels = ['delta', 'sigma', 'h', 'rho_p', 'm_l', 'm_d', 'S_T', 'U', 'tan_phi', 'P']
descriptions = [
    'Fuel depth (cm)',
    'Fuel particle area-to-volume ratio (cm⁻¹)',
    'Fuel particle low heat content (Kcal·kg⁻¹)',
    'Oven-dry particle density (D.W. g·cm⁻³)',
    'Moisture content of live fuel (H₂O g D.W. g⁻¹)',
    'Moisture content of dead fuel (H₂O g D.W. g⁻¹)',
    'Fuel particle total mineral content (MIN. g D.W. g⁻¹)',
    'Wind speed at midflame height (km·h⁻¹)',
    'Slope (tan φ)',
    'Dead fuel loading / total fuel loading',
]

# Marginal parameters: (dist_type, mu, sigma) for the UNDERLYING normal
# For LogNormal: mu and sigma are for log(X) ~ N(mu, sigma)
# For Normal: mu and sigma are mean and std
marginals = {
    'delta':   ('lognormal', 2.19, 0.517),
    'sigma':   ('lognormal', 3.31, 0.294),
    'h':       ('lognormal', 8.48, 0.063),
    'rho_p':   ('lognormal', -0.592, 0.219),
    'm_l':     ('normal',    1.18, 0.377),
    'm_d':     ('normal',    0.19, 0.047),
    'S_T':     ('normal',    0.049, 0.011),
    'U':       ('lognormal_scaled', 1.0174, 0.5569, 6.9),
    'tan_phi': ('normal',    0.38, 0.186),
    'P':       ('lognormal', -2.19, 0.64),
}

# Correlation matrix (Pearson) — only md and U are correlated
# Indices: 0=delta, 1=sigma, 2=h, 3=rho_p, 4=m_l, 5=m_d, 6=S_T, 7=U, 8=tan_phi, 9=P
corr_pearson = np.eye(d)
corr_pearson[5, 7] = corr_pearson[7, 5] = -0.8  # ρ(m_d, U) = -0.8

print("Pearson correlation matrix (non-zero off-diagonals):")
for i in range(d):
    for j in range(i+1, d):
        if abs(corr_pearson[i, j]) > 1e-10:
            print(f"  ρ({labels[i]}, {labels[j]}) = {corr_pearson[i, j]:.1f}")

# Summary table
print(f"\n{'#':>2s}  {'Var':>8s}  {'Type':>16s}  {'Params':>30s}  {'Description'}")
print("-" * 95)
for i, name in enumerate(labels):
    dist_type = marginals[name][0]
    if dist_type == 'lognormal':
        params = f"LogN(mu={marginals[name][1]:.3f}, sigma={marginals[name][2]:.3f})"
    elif dist_type == 'lognormal_scaled':
        params = f"{marginals[name][3]:.1f} × LogN(mu={marginals[name][1]:.4f}, sigma={marginals[name][2]:.4f})"
    else:
        params = f"N(mu={marginals[name][1]:.3f}, sigma={marginals[name][2]:.3f})"
    print(f"{i+1:>2d}  {name:>8s}  {dist_type:>16s}  {params:>30s}  {descriptions[i]}")
d = 10

labels = ['delta', 'sigma', 'h', 'rho_p', 'm_l', 'm_d', 'S_T', 'U', 'tan_phi', 'P']
descriptions = [
    'Fuel depth (cm)',
    'Fuel particle area-to-volume ratio (cm⁻¹)',
    'Fuel particle low heat content (Kcal·kg⁻¹)',
    'Oven-dry particle density (D.W. g·cm⁻³)',
    'Moisture content of live fuel (H₂O g D.W. g⁻¹)',
    'Moisture content of dead fuel (H₂O g D.W. g⁻¹)',
    'Fuel particle total mineral content (MIN. g D.W. g⁻¹)',
    'Wind speed at midflame height (km·h⁻¹)',
    'Slope (tan φ)',
    'Dead fuel loading / total fuel loading',
]

# Marginal parameters: (dist_type, mu, sigma) for the UNDERLYING normal
# For LogNormal: mu and sigma are for log(X) ~ N(mu, sigma)
# For Normal: mu and sigma are mean and std
marginals = {
    'delta':   ('lognormal', 2.19, 0.517),
    'sigma':   ('lognormal', 3.31, 0.294),
    'h':       ('lognormal', 8.48, 0.063),
    'rho_p':   ('lognormal', -0.592, 0.219),
    'm_l':     ('normal',    1.18, 0.377),
    'm_d':     ('normal',    0.19, 0.047),
    'S_T':     ('normal',    0.049, 0.011),
    'U':       ('lognormal_scaled', 1.0174, 0.5569, 6.9),
    'tan_phi': ('normal',    0.38, 0.186),
    'P':       ('lognormal', -2.19, 0.64),
}

# Correlation matrix (Pearson) — only md and U are correlated
# Indices: 0=delta, 1=sigma, 2=h, 3=rho_p, 4=m_l, 5=m_d, 6=S_T, 7=U, 8=tan_phi, 9=P
corr_pearson = np.eye(d)
corr_pearson[5, 7] = corr_pearson[7, 5] = -0.8  # ρ(m_d, U) = -0.8

print("Pearson correlation matrix (non-zero off-diagonals):")
for i in range(d):
    for j in range(i+1, d):
        if abs(corr_pearson[i, j]) > 1e-10:
            print(f"  ρ({labels[i]}, {labels[j]}) = {corr_pearson[i, j]:.1f}")

# Summary table
print(f"\n{'#':>2s}  {'Var':>8s}  {'Type':>16s}  {'Params':>30s}  {'Description'}")
print("-" * 95)
for i, name in enumerate(labels):
    dist_type = marginals[name][0]
    if dist_type == 'lognormal':
        params = f"LogN(mu={marginals[name][1]:.3f}, sigma={marginals[name][2]:.3f})"
    elif dist_type == 'lognormal_scaled':
        params = f"{marginals[name][3]:.1f} × LogN(mu={marginals[name][1]:.4f}, sigma={marginals[name][2]:.4f})"
    else:
        params = f"N(mu={marginals[name][1]:.3f}, sigma={marginals[name][2]:.3f})"
    print(f"{i+1:>2d}  {name:>8s}  {dist_type:>16s}  {params:>30s}  {descriptions[i]}")

Pearson correlation matrix (non-zero off-diagonals):
  ρ(m_d, U) = -0.8

 #       Var              Type                          Params  Description
-----------------------------------------------------------------------------------------------
 1     delta         lognormal     LogN(mu=2.190, sigma=0.517)  Fuel depth (cm)
 2     sigma         lognormal     LogN(mu=3.310, sigma=0.294)  Fuel particle area-to-volume ratio (cm⁻¹)
 3         h         lognormal     LogN(mu=8.480, sigma=0.063)  Fuel particle low heat content (Kcal·kg⁻¹)
 4     rho_p         lognormal    LogN(mu=-0.592, sigma=0.219)  Oven-dry particle density (D.W. g·cm⁻³)
 5       m_l            normal        N(mu=1.180, sigma=0.377)  Moisture content of live fuel (H₂O g D.W. g⁻¹)
 6       m_d            normal        N(mu=0.190, sigma=0.047)  Moisture content of dead fuel (H₂O g D.W. g⁻¹)
 7       S_T            normal        N(mu=0.049, sigma=0.011)  Fuel particle total mineral content (MIN. g D.W. g⁻¹)
 8         U  lognormal_scaled  6.9 × LogN(mu=1.0174, sigma=0.5569)  Wind speed at midflame height (km·h⁻¹)
 9   tan_phi            normal        N(mu=0.380, sigma=0.186)  Slope (tan φ)
10         P         lognormal    LogN(mu=-2.190, sigma=0.640)  Dead fuel loading / total fuel loading

Gaussian Copula Distribution Class¶

Handles mixed LogNormal + Normal + scaled LogNormal marginals with a latent Gaussian copula and truncation/rejection rules. Conditional sampling in latent normal space, mapped back through inverse CDFs.

In [3]:

Copied!





class GaussianCopulaFire:
    """Gaussian copula for the fire spread model's mixed distributions.

    Supports LogNormal, Normal, and scaled LogNormal marginals
    with truncation/rejection rules matching Demange-Chryst et al.
    Conditional sampling is performed in latent normal space.
    """
    def __init__(self, marginals, latent_corr):
        self.d = len(marginals)
        self.labels = list(marginals.keys())
        self._marginals = marginals
        self._latent_corr = np.asarray(latent_corr)

        # Underlying multivariate normal for latent-space sampling
        self._mvn = MultivariateNormal(
            mean=np.zeros(self.d), cov=self._latent_corr
        )

    @staticmethod
    def _marginal_to_latent(x, dist_type, *params):
        """Transform original -> N(0,1)."""
        if dist_type == 'lognormal':
            mu, sigma = params
            x = np.clip(np.asarray(x), 1e-15, None)
            return (np.log(x) - mu) / sigma
        elif dist_type == 'lognormal_scaled':
            mu, sigma, scale = params
            x = np.clip(np.asarray(x), 1e-15, None)
            return (np.log(x / scale) - mu) / sigma
        else:  # normal
            mu, sigma = params
            return (np.asarray(x) - mu) / sigma

    @staticmethod
    def _marginal_from_latent(z, dist_type, *params):
        """Transform N(0,1) -> original."""
        if dist_type == 'lognormal':
            mu, sigma = params
            return np.exp(mu + sigma * np.asarray(z))
        elif dist_type == 'lognormal_scaled':
            mu, sigma, scale = params
            return scale * np.exp(mu + sigma * np.asarray(z))
        else:  # normal
            mu, sigma = params
            return mu + sigma * np.asarray(z)

    def _to_original(self, Z):
        """Transform latent normal -> original space."""
        X = np.zeros_like(Z)
        for j in range(self.d):
            name = self.labels[j]
            dist_type, *params = self._marginals[name]
            X[:, j] = self._marginal_from_latent(Z[:, j], dist_type, *params)
        return X

    def _apply_rejections(self, X):
        """Apply truncation/rejection rules. Returns valid-only rows."""
        mask = np.ones(len(X), dtype=bool)
        # All negative values rejected
        mask &= np.all(X > 0, axis=1)
        # S_T <= 1 (variable index 6)
        mask &= X[:, 6] <= 1.0
        # P <= 1 (variable index 9)
        mask &= X[:, 9] <= 1.0
        # sigma >= 3/0.6 (variable index 1)
        mask &= X[:, 1] >= 3.0 / 0.6
        return X[mask]

    def sample_joint(self, n, max_rejections=20):
        """Draw n joint samples, applying rejections."""
        # Oversample and reject to account for truncation
        X_all = []
        remaining = n
        for _ in range(max_rejections):
            Z = self._mvn.sample_joint(int(remaining * 1.2 + 100))
            X_candidate = self._to_original(Z)
            X_valid = self._apply_rejections(X_candidate)
            X_all.append(X_valid)
            remaining = n - sum(len(x) for x in X_all)
            if remaining <= 0:
                break
        return np.vstack(X_all)[:n]

    def sample_conditional_batch(self, u_indices, fixed_X):
        """Draw N conditional samples (one per row of fixed_X)."""
        u = np.asarray(u_indices)
        N = fixed_X.shape[0]
        fixed_X = np.asarray(fixed_X, dtype=float)

        if len(u) == 0:
            return self.sample_joint(N)

        # Map fixed variables to latent space
        Z_fixed = np.zeros((N, len(u)))
        for k, idx in enumerate(u):
            name = self.labels[idx]
            dist_type, *params = self._marginals[name]
            Z_fixed[:, k] = self._marginal_to_latent(
                fixed_X[:, k], dist_type, *params
            )

        # Condition in latent normal space
        Z_cond = self._mvn.sample_conditional_batch(u, Z_fixed)

        # Map back to original space
        return self._to_original(Z_cond)


joint = GaussianCopulaFire(marginals, corr_pearson)
print("GaussianCopulaFire distribution ready.")
print(f"Dimension: {joint.d}")
class GaussianCopulaFire:
    """Gaussian copula for the fire spread model's mixed distributions.

    Supports LogNormal, Normal, and scaled LogNormal marginals
    with truncation/rejection rules matching Demange-Chryst et al.
    Conditional sampling is performed in latent normal space.
    """
    def __init__(self, marginals, latent_corr):
        self.d = len(marginals)
        self.labels = list(marginals.keys())
        self._marginals = marginals
        self._latent_corr = np.asarray(latent_corr)

        # Underlying multivariate normal for latent-space sampling
        self._mvn = MultivariateNormal(
            mean=np.zeros(self.d), cov=self._latent_corr
        )

    @staticmethod
    def _marginal_to_latent(x, dist_type, *params):
        """Transform original -> N(0,1)."""
        if dist_type == 'lognormal':
            mu, sigma = params
            x = np.clip(np.asarray(x), 1e-15, None)
            return (np.log(x) - mu) / sigma
        elif dist_type == 'lognormal_scaled':
            mu, sigma, scale = params
            x = np.clip(np.asarray(x), 1e-15, None)
            return (np.log(x / scale) - mu) / sigma
        else:  # normal
            mu, sigma = params
            return (np.asarray(x) - mu) / sigma

    @staticmethod
    def _marginal_from_latent(z, dist_type, *params):
        """Transform N(0,1) -> original."""
        if dist_type == 'lognormal':
            mu, sigma = params
            return np.exp(mu + sigma * np.asarray(z))
        elif dist_type == 'lognormal_scaled':
            mu, sigma, scale = params
            return scale * np.exp(mu + sigma * np.asarray(z))
        else:  # normal
            mu, sigma = params
            return mu + sigma * np.asarray(z)

    def _to_original(self, Z):
        """Transform latent normal -> original space."""
        X = np.zeros_like(Z)
        for j in range(self.d):
            name = self.labels[j]
            dist_type, *params = self._marginals[name]
            X[:, j] = self._marginal_from_latent(Z[:, j], dist_type, *params)
        return X

    def _apply_rejections(self, X):
        """Apply truncation/rejection rules. Returns valid-only rows."""
        mask = np.ones(len(X), dtype=bool)
        # All negative values rejected
        mask &= np.all(X > 0, axis=1)
        # S_T <= 1 (variable index 6)
        mask &= X[:, 6] <= 1.0
        # P <= 1 (variable index 9)
        mask &= X[:, 9] <= 1.0
        # sigma >= 3/0.6 (variable index 1)
        mask &= X[:, 1] >= 3.0 / 0.6
        return X[mask]

    def sample_joint(self, n, max_rejections=20):
        """Draw n joint samples, applying rejections."""
        # Oversample and reject to account for truncation
        X_all = []
        remaining = n
        for _ in range(max_rejections):
            Z = self._mvn.sample_joint(int(remaining * 1.2 + 100))
            X_candidate = self._to_original(Z)
            X_valid = self._apply_rejections(X_candidate)
            X_all.append(X_valid)
            remaining = n - sum(len(x) for x in X_all)
            if remaining <= 0:
                break
        return np.vstack(X_all)[:n]

    def sample_conditional_batch(self, u_indices, fixed_X):
        """Draw N conditional samples (one per row of fixed_X)."""
        u = np.asarray(u_indices)
        N = fixed_X.shape[0]
        fixed_X = np.asarray(fixed_X, dtype=float)

        if len(u) == 0:
            return self.sample_joint(N)

        # Map fixed variables to latent space
        Z_fixed = np.zeros((N, len(u)))
        for k, idx in enumerate(u):
            name = self.labels[idx]
            dist_type, *params = self._marginals[name]
            Z_fixed[:, k] = self._marginal_to_latent(
                fixed_X[:, k], dist_type, *params
            )

        # Condition in latent normal space
        Z_cond = self._mvn.sample_conditional_batch(u, Z_fixed)

        # Map back to original space
        return self._to_original(Z_cond)


joint = GaussianCopulaFire(marginals, corr_pearson)
print("GaussianCopulaFire distribution ready.")
print(f"Dimension: {joint.d}")

GaussianCopulaFire distribution ready.
Dimension: 10

In [4]:

Copied!





# Validate: check acceptance rate and basic distribution properties
X_check = joint.sample_joint(100000)
print(f"Acceptance rate: {len(X_check) / (len(X_check) * 1.0):.2f}")
print(f"\n{'Variable':>8s}  {'Mean':>12s}  {'Std':>12s}  {'Min':>12s}  {'Max':>12s}")
print("-" * 56)
for i, name in enumerate(labels):
    print(f"{name:>8s}  {X_check[:, i].mean():12.4g}  {X_check[:, i].std():12.4g}"
          f"  {X_check[:, i].min():12.4g}  {X_check[:, i].max():12.4g}")

# Check empirical correlation between m_d and U
emp_corr = np.corrcoef(X_check.T)
print(f"\nEmpirical ρ(m_d, U) = {emp_corr[5, 7]:.4f}  (target: -0.8)")
# Validate: check acceptance rate and basic distribution properties
X_check = joint.sample_joint(100000)
print(f"Acceptance rate: {len(X_check) / (len(X_check) * 1.0):.2f}")
print(f"\n{'Variable':>8s}  {'Mean':>12s}  {'Std':>12s}  {'Min':>12s}  {'Max':>12s}")
print("-" * 56)
for i, name in enumerate(labels):
    print(f"{name:>8s}  {X_check[:, i].mean():12.4g}  {X_check[:, i].std():12.4g}"
          f"  {X_check[:, i].min():12.4g}  {X_check[:, i].max():12.4g}")

# Check empirical correlation between m_d and U
emp_corr = np.corrcoef(X_check.T)
print(f"\nEmpirical ρ(m_d, U) = {emp_corr[5, 7]:.4f}  (target: -0.8)")

Acceptance rate: 1.00

Variable          Mean           Std           Min           Max
--------------------------------------------------------
   delta         10.22         5.695        0.7174         71.27
   sigma         28.61         8.625          7.55         96.56
       h          4827         302.5          3718          6217
   rho_p        0.5669        0.1257        0.2232         1.472
     m_l          1.18        0.3771      0.000147         2.796
     m_d        0.1902       0.04694      0.002498        0.4283
     S_T       0.04904       0.01097      0.001476        0.0974
       U         22.25         13.54         1.528         235.4
 tan_phi        0.3892        0.1758     0.0001029         1.142
       P        0.1368       0.09634      0.005104        0.9947

Empirical ρ(m_d, U) = -0.7366  (target: -0.8)

Rothermel Fire Spread Model¶

The rate of spread $R$ (cm $\cdot$ s$^{-1}$) computed from Rothermel's equations with Albini and Catchpole modifications. Internal unit conversions from metric to imperial follow the Demange-Chryst et al. implementation exactly.

In [5]:

Copied!





# --- Unit conversion helpers ---
def kg2lb(x):
    return 2.204622 * x

def m2ft(x):
    return 3.28083 * x

def ft2m(x):
    return 0.3047995 * x


def fire_spread_rate(x):
    """Rate of fire spread (cm/s) from Rothermel's equations.

    x = [delta, sigma, h, rho_p, m_l, m_d, S_T, U, tan_phi, P]

    Returns 0 for invalid (rejected) inputs.
    """
    delta, sigma, h, rho_p, m_l, m_d, S_T, U, tan_phi, P = x

    # Rejection rules
    if delta <= 0 or h <= 0 or rho_p <= 0 or U <= 0 or P <= 0:
        return 0.0
    if S_T >= 1:
        return 0.0

    # Fuel loading w_0 (kg/m²)
    w_0 = 4.8 / 4.8824 * 1.0 / (1.0 + np.exp((15.0 - delta) / 3.5))
    # Convert to lb/ft²
    w_0 = kg2lb(w_0) * 0.09290272

    # Convert metric inputs to imperial
    delta_ft = m2ft(delta / 100.0)
    sigma_ft = 1.0 / m2ft(0.01 / sigma)
    rho_p_lb = kg2lb(rho_p) * 28.3167
    h_btu = 4184.0 * h / 1055.87 * 0.4235924
    U_ftmin = m2ft(U * 1000.0) / 60.0

    # Rothermel equations
    Gamma_max = sigma_ft**(3.0/2.0) / (495.0 + 0.0594 * sigma_ft**(3.0/2.0))
    beta_op = 3.348 * sigma_ft**(-0.8189)
    A_val = 133.0 * sigma_ft**(-0.7913)
    theta_star = (301.4 - 305.87 * (m_l - m_d) + 2260.0 * m_d) / m_l / 2260.0
    theta = min(1.0, max(theta_star, 0.0))
    mu_m = np.exp(-7.3 * P * m_d - (7.3 * theta + 2.13) * (1.0 - P) * m_l)
    mu_S = min(0.174 * S_T**(-0.19), 1.0)
    C = 7.47 * np.exp(-0.133 * sigma_ft**(0.55))
    B_val = 0.02526 * sigma_ft**(0.54)
    E_val = 0.715 * np.exp(-3.59e-4 * sigma_ft)
    w_n = w_0 * (1.0 - S_T)
    rho_b = w_0 / delta_ft
    epsilon = np.exp(-138.0 / sigma_ft)
    Q_ig = 130.87 + 1054.43 * m_d
    beta = rho_b / rho_p_lb
    Gamma = (Gamma_max * (beta / beta_op)**A_val
             * np.exp(A_val * (1.0 - beta / beta_op)))

    in_exp = (0.792 + 0.681 * sigma_ft**(0.5)) * (beta + 0.1)
    in_exp = min(in_exp, 709.78)

    ksi = np.exp(in_exp) / (192.0 + 0.2595 * sigma_ft)
    phi_W = C * U_ftmin**B_val * (beta / beta_op)**(-E_val)
    phi_S = 5.275 * beta**(-0.3) * tan_phi**2
    I_R = Gamma * w_n * h_btu * mu_m * mu_S

    # Rate of spread in ft/min -> cm/s
    R_ft_min = (I_R * ksi * (1.0 + phi_W + phi_S)) / (rho_b * epsilon * Q_ig)
    return ft2m(R_ft_min) * 100.0 / 60.0


def fire_spread_batch(X):
    """Vectorised fire spread rate for batch evaluation."""
    delta = X[:, 0]; sigma = X[:, 1]; h = X[:, 2]
    rho_p = X[:, 3]; m_l = X[:, 4]; m_d = X[:, 5]
    S_T = X[:, 6]; U = X[:, 7]; tan_phi = X[:, 8]; P = X[:, 9]

    # Rejection mask
    mask = ((delta > 0) & (h > 0) & (rho_p > 0) & (U > 0) & (P > 0)
            & (S_T < 1))
    R = np.zeros(len(X))

    dm, sm, hm, rm, mm_l, mm_d, sm_T, Um, tm, Pm = [x[mask] for x in
        [delta, sigma, h, rho_p, m_l, m_d, S_T, U, tan_phi, P]]

    # Fuel loading
    w_0 = 4.8 / 4.8824 / (1.0 + np.exp((15.0 - dm) / 3.5))
    w_0 = kg2lb(w_0) * 0.09290272

    # Convert to imperial
    dm_ft = m2ft(dm / 100.0)
    sm_ft = 1.0 / m2ft(0.01 / sm)
    rm_lb = kg2lb(rm) * 28.3167
    hm_btu = 4184.0 * hm / 1055.87 * 0.4235924
    Um_ftmin = m2ft(Um * 1000.0) / 60.0

    # Rothermel equations
    Gamma_max = sm_ft**1.5 / (495.0 + 0.0594 * sm_ft**1.5)
    beta_op = 3.348 * sm_ft**(-0.8189)
    A_val = 133.0 * sm_ft**(-0.7913)
    theta_star = (301.4 - 305.87 * (mm_l - mm_d) + 2260.0 * mm_d) / mm_l / 2260.0
    theta = np.clip(theta_star, 0.0, 1.0)
    mu_m = np.exp(-7.3 * Pm * mm_d - (7.3 * theta + 2.13) * (1.0 - Pm) * mm_l)
    mu_S = np.minimum(0.174 * sm_T**(-0.19), 1.0)
    C = 7.47 * np.exp(-0.133 * sm_ft**(0.55))
    B_val = 0.02526 * sm_ft**(0.54)
    E_val = 0.715 * np.exp(-3.59e-4 * sm_ft)
    w_n = w_0 * (1.0 - sm_T)
    rho_b = w_0 / dm_ft
    epsilon = np.exp(-138.0 / sm_ft)
    Q_ig = 130.87 + 1054.43 * mm_d
    beta = rho_b / rm_lb
    Gamma = (Gamma_max * (beta / beta_op)**A_val
             * np.exp(A_val * (1.0 - beta / beta_op)))

    in_exp = np.minimum((0.792 + 0.681 * sm_ft**0.5) * (beta + 0.1), 709.78)
    ksi = np.exp(in_exp) / (192.0 + 0.2595 * sm_ft)
    phi_W = C * Um_ftmin**B_val * (beta / beta_op)**(-E_val)
    phi_S = 5.275 * beta**(-0.3) * tm**2
    I_R = Gamma * w_n * hm_btu * mu_m * mu_S

    R_ft_min = (I_R * ksi * (1.0 + phi_W + phi_S)) / (rho_b * epsilon * Q_ig)
    R[mask] = ft2m(R_ft_min) * 100.0 / 60.0
    return R


t = 60.0  # failure threshold (cm/s)
print(f"Failure threshold: R > {t} cm/s")
# --- Unit conversion helpers ---
def kg2lb(x):
    return 2.204622 * x

def m2ft(x):
    return 3.28083 * x

def ft2m(x):
    return 0.3047995 * x


def fire_spread_rate(x):
    """Rate of fire spread (cm/s) from Rothermel's equations.

    x = [delta, sigma, h, rho_p, m_l, m_d, S_T, U, tan_phi, P]

    Returns 0 for invalid (rejected) inputs.
    """
    delta, sigma, h, rho_p, m_l, m_d, S_T, U, tan_phi, P = x

    # Rejection rules
    if delta <= 0 or h <= 0 or rho_p <= 0 or U <= 0 or P <= 0:
        return 0.0
    if S_T >= 1:
        return 0.0

    # Fuel loading w_0 (kg/m²)
    w_0 = 4.8 / 4.8824 * 1.0 / (1.0 + np.exp((15.0 - delta) / 3.5))
    # Convert to lb/ft²
    w_0 = kg2lb(w_0) * 0.09290272

    # Convert metric inputs to imperial
    delta_ft = m2ft(delta / 100.0)
    sigma_ft = 1.0 / m2ft(0.01 / sigma)
    rho_p_lb = kg2lb(rho_p) * 28.3167
    h_btu = 4184.0 * h / 1055.87 * 0.4235924
    U_ftmin = m2ft(U * 1000.0) / 60.0

    # Rothermel equations
    Gamma_max = sigma_ft**(3.0/2.0) / (495.0 + 0.0594 * sigma_ft**(3.0/2.0))
    beta_op = 3.348 * sigma_ft**(-0.8189)
    A_val = 133.0 * sigma_ft**(-0.7913)
    theta_star = (301.4 - 305.87 * (m_l - m_d) + 2260.0 * m_d) / m_l / 2260.0
    theta = min(1.0, max(theta_star, 0.0))
    mu_m = np.exp(-7.3 * P * m_d - (7.3 * theta + 2.13) * (1.0 - P) * m_l)
    mu_S = min(0.174 * S_T**(-0.19), 1.0)
    C = 7.47 * np.exp(-0.133 * sigma_ft**(0.55))
    B_val = 0.02526 * sigma_ft**(0.54)
    E_val = 0.715 * np.exp(-3.59e-4 * sigma_ft)
    w_n = w_0 * (1.0 - S_T)
    rho_b = w_0 / delta_ft
    epsilon = np.exp(-138.0 / sigma_ft)
    Q_ig = 130.87 + 1054.43 * m_d
    beta = rho_b / rho_p_lb
    Gamma = (Gamma_max * (beta / beta_op)**A_val
             * np.exp(A_val * (1.0 - beta / beta_op)))

    in_exp = (0.792 + 0.681 * sigma_ft**(0.5)) * (beta + 0.1)
    in_exp = min(in_exp, 709.78)

    ksi = np.exp(in_exp) / (192.0 + 0.2595 * sigma_ft)
    phi_W = C * U_ftmin**B_val * (beta / beta_op)**(-E_val)
    phi_S = 5.275 * beta**(-0.3) * tan_phi**2
    I_R = Gamma * w_n * h_btu * mu_m * mu_S

    # Rate of spread in ft/min -> cm/s
    R_ft_min = (I_R * ksi * (1.0 + phi_W + phi_S)) / (rho_b * epsilon * Q_ig)
    return ft2m(R_ft_min) * 100.0 / 60.0


def fire_spread_batch(X):
    """Vectorised fire spread rate for batch evaluation."""
    delta = X[:, 0]; sigma = X[:, 1]; h = X[:, 2]
    rho_p = X[:, 3]; m_l = X[:, 4]; m_d = X[:, 5]
    S_T = X[:, 6]; U = X[:, 7]; tan_phi = X[:, 8]; P = X[:, 9]

    # Rejection mask
    mask = ((delta > 0) & (h > 0) & (rho_p > 0) & (U > 0) & (P > 0)
            & (S_T < 1))
    R = np.zeros(len(X))

    dm, sm, hm, rm, mm_l, mm_d, sm_T, Um, tm, Pm = [x[mask] for x in
        [delta, sigma, h, rho_p, m_l, m_d, S_T, U, tan_phi, P]]

    # Fuel loading
    w_0 = 4.8 / 4.8824 / (1.0 + np.exp((15.0 - dm) / 3.5))
    w_0 = kg2lb(w_0) * 0.09290272

    # Convert to imperial
    dm_ft = m2ft(dm / 100.0)
    sm_ft = 1.0 / m2ft(0.01 / sm)
    rm_lb = kg2lb(rm) * 28.3167
    hm_btu = 4184.0 * hm / 1055.87 * 0.4235924
    Um_ftmin = m2ft(Um * 1000.0) / 60.0

    # Rothermel equations
    Gamma_max = sm_ft**1.5 / (495.0 + 0.0594 * sm_ft**1.5)
    beta_op = 3.348 * sm_ft**(-0.8189)
    A_val = 133.0 * sm_ft**(-0.7913)
    theta_star = (301.4 - 305.87 * (mm_l - mm_d) + 2260.0 * mm_d) / mm_l / 2260.0
    theta = np.clip(theta_star, 0.0, 1.0)
    mu_m = np.exp(-7.3 * Pm * mm_d - (7.3 * theta + 2.13) * (1.0 - Pm) * mm_l)
    mu_S = np.minimum(0.174 * sm_T**(-0.19), 1.0)
    C = 7.47 * np.exp(-0.133 * sm_ft**(0.55))
    B_val = 0.02526 * sm_ft**(0.54)
    E_val = 0.715 * np.exp(-3.59e-4 * sm_ft)
    w_n = w_0 * (1.0 - sm_T)
    rho_b = w_0 / dm_ft
    epsilon = np.exp(-138.0 / sm_ft)
    Q_ig = 130.87 + 1054.43 * mm_d
    beta = rho_b / rm_lb
    Gamma = (Gamma_max * (beta / beta_op)**A_val
             * np.exp(A_val * (1.0 - beta / beta_op)))

    in_exp = np.minimum((0.792 + 0.681 * sm_ft**0.5) * (beta + 0.1), 709.78)
    ksi = np.exp(in_exp) / (192.0 + 0.2595 * sm_ft)
    phi_W = C * Um_ftmin**B_val * (beta / beta_op)**(-E_val)
    phi_S = 5.275 * beta**(-0.3) * tm**2
    I_R = Gamma * w_n * hm_btu * mu_m * mu_S

    R_ft_min = (I_R * ksi * (1.0 + phi_W + phi_S)) / (rho_b * epsilon * Q_ig)
    R[mask] = ft2m(R_ft_min) * 100.0 / 60.0
    return R


t = 60.0  # failure threshold (cm/s)
print(f"Failure threshold: R > {t} cm/s")

Failure threshold: R > 60.0 cm/s

In [24]:

Copied!





# Quick histogram of rate of spread
R_sample = fire_spread_batch(X_check[:100000])
R_sample_valid = R_sample[R_sample > 0]

fig, ax = plt.subplots(figsize=(8, 4))
ax.hist(R_sample_valid, bins=100, density=True, color='darkorange', alpha=0.85)
ax.axvline(t, color='crimson', linestyle='--', linewidth=2, label=f'Threshold t = {t} cm/s')
ax.set_xlabel('Rate of Spread R (cm/s)')
ax.set_ylabel('Density')
ax.set_title('Fire Spread Rate Distribution (Rothermel Model)')
ax.legend()
ax.set_xlim(0, min(200, R_sample_valid.max()))
plt.tight_layout()
plt.show()

pf_est = np.mean(R_sample > t)
print(f"Estimated failure probability: {pf_est:.2e}")
print(f"  (paper: ~1.4e-4 — with N=100000 this estimate is noisy)")
print(f"Expected failures in 100000 samples: {100000 * pf_est:.0f}")
# Quick histogram of rate of spread
R_sample = fire_spread_batch(X_check[:100000])
R_sample_valid = R_sample[R_sample > 0]

fig, ax = plt.subplots(figsize=(8, 4))
ax.hist(R_sample_valid, bins=100, density=True, color='darkorange', alpha=0.85)
ax.axvline(t, color='crimson', linestyle='--', linewidth=2, label=f'Threshold t = {t} cm/s')
ax.set_xlabel('Rate of Spread R (cm/s)')
ax.set_ylabel('Density')
ax.set_title('Fire Spread Rate Distribution (Rothermel Model)')
ax.legend()
ax.set_xlim(0, min(200, R_sample_valid.max()))
plt.tight_layout()
plt.show()

pf_est = np.mean(R_sample > t)
print(f"Estimated failure probability: {pf_est:.2e}")
print(f"  (paper: ~1.4e-4 — with N=100000 this estimate is noisy)")
print(f"Expected failures in 100000 samples: {100000 * pf_est:.0f}")

No description has been provided for this image

Estimated failure probability: 1.60e-04
  (paper: ~1.4e-4 — with N=100000 this estimate is noisy)
Expected failures in 100000 samples: 16

Target Shapley Effects (on Failure Indicator $\mathbf{1}_{R>60}$)¶

The paper computes target (reliability-oriented) Shapley effects on the binary failure indicator. With $p_f \approx 1.4 \times 10^{-4}$, standard MC is very inefficient — most samples are non-failures. The paper uses importance sampling to address this.

Here we attempt standard MC with a large $N$ to see what's achievable without IS, and compare against the paper's reference values (computed with IS and $N = 10^7$, $N_O = 10^6$, $N_I = 3$).

In [8]:

Copied!





# Estimate pf more accurately
X_pf = joint.sample_joint(500000)
pf_est2 = np.mean(fire_spread_batch(X_pf) > t)
print(f"Estimated failure probability (500k samples): {pf_est2:.2e}")
print(f"  Paper: ~1.4e-4")
print(f"Expected failures in N=20000: {20000 * pf_est2:.0f}")
# Estimate pf more accurately
X_pf = joint.sample_joint(500000)
pf_est2 = np.mean(fire_spread_batch(X_pf) > t)
print(f"Estimated failure probability (500k samples): {pf_est2:.2e}")
print(f"  Paper: ~1.4e-4")
print(f"Expected failures in N=20000: {20000 * pf_est2:.0f}")

Estimated failure probability (500k samples): 1.28e-04
  Paper: ~1.4e-4
Expected failures in N=20000: 3

In [13]:

Copied!





def failure_indicator(x):
    """Binary failure: 1 if rate of spread exceeds threshold, else 0."""
    R = fire_spread_rate(x)
    return 1.0 if R > t else 0.0


def failure_indicator_batch(X):
    """Vectorised failure indicator."""
    R = fire_spread_batch(X)
    return (R > t).astype(float)


# Target Shapley effects (standard MC — limited by pf)
# With pf ~ 1.4e-4, only ~14 failures per 100,000 samples.
# This is a stress test of standard MC for rare events.
eff_t, sh_t, var_t, lower_t, upper_t = shapley_effects(
    failure_indicator,
    joint,
    N=300000,
    method='exhaustive',
    B=200,
    alpha=0.05,
    predict_batch=failure_indicator_batch,
    random_state=42,
    progress=True,
)

results_target = pd.DataFrame({
    'Variable': labels,
    'Target Shapley': eff_t,
    'Lower': lower_t,
    'Upper': upper_t,
})
results_target
def failure_indicator(x):
    """Binary failure: 1 if rate of spread exceeds threshold, else 0."""
    R = fire_spread_rate(x)
    return 1.0 if R > t else 0.0


def failure_indicator_batch(X):
    """Vectorised failure indicator."""
    R = fire_spread_batch(X)
    return (R > t).astype(float)


# Target Shapley effects (standard MC — limited by pf)
# With pf ~ 1.4e-4, only ~14 failures per 100,000 samples.
# This is a stress test of standard MC for rare events.
eff_t, sh_t, var_t, lower_t, upper_t = shapley_effects(
    failure_indicator,
    joint,
    N=300000,
    method='exhaustive',
    B=200,
    alpha=0.05,
    predict_batch=failure_indicator_batch,
    random_state=42,
    progress=True,
)

results_target = pd.DataFrame({
    'Variable': labels,
    'Target Shapley': eff_t,
    'Lower': lower_t,
    'Upper': upper_t,
})
results_target

MC Shapley: 100%|██████████| 613500000/613500000 [16:22<00:00, 624333.48evals/s]

Out[13]:

	Variable	Target Shapley	Lower	Upper
0	delta	0.157281	0.138779	0.184766
1	sigma	0.201809	0.176075	0.245215
2	h	0.023252	-0.019919	0.047672
3	rho_p	0.023132	-0.018595	0.049574
4	m_l	0.143506	0.124052	0.172470
5	m_d	0.134637	0.113433	0.158789
6	S_T	0.030478	-0.012900	0.056654
7	U	0.167145	0.145114	0.197286
8	tan_phi	0.041205	0.011414	0.065695
9	P	0.077556	0.054687	0.097649

Comparison with Demange-Chryst Reference Values¶

The paper's Table 4 provides reference target Shapley effects computed with importance sampling ($N = 10^7$, $N_O = 10^6$, $N_I = 3$). We compare our standard-MC estimates against these values.

Note: With $p_f \approx 1.4 \times 10^{-4}$, standard MC target Shapley estimates are expected to have high variance. The Demange-Chryst paper notes that existing non-IS estimators "return a value very close to 0 almost every time because there are too few failure points." This comparison illustrates exactly why importance sampling is essential for ROSA with rare events.

In [16]:

Copied!





# Reference target Shapley effects from the paper (Table 4)
ref_target_shap = np.array([0.152, 0.247, 0.011, 0.003, 0.162,
                            0.145, 0.016, 0.182, 0.009, 0.073])

comparison = pd.DataFrame({
    'Variable': labels,
    'Target Shapley (MC)': eff_t.round(4),
    'Paper target (IS ref.)': ref_target_shap.round(3),
})
comparison
# Reference target Shapley effects from the paper (Table 4)
ref_target_shap = np.array([0.152, 0.247, 0.011, 0.003, 0.162,
                            0.145, 0.016, 0.182, 0.009, 0.073])

comparison = pd.DataFrame({
    'Variable': labels,
    'Target Shapley (MC)': eff_t.round(4),
    'Paper target (IS ref.)': ref_target_shap.round(3),
})
comparison

Out[16]:

	Variable	Target Shapley (MC)	Paper target (IS ref.)
0	delta	0.1573	0.152
1	sigma	0.2018	0.247
2	h	0.0233	0.011
3	rho_p	0.0231	0.003
4	m_l	0.1435	0.162
5	m_d	0.1346	0.145
6	S_T	0.0305	0.016
7	U	0.1671	0.182
8	tan_phi	0.0412	0.009
9	P	0.0776	0.073

In [20]:

Copied!





# Bar chart comparison
fig, ax = plt.subplots(figsize=(16, 7))
bar_width = 0.22
x = np.arange(d)

ax.bar(x- bar_width/2, eff_t, bar_width,
       yerr=[eff_t - lower_t, upper_t - eff_t],
       capsize=2, color='darkorange', alpha=0.85,
       label='Target Shapley (MC, no IS)')
ax.bar(x + bar_width/2, ref_target_shap, bar_width,
       color='crimson', alpha=0.85,
       label='Paper target (IS ref.)')

ax.set_xlabel('Input Variable')
ax.set_ylabel('Sensitivity Index')
ax.set_title('Fire Spread Model — Shapley Effects\\n'
             'Variance-Based vs. Target (with/without Importance Sampling)')
ax.set_xticks(x)
ax.set_xticklabels(labels)
ax.legend(fontsize=9)
ax.set_ylim(0, None)
plt.tight_layout()
plt.show()
# Bar chart comparison
fig, ax = plt.subplots(figsize=(16, 7))
bar_width = 0.22
x = np.arange(d)

ax.bar(x- bar_width/2, eff_t, bar_width,
       yerr=[eff_t - lower_t, upper_t - eff_t],
       capsize=2, color='darkorange', alpha=0.85,
       label='Target Shapley (MC, no IS)')
ax.bar(x + bar_width/2, ref_target_shap, bar_width,
       color='crimson', alpha=0.85,
       label='Paper target (IS ref.)')

ax.set_xlabel('Input Variable')
ax.set_ylabel('Sensitivity Index')
ax.set_title('Fire Spread Model — Shapley Effects\\n'
             'Variance-Based vs. Target (with/without Importance Sampling)')
ax.set_xticks(x)
ax.set_xticklabels(labels)
ax.legend(fontsize=9)
ax.set_ylim(0, None)
plt.tight_layout()
plt.show()

Key Takeaways¶

Target Shapley effects without importance sampling are challenging to estimate when $p_f \approx 10^{-4}$. The Demange-Chryst paper's IS-based approach is essential for reliable estimation in this regime.
The Gaussian copula with mixed marginals successfully handles the complex distribution specification: LogNormal, scaled LogNormal, Normal, truncated distributions, and a strong correlation ($\rho = -0.8$) between $m_d$ and $U$.
This notebook demonstrates the limit of standard MC for ROSA with rare events — a key motivation for the importance sampling estimators developed by Demange-Chryst et al. ShapleyX does not currently implement IS, making this a natural area for future development.
Rejection sampling in the distribution class handles the truncation rules, though with some efficiency loss for strongly truncated marginals.
The 5 most influential inputs for target Shapley effects are the same ($\delta$, $\sigma$, $m_l$, $m_d$, $U$), as noted in the paper.

In [ ]:

Copied!

%load_ext watermark
%watermark -n -u -v -iv -w
%load_ext watermark
%watermark -n -u -v -iv -w