Tutorial: Supply Chain Reverse Stress Testing via Differentiable ADTs¶

Supply chain risk and financial contagion look like completely different problems. But they share a single algebraic spine — and once you see it, you understand both.

The duality:

	Financial contagion	Supply chain
Logic	OR — any counterparty fails	AND — every tier must deliver
Math	Sum (accumulation)	Min (bottleneck)
Type theory	Sum type (Union)	Product type (Record)
PCL combinator	`parallel(flow_A, flow_B)`	`fold(β→∞, [tier_1, ..., tier_n])`
Failure mode	Contagion	Starvation

This is the Curry-Howard correspondence applied to economic risk: the same algebraic structure that connects OR to Sum types in type theory connects financial contagion to additive aggregation in risk modelling.

What this tutorial demonstrates:

Network topology: the copper supply chain as a Pacioli manifold — 8 nodes, 4 tiers, Laplacian spectrum
Forward simulation: how a mine disruption propagates to the OEM via the AND layer
Shock scenarios: CHL_Mine at 50%, 20%, 0% — and the dual-mine shock
Reverse Stress Testing (RST): JAX gradient descent finds the minimum inventory buffer to survive a given shock
Criticality: the gradient $\partial\text{output}/\partial\text{capacity}_i$ identifies the binding constraint — the 'straw that breaks the camel's back'
PCL connection: every supply chain is a PCL computation tree; fold is AND, parallel is OR

References: Smith et al. (2025) arXiv:2511.07289; Buckley (2026) PCL doi:10.5281/zenodo.20262070

In [ ]:

Copied!

# Uncomment to install in Colab
# !pip install econiac jax[cpu] matplotlib
# Uncomment to install in Colab
# !pip install econiac jax[cpu] matplotlib

In [ ]:

Copied!





import numpy as np
import matplotlib.pyplot as plt
import matplotlib.colors as mcolors

from econiac.economics.supply_chain import (
    COPPER_CHAIN,
    COPPER_NODES,
    simulate_chain,
    reverse_stress_test,
    laplacian_spectrum,
    apply_shock,
    to_pcl_description,
)

plt.rcParams.update({'figure.dpi': 120, 'font.size': 11})
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.colors as mcolors

from econiac.economics.supply_chain import (
    COPPER_CHAIN,
    COPPER_NODES,
    simulate_chain,
    reverse_stress_test,
    laplacian_spectrum,
    apply_shock,
    to_pcl_description,
)

plt.rcParams.update({'figure.dpi': 120, 'font.size': 11})

1. Network topology: the copper supply chain¶

The copper wire supply chain (after Smith et al. 2025) has four tiers:

Tier 0 (demand):    US_OEM
Tier 1 (refining):  US_Refinery  CAN_Refinery  CHL_Refinery
Tier 2 (smelting):  CAN_Smelter  CHL_Smelter
Tier 3 (mining):    CHL_Mine     DRC_Mine

The OEM (Tier 0) aggregates supply from three refineries via market weights — this is the OR/Sum layer (financial-risk style: more sources = more resilient).

The refineries and smelters depend on upstream tiers via Leontief constraints — this is the AND/Min layer (supply-chain style: you need every required input).

The hybrid structure — OR at the top, AND at the bottom — is typical of real supply chains.

In [ ]:

Copied!

print(to_pcl_description(COPPER_CHAIN))
print(to_pcl_description(COPPER_CHAIN))

Laplacian spectrum: where is the bottleneck?¶

The Laplacian spectrum of the supply chain graph gives the same topological information as ModelLG.laplacian_spectrum() does for the financial circuit:

λ₂ (Fiedler value) = algebraic connectivity — how quickly a shock propagates
Fiedler vector = the cut that divides the graph at its weakest link

Nodes on the negative side of the Fiedler vector are 'upstream' of the most vulnerable cut. A disruption there propagates farthest.

In [ ]:

Copied!





spec = laplacian_spectrum(COPPER_CHAIN)

fig, axes = plt.subplots(1, 2, figsize=(12, 4))

# Eigenvalue spectrum
ax = axes[0]
eigvals = spec['eigenvalues']
ax.stem(range(len(eigvals)), eigvals, basefmt='k-')
ax.axhline(spec['fiedler_value'], color='crimson', linestyle='--', label=f"λ₂ = {spec['fiedler_value']:.3f}")
ax.set_xlabel('Index')
ax.set_ylabel('Eigenvalue')
ax.set_title('Laplacian spectrum')
ax.legend()

# Fiedler vector
ax = axes[1]
fv = spec['fiedler_vector']
order = np.argsort(fv)
colors = ['#d62728' if v < 0 else '#1f77b4' for v in fv[order]]
ax.barh([COPPER_NODES[i] for i in order], fv[order], color=colors)
ax.axvline(0, color='black', linewidth=0.8)
ax.set_xlabel('Fiedler component')
ax.set_title('Fiedler vector\n(red = upstream of weakest cut)')

plt.tight_layout()
plt.show()

print(f"Connected components: {spec['n_components']}")
print(f"Fiedler value λ₂ = {spec['fiedler_value']:.4f}  (algebraic connectivity)")
spec = laplacian_spectrum(COPPER_CHAIN)

fig, axes = plt.subplots(1, 2, figsize=(12, 4))

# Eigenvalue spectrum
ax = axes[0]
eigvals = spec['eigenvalues']
ax.stem(range(len(eigvals)), eigvals, basefmt='k-')
ax.axhline(spec['fiedler_value'], color='crimson', linestyle='--', label=f"λ₂ = {spec['fiedler_value']:.3f}")
ax.set_xlabel('Index')
ax.set_ylabel('Eigenvalue')
ax.set_title('Laplacian spectrum')
ax.legend()

# Fiedler vector
ax = axes[1]
fv = spec['fiedler_vector']
order = np.argsort(fv)
colors = ['#d62728' if v < 0 else '#1f77b4' for v in fv[order]]
ax.barh([COPPER_NODES[i] for i in order], fv[order], color=colors)
ax.axvline(0, color='black', linewidth=0.8)
ax.set_xlabel('Fiedler component')
ax.set_title('Fiedler vector\n(red = upstream of weakest cut)')

plt.tight_layout()
plt.show()

print(f"Connected components: {spec['n_components']}")
print(f"Fiedler value λ₂ = {spec['fiedler_value']:.4f}  (algebraic connectivity)")

2. Forward simulation: capacity propagation¶

The simulation unrolls n_steps of the supply chain dynamics:

capacity[t+1] = AND_step(capacity[t] + buffers, BOM)   # Leontief constraint
capacity[t+1] = OR_step(capacity[t], market_weights)    # Market aggregation

Starting from full capacity, the SoftMin (temperature=0.1) introduces a small approximation gap — the network settles to a steady state just below 1.0 at nodes that aggregate from multiple suppliers (the CHL_Smelter requires two mines; SoftMin of two 1.0s at T=0.1 is $1 - \ln 2 / 10 \approx 0.931$).

In [ ]:

Copied!





baseline = np.ones(len(COPPER_NODES))
result_base = simulate_chain(COPPER_CHAIN, baseline)

fig, ax = plt.subplots(figsize=(10, 4))
history = result_base['capacity_history']

colors = plt.cm.tab10(np.linspace(0, 1, len(COPPER_NODES)))
for i, name in enumerate(COPPER_NODES):
    ax.plot(history[:, i], label=name, color=colors[i])

ax.set_xlabel('Simulation step')
ax.set_ylabel('Capacity')
ax.set_title('Baseline capacity propagation (full capacity, no shock)')
ax.set_ylim(0, 1.05)
ax.legend(bbox_to_anchor=(1.01, 1), loc='upper left', fontsize=9)
plt.tight_layout()
plt.show()

print("Final capacities:")
for name, cap in zip(COPPER_NODES, result_base['final_capacity']):
    bar = '█' * int(cap * 30)
    print(f"  {name:22s}  {cap:.3f}  {bar}")
baseline = np.ones(len(COPPER_NODES))
result_base = simulate_chain(COPPER_CHAIN, baseline)

fig, ax = plt.subplots(figsize=(10, 4))
history = result_base['capacity_history']

colors = plt.cm.tab10(np.linspace(0, 1, len(COPPER_NODES)))
for i, name in enumerate(COPPER_NODES):
    ax.plot(history[:, i], label=name, color=colors[i])

ax.set_xlabel('Simulation step')
ax.set_ylabel('Capacity')
ax.set_title('Baseline capacity propagation (full capacity, no shock)')
ax.set_ylim(0, 1.05)
ax.legend(bbox_to_anchor=(1.01, 1), loc='upper left', fontsize=9)
plt.tight_layout()
plt.show()

print("Final capacities:")
for name, cap in zip(COPPER_NODES, result_base['final_capacity']):
    bar = '█' * int(cap * 30)
    print(f"  {name:22s}  {cap:.3f}  {bar}")

3. Shock scenarios: disruption at CHL_Mine¶

Chile dominates global copper ore production. What happens when CHL_Mine is disrupted — whether by a strike, a tailing collapse, or geopolitical action?

The AND layer propagates the shortfall upward:

CHL_Mine ↓  →  CHL_Smelter ↓  →  CHL_Refinery ↓  →  US_OEM ↓
             (also)  CAN_Smelter ↓  →  US_Refinery, CAN_Refinery ↓

The OEM has partial resilience via the market-weights OR layer — it can substitute toward US_Refinery and CAN_Refinery. But if both CHL_Mine and DRC_Mine are disrupted, the CAN_Smelter (which sources from CHL_Mine) also fails, collapsing all three refineries.

In [ ]:

Copied!





shock_levels = [1.0, 0.75, 0.5, 0.25, 0.0]
oem_outputs = []
node_finals = []

for sl in shock_levels:
    shocked = apply_shock(COPPER_CHAIN, shocked_nodes=['CHL_Mine'], shock_fraction=sl)
    result = simulate_chain(COPPER_CHAIN, shocked)
    oem_outputs.append(result['final_capacity'][0])
    node_finals.append(result['final_capacity'])

fig, axes = plt.subplots(1, 2, figsize=(13, 5))

# OEM output vs shock level
ax = axes[0]
ax.plot([1 - sl for sl in shock_levels], oem_outputs, 'o-', color='crimson', linewidth=2)
ax.axhline(0.8, color='gray', linestyle='--', label='80% threshold')
ax.set_xlabel('CHL_Mine disruption (fraction lost)')
ax.set_ylabel('US_OEM final capacity')
ax.set_title('OEM output vs. CHL_Mine shock\n(AND layer propagation)')
ax.legend()
ax.set_xlim(0, 1)
ax.set_ylim(0, 1.05)

# Heatmap: all nodes vs shock levels
ax = axes[1]
matrix = np.array(node_finals)
im = ax.imshow(matrix.T, aspect='auto', vmin=0, vmax=1, cmap='RdYlGn')
ax.set_xticks(range(len(shock_levels)))
ax.set_xticklabels([f'{int((1-sl)*100)}%' for sl in shock_levels])
ax.set_xlabel('CHL_Mine disruption (%)')
ax.set_yticks(range(len(COPPER_NODES)))
ax.set_yticklabels(COPPER_NODES, fontsize=9)
ax.set_title('Capacity heatmap by node and shock')
plt.colorbar(im, ax=ax, label='Capacity')

plt.tight_layout()
plt.show()
shock_levels = [1.0, 0.75, 0.5, 0.25, 0.0]
oem_outputs = []
node_finals = []

for sl in shock_levels:
    shocked = apply_shock(COPPER_CHAIN, shocked_nodes=['CHL_Mine'], shock_fraction=sl)
    result = simulate_chain(COPPER_CHAIN, shocked)
    oem_outputs.append(result['final_capacity'][0])
    node_finals.append(result['final_capacity'])

fig, axes = plt.subplots(1, 2, figsize=(13, 5))

# OEM output vs shock level
ax = axes[0]
ax.plot([1 - sl for sl in shock_levels], oem_outputs, 'o-', color='crimson', linewidth=2)
ax.axhline(0.8, color='gray', linestyle='--', label='80% threshold')
ax.set_xlabel('CHL_Mine disruption (fraction lost)')
ax.set_ylabel('US_OEM final capacity')
ax.set_title('OEM output vs. CHL_Mine shock\n(AND layer propagation)')
ax.legend()
ax.set_xlim(0, 1)
ax.set_ylim(0, 1.05)

# Heatmap: all nodes vs shock levels
ax = axes[1]
matrix = np.array(node_finals)
im = ax.imshow(matrix.T, aspect='auto', vmin=0, vmax=1, cmap='RdYlGn')
ax.set_xticks(range(len(shock_levels)))
ax.set_xticklabels([f'{int((1-sl)*100)}%' for sl in shock_levels])
ax.set_xlabel('CHL_Mine disruption (%)')
ax.set_yticks(range(len(COPPER_NODES)))
ax.set_yticklabels(COPPER_NODES, fontsize=9)
ax.set_title('Capacity heatmap by node and shock')
plt.colorbar(im, ax=ax, label='Capacity')

plt.tight_layout()
plt.show()

In [ ]:

Copied!





# Dual-mine shock: CHL_Mine + DRC_Mine both disrupted
print("Dual mining shock scenarios (CHL_Mine = DRC_Mine = shock_fraction):")
print(f"  {'Shock':10s}  {'OEM output':>11}  {'vs. single':>12}")
print(f"  {'-'*36}")

for sl in [0.5, 0.2, 0.0]:
    single = apply_shock(COPPER_CHAIN, shocked_nodes=['CHL_Mine'], shock_fraction=sl)
    dual = apply_shock(COPPER_CHAIN, shocked_nodes=['CHL_Mine', 'DRC_Mine'], shock_fraction=sl)
    r_single = simulate_chain(COPPER_CHAIN, single)['final_capacity'][0]
    r_dual = simulate_chain(COPPER_CHAIN, dual)['final_capacity'][0]
    print(f"  {int((1-sl)*100):3d}% loss  {r_single:10.3f}  {r_dual:10.3f}  (dual {r_dual-r_single:+.3f})")

print()
print("Insight: the DRC_Mine matters to CHL_Smelter (AND gate),")
print("so dual disruption is super-additive — worse than the sum of parts.")
# Dual-mine shock: CHL_Mine + DRC_Mine both disrupted
print("Dual mining shock scenarios (CHL_Mine = DRC_Mine = shock_fraction):")
print(f"  {'Shock':10s}  {'OEM output':>11}  {'vs. single':>12}")
print(f"  {'-'*36}")

for sl in [0.5, 0.2, 0.0]:
    single = apply_shock(COPPER_CHAIN, shocked_nodes=['CHL_Mine'], shock_fraction=sl)
    dual = apply_shock(COPPER_CHAIN, shocked_nodes=['CHL_Mine', 'DRC_Mine'], shock_fraction=sl)
    r_single = simulate_chain(COPPER_CHAIN, single)['final_capacity'][0]
    r_dual = simulate_chain(COPPER_CHAIN, dual)['final_capacity'][0]
    print(f"  {int((1-sl)*100):3d}% loss  {r_single:10.3f}  {r_dual:10.3f}  (dual {r_dual-r_single:+.3f})")

print()
print("Insight: the DRC_Mine matters to CHL_Smelter (AND gate),")
print("so dual disruption is super-additive — worse than the sum of parts.")

4. Reverse Stress Testing (RST)¶

The RST question: given a shock, what is the minimum inventory buffer at each node that ensures the OEM achieves a required output level?

This is the 'straw that breaks the camel's back' problem in reverse: instead of asking 'what shock breaks the chain?', we ask 'what buffer saves it?'

The solution is a differentiable optimisation:

$$\text{buffers}^* = \arg\min_{b \geq 0} \; \lambda \sum_i b_i + 10 \cdot \text{ReLU}(q^* - \hat{y}(b))$$

where $q^*$ is the required output and $\hat{y}(b)$ is the simulated OEM capacity with buffers $b$. JAX grad() propagates through the entire SoftMin chain, identifying which nodes are binding constraints (high gradient = high criticality).

In [ ]:

Copied!





shock_50 = apply_shock(COPPER_CHAIN, shocked_nodes=['CHL_Mine'], shock_fraction=0.5)

rst_results = {}
for required in [0.70, 0.85, 0.95]:
    rst_results[required] = reverse_stress_test(
        COPPER_CHAIN,
        shock_capacity=shock_50,
        required_output=required,
        output_node=0,
        budget_weight=0.1,
        n_epochs=400,
    )

fig, axes = plt.subplots(1, 3, figsize=(14, 5))

for ax, (req, rst) in zip(axes, rst_results.items()):
    # Loss convergence
    ax.plot(rst['loss_history'], color='steelblue', linewidth=1.5)
    title = f'Required ≥ {req:.0%}\n'
    title += '✓ survived' if rst['survived'] else '✗ not achieved'
    title += f'  (output={rst["final_output"]:.3f})'
    ax.set_title(title, fontsize=10)
    ax.set_xlabel('Epoch')
    ax.set_ylabel('RST loss')
    ax.set_yscale('log')

plt.suptitle('RST convergence: CHL_Mine shocked to 50%', y=1.02)
plt.tight_layout()
plt.show()
shock_50 = apply_shock(COPPER_CHAIN, shocked_nodes=['CHL_Mine'], shock_fraction=0.5)

rst_results = {}
for required in [0.70, 0.85, 0.95]:
    rst_results[required] = reverse_stress_test(
        COPPER_CHAIN,
        shock_capacity=shock_50,
        required_output=required,
        output_node=0,
        budget_weight=0.1,
        n_epochs=400,
    )

fig, axes = plt.subplots(1, 3, figsize=(14, 5))

for ax, (req, rst) in zip(axes, rst_results.items()):
    # Loss convergence
    ax.plot(rst['loss_history'], color='steelblue', linewidth=1.5)
    title = f'Required ≥ {req:.0%}\n'
    title += '✓ survived' if rst['survived'] else '✗ not achieved'
    title += f'  (output={rst["final_output"]:.3f})'
    ax.set_title(title, fontsize=10)
    ax.set_xlabel('Epoch')
    ax.set_ylabel('RST loss')
    ax.set_yscale('log')

plt.suptitle('RST convergence: CHL_Mine shocked to 50%', y=1.02)
plt.tight_layout()
plt.show()

In [ ]:

Copied!





# Buffer allocation and criticality for 85% target
rst = rst_results[0.85]

fig, axes = plt.subplots(1, 2, figsize=(12, 5))

# Buffer allocation
ax = axes[0]
order = np.argsort(-rst['buffers'])
colors = ['#2ca02c' if b > 0.001 else '#aec7e8' for b in rst['buffers'][order]]
ax.barh([COPPER_NODES[i] for i in order], rst['buffers'][order], color=colors)
ax.set_xlabel('Buffer size')
ax.set_title(f'Optimal buffer allocation\n(total cost = {rst["buffers"].sum():.3f})')

# Criticality
ax = axes[1]
crit_order = np.argsort(-rst['criticality'])
colors = ['#d62728' if c > 0.05 else '#aec7e8' for c in rst['criticality'][crit_order]]
ax.barh([COPPER_NODES[i] for i in crit_order], rst['criticality'][crit_order], color=colors)
ax.set_xlabel('|∂output / ∂capacity|')
ax.set_title('Node criticality\n(gradient of OEM output w.r.t. node capacity)')

plt.suptitle('RST result: 85% OEM output target, CHL_Mine at 50%', y=1.02)
plt.tight_layout()
plt.show()

print("Interpretation:")
top_crit = COPPER_NODES[np.argmax(rst['criticality'])]
top_buf = COPPER_NODES[np.argmax(rst['buffers'])]
print(f"  Most critical node (by gradient): {top_crit}")
print(f"  Highest buffered node:            {top_buf}")
print(f"  These should agree — the gradient guides the optimizer to the right node.")
# Buffer allocation and criticality for 85% target
rst = rst_results[0.85]

fig, axes = plt.subplots(1, 2, figsize=(12, 5))

# Buffer allocation
ax = axes[0]
order = np.argsort(-rst['buffers'])
colors = ['#2ca02c' if b > 0.001 else '#aec7e8' for b in rst['buffers'][order]]
ax.barh([COPPER_NODES[i] for i in order], rst['buffers'][order], color=colors)
ax.set_xlabel('Buffer size')
ax.set_title(f'Optimal buffer allocation\n(total cost = {rst["buffers"].sum():.3f})')

# Criticality
ax = axes[1]
crit_order = np.argsort(-rst['criticality'])
colors = ['#d62728' if c > 0.05 else '#aec7e8' for c in rst['criticality'][crit_order]]
ax.barh([COPPER_NODES[i] for i in crit_order], rst['criticality'][crit_order], color=colors)
ax.set_xlabel('|∂output / ∂capacity|')
ax.set_title('Node criticality\n(gradient of OEM output w.r.t. node capacity)')

plt.suptitle('RST result: 85% OEM output target, CHL_Mine at 50%', y=1.02)
plt.tight_layout()
plt.show()

print("Interpretation:")
top_crit = COPPER_NODES[np.argmax(rst['criticality'])]
top_buf = COPPER_NODES[np.argmax(rst['buffers'])]
print(f"  Most critical node (by gradient): {top_crit}")
print(f"  Highest buffered node:            {top_buf}")
print(f"  These should agree — the gradient guides the optimizer to the right node.")

5. Dual-mine RST: when the Fiedler cut is breached¶

The Laplacian analysis in Section 1 identified DRC_Mine and CHL_Smelter on the negative (upstream) side of the Fiedler cut. When both mines are disrupted, the cut is breached — CHL_Smelter loses both inputs simultaneously.

How much buffering does it take to restore 80% OEM output after both mines fail?

In [ ]:

Copied!





shock_dual = apply_shock(
    COPPER_CHAIN,
    shocked_nodes=['CHL_Mine', 'DRC_Mine'],
    shock_fraction=0.2,
)

rst_dual = reverse_stress_test(
    COPPER_CHAIN,
    shock_capacity=shock_dual,
    required_output=0.80,
    output_node=0,
    budget_weight=0.05,
    n_epochs=500,
)

fig, axes = plt.subplots(1, 2, figsize=(12, 5))

ax = axes[0]
ax.plot(rst_dual['loss_history'], color='crimson')
ax.set_xlabel('Epoch')
ax.set_ylabel('RST loss')
ax.set_yscale('log')
result_str = '✓' if rst_dual['survived'] else '✗'
ax.set_title(f'Convergence: dual mine shock\n{result_str} final output = {rst_dual["final_output"]:.3f}')

ax = axes[1]
order = np.argsort(-rst_dual['criticality'])
ax.barh(
    [COPPER_NODES[i] for i in order],
    rst_dual['criticality'][order],
    color='#d62728'
)
ax.set_xlabel('Criticality')
ax.set_title(f'Node criticality: dual mine shock\n(total buffer cost = {rst_dual["buffers"].sum():.3f})')

plt.tight_layout()
plt.show()

print(f"Dual-mine RST summary:")
print(f"  Both mines shocked to 20% capacity")
print(f"  Required OEM output: 80%")
print(f"  Achieved: {rst_dual['final_output']:.3f}  ({'survived' if rst_dual['survived'] else 'failed'})")
print(f"  Total buffer cost: {rst_dual['buffers'].sum():.4f}")
print()
print("  Buffer allocation:")
for i in np.argsort(-rst_dual['buffers'])[:4]:
    if rst_dual['buffers'][i] > 0.001:
        print(f"    {COPPER_NODES[i]:22s}  {rst_dual['buffers'][i]:.4f}")
shock_dual = apply_shock(
    COPPER_CHAIN,
    shocked_nodes=['CHL_Mine', 'DRC_Mine'],
    shock_fraction=0.2,
)

rst_dual = reverse_stress_test(
    COPPER_CHAIN,
    shock_capacity=shock_dual,
    required_output=0.80,
    output_node=0,
    budget_weight=0.05,
    n_epochs=500,
)

fig, axes = plt.subplots(1, 2, figsize=(12, 5))

ax = axes[0]
ax.plot(rst_dual['loss_history'], color='crimson')
ax.set_xlabel('Epoch')
ax.set_ylabel('RST loss')
ax.set_yscale('log')
result_str = '✓' if rst_dual['survived'] else '✗'
ax.set_title(f'Convergence: dual mine shock\n{result_str} final output = {rst_dual["final_output"]:.3f}')

ax = axes[1]
order = np.argsort(-rst_dual['criticality'])
ax.barh(
    [COPPER_NODES[i] for i in order],
    rst_dual['criticality'][order],
    color='#d62728'
)
ax.set_xlabel('Criticality')
ax.set_title(f'Node criticality: dual mine shock\n(total buffer cost = {rst_dual["buffers"].sum():.3f})')

plt.tight_layout()
plt.show()

print(f"Dual-mine RST summary:")
print(f"  Both mines shocked to 20% capacity")
print(f"  Required OEM output: 80%")
print(f"  Achieved: {rst_dual['final_output']:.3f}  ({'survived' if rst_dual['survived'] else 'failed'})")
print(f"  Total buffer cost: {rst_dual['buffers'].sum():.4f}")
print()
print("  Buffer allocation:")
for i in np.argsort(-rst_dual['buffers'])[:4]:
    if rst_dual['buffers'][i] > 0.001:
        print(f"    {COPPER_NODES[i]:22s}  {rst_dual['buffers'][i]:.4f}")

6. The PCL connection: every supply chain is a computation tree¶

The deepest insight: a supply chain network is a PCL computation tree.

Each tier maps to a PCL combinator:

AND tier (Leontief BOM):   fold(β→∞, [supplier_1, supplier_2, ...])
OR tier (market weights):  parallel(flow_supplier_A, flow_supplier_B, ...)
Sequential tiers:          sequence(tier_3, tier_2, tier_1, tier_0)
Buffer injection:          flow(mine, smelter, "inventory", qty)

The RST problem is exactly conservation_loss calibration on the Pacioli manifold: find the minimum flow(mine, smelter, ...) transfers that maintain output above threshold.

The Curry-Howard table¶

Economic concept	Type theory	Logic	Math	PCL
Supply chain dependency	Product type (Record)	AND (∧)	min	`fold(β→∞, ...)`
Market sourcing	Sum type (Union)	OR (∨)	sum	`parallel(...)`
Inventory buffer	Term (value)	Proof	constant	`flow(...)`
Supply chain network	Function type	Implication	composition	`sequence(...)`
RST optimisation	Type inhabitation	Proof search	gradient descent	`conservation_loss`

The rightmost column shows that RST is proof search in PCL: finding the minimum set of flow terms that make the supply chain computation well-typed (output ≥ threshold).

In [ ]:

Copied!





# Demonstrate that temperature controls the AND/OR sharpness
# Lower T = sharper min (true Leontief); higher T = smoother (average)
import jax.numpy as jnp
from econiac.economics.supply_chain import soft_min_axis

x = jnp.array([0.5, 0.9, 1.0])  # Three suppliers, one at 50%
temperatures = [0.01, 0.05, 0.1, 0.5, 1.0, 5.0]

print("SoftMin([0.5, 0.9, 1.0]) at different temperatures:")
print(f"  True min = {float(x.min()):.3f}   Mean = {float(x.mean()):.3f}")
print()
print(f"  {'Temperature':>12}  {'SoftMin':>8}  {'Interpretation'}")
print(f"  {'-'*50}")
for T in temperatures:
    val = float(soft_min_axis(x, temperature=T))
    if T <= 0.05:
        interp = "≈ true min (strict Leontief)"
    elif T >= 2.0:
        interp = "≈ mean (full substitutability)"
    else:
        interp = "smooth interpolation"
    print(f"  {T:>12.2f}  {val:>8.4f}  {interp}")

print()
print("Temperature is the interpolation parameter between:")
print("  T→0: pure AND (Leontief, no substitution possible)")
print("  T→∞: pure average (perfect substitutability)")
print()
print("Real supply chains are in between: partial substitution exists,")
print("but switching suppliers has cost and lead time (finite T).")
# Demonstrate that temperature controls the AND/OR sharpness
# Lower T = sharper min (true Leontief); higher T = smoother (average)
import jax.numpy as jnp
from econiac.economics.supply_chain import soft_min_axis

x = jnp.array([0.5, 0.9, 1.0])  # Three suppliers, one at 50%
temperatures = [0.01, 0.05, 0.1, 0.5, 1.0, 5.0]

print("SoftMin([0.5, 0.9, 1.0]) at different temperatures:")
print(f"  True min = {float(x.min()):.3f}   Mean = {float(x.mean()):.3f}")
print()
print(f"  {'Temperature':>12}  {'SoftMin':>8}  {'Interpretation'}")
print(f"  {'-'*50}")
for T in temperatures:
    val = float(soft_min_axis(x, temperature=T))
    if T <= 0.05:
        interp = "≈ true min (strict Leontief)"
    elif T >= 2.0:
        interp = "≈ mean (full substitutability)"
    else:
        interp = "smooth interpolation"
    print(f"  {T:>12.2f}  {val:>8.4f}  {interp}")

print()
print("Temperature is the interpolation parameter between:")
print("  T→0: pure AND (Leontief, no substitution possible)")
print("  T→∞: pure average (perfect substitutability)")
print()
print("Real supply chains are in between: partial substitution exists,")
print("but switching suppliers has cost and lead time (finite T).")

Summary¶

What we built:

Hybrid AND/OR supply chain model — Leontief constraints at the physical tier, market-weighted aggregation at the demand tier
Laplacian spectrum — identifies the Fiedler cut (weakest link) before any shock
Shock propagation — CHL_Mine disruption propagates via the AND layer, partially mitigated by the OR layer at the OEM
RST via JAX gradient descent — finds minimum inventory buffers to survive a shock
Criticality — the gradient $\partial\text{output}/\partial\text{capacity}_i$ is the supply chain analogue of DebtRank in financial contagion

The algebraic spine:

Supply chain  =  fold(β→∞, tiers)              # AND / Product type
Contagion     =  parallel(flows)               # OR  / Sum type
RST           =  argmin conservation_loss      # Proof search / type inhabitation

Connections to other econiac tutorials:

LowGrow SSE: laplacian_spectrum() finds the bond-spiral cycle; the same Fiedler analysis identifies fragile financial circuits
GL Model PC: gibbs_weights is the continuous limit of the AND/OR duality — at β→∞ it becomes the argmax/argmin; at β=0 it is the mean
GEMMES: CurvedBalanceSheet curvature from climate damage is physically analogous to supply chain capacity loss — both are field-strength terms on the Pacioli manifold that no sector can claim as an asset

References: Smith et al. (2025) arXiv:2511.07289 — Buckley (2026) PCL doi:10.5281/zenodo.20262070