Perturbation-space heatmaps

This notebook illustrates how to use demovuln outputs to visualize population reduction across slices of the perturbation space.

The example uses a simple three-stage model with one juvenile stage, one subadult stage, and one adult stage. The model is intentionally minimal, but the same workflow applies to empirical projection matrices.

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

from demovuln import MatrixPopulationModel, PerturbationGrid, run_grid

Define a three-stage projection matrix

The model contains one fecundity row and explicit stage definitions. Here, the adult source stage is the third column.

A = np.array([
    [0.0, 0.0, 1.2],
    [0.35, 0.55, 0.0],
    [0.0, 0.20, 0.85],
])

model = MatrixPopulationModel(
    A,
    adult_stages=[2],
    juvenile_stages=[0, 1],
    name="three_stage_example",
)

print(f"Dominant eigenvalue: {model.lambda_:.3f}")
print("Stable stage distribution:", model.stable_distribution())

Helper function for heatmaps

The table returned by run_grid is transformed into a two-dimensional matrix using pandas.pivot_table. The plotting function below keeps each heatmap in its own figure.

def plot_heatmap(table, *, x, y, value="population_reduction", title=None):
    pivot = table.pivot_table(index=y, columns=x, values=value, aggfunc="mean")

    fig, ax = plt.subplots(figsize=(7, 4.5))
    image = ax.imshow(
        pivot.values,
        origin="lower",
        aspect="auto",
        extent=[
            pivot.columns.min(),
            pivot.columns.max(),
            pivot.index.min(),
            pivot.index.max(),
        ],
    )
    ax.set_xlabel(x)
    ax.set_ylabel(y)
    if title is not None:
        ax.set_title(title)
    cbar = fig.colorbar(image, ax=ax)
    cbar.set_label("Population reduction (%)")
    fig.tight_layout()
    plt.show()

Magnitude-duration slices at fixed period

We first fix the perturbation period and vary magnitude and duration. This is analogous to taking a two-dimensional slice through the full perturbation space.

fixed_period_grid = PerturbationGrid(
    magnitudes=np.linspace(0, 0.8, 17),
    durations=range(0, 7),
    periods=[8],
)

fixed_period_results = {}

for target in ["adult_survival", "juvenile_survival", "fecundity"]:
    fixed_period_results[target] = run_grid(
        model,
        target=target,
        grid=fixed_period_grid,
        t_max=48,
        recovery_steps=12,
    )
    print(f"{target}: Φ = {fixed_period_results[target].vulnerability:.2f}%")

for target, result in fixed_period_results.items():
    plot_heatmap(
        result.table,
        x="magnitude",
        y="duration",
        title=f"{target}: magnitude-duration slice at period = 8",
    )

Duration-period slices at fixed magnitude

We can also fix perturbation magnitude and examine how population reduction depends on duration and recurrence.

fixed_magnitude = 0.25
duration_period_grid = PerturbationGrid(
    magnitudes=[fixed_magnitude],
    durations=range(0, 9),
    periods=range(1, 13),
)

adult_duration_period = run_grid(
    model,
    target="adult_survival",
    grid=duration_period_grid,
    t_max=60,
    recovery_steps=12,
    skip_infeasible=True,
)

plot_heatmap(
    adult_duration_period.table,
    x="period",
    y="duration",
    title=f"Adult survival: duration-period slice at magnitude = {fixed_magnitude}",
)

Integrated vulnerability by demographic target

Finally, we compute integrated vulnerability for several demographic targets over the same perturbation grid.

full_grid = PerturbationGrid(
    magnitudes=np.linspace(0, 0.8, 17),
    durations=range(0, 7),
    periods=range(1, 13),
)

target_summary = []

for target in ["adult_survival", "juvenile_survival", "fecundity", "all"]:
    result = run_grid(
        model,
        target=target,
        grid=full_grid,
        t_max=60,
        recovery_steps=12,
    )
    target_summary.append({"target": target, "vulnerability": result.vulnerability})

target_summary = pd.DataFrame(target_summary)
target_summary

fig, ax = plt.subplots(figsize=(7, 4))
ax.bar(target_summary["target"], target_summary["vulnerability"])
ax.set_xlabel("Perturbed demographic target")
ax.set_ylabel("Integrated vulnerability Φ (%)")
ax.set_title("Integrated vulnerability across perturbation targets")
ax.tick_params(axis="x", rotation=30)
fig.tight_layout()
plt.show()

Notes

These examples are deliberately small enough to run quickly. For empirical analyses, the same functions can be used with larger projection matrices and denser perturbation grids. The computational cost scales with the number of simulated perturbation regimes and the projection horizon.