un_insure.md: harden fsolve calibration against numpy array→scalar strictness (defensive)

[mmcky]

Summary

This is a defensive / robustness suggestion, not a live bug. The calibration cell in lectures/un_insure.md passes a bare Python scalar as the initial guess to scipy.optimize.fsolve, and its residual function returns array-valued expressions when fsolve hands it a 1-element array. That pattern executes fine on the current scientific stack, but it briefly broke notebook execution downstream under a dependency-version drift, so it is worth hardening.

Where

In lectures/un_insure.md, the autarky calibration uses (around lines 613, 632 and 635):

# inside r_error(self, r)
Vu_star = sp.optimize.fsolve(Vu_error_Λ, 15000, args=(r,))[0]
...
r_calibrated = sp.optimize.brentq(r_error_Λ, 1e-10, 1 - 1e-10)
Vu_aut = sp.optimize.fsolve(Vu_error_Λ, 15000, args=(r_calibrated,))[0]

fsolve internally does x0 = asarray(x0).flatten(), so the iterate Vu arrives at Vu_error as a shape-(1,) array, and the residual is computed as a 1-element array.

Background

The QuantEcon build container (quantecon-build) had its scientific stack unpinned, so it drifted ahead of the anaconda=2025.12 baseline this repo pins. During that drift (≈ March 2026) un_insure.md failed with TypeError: only 0-dimensional arrays can be converted to Python scalars in the brentq → fsolve path, on repos building against the lean container. That is tracked in QuantEcon/actions#28 and the container has now been pinned back to the 2025.12 baseline in QuantEcon/actions#84, which resolves the CI symptom.

The underlying trigger class is real: numpy 2.4.0 promoted converting an ndim > 0 array to a Python scalar (via float() / int()) from a deprecation to a hard TypeError.

Why it is not a live bug

I re-ran the actual calibration cell across four numpy/scipy combinations and it passes on all of them:

numpy	scipy	result
2.3.5	1.16.3 (anaconda 2025.12)	✅
2.4.0	1.15.3	✅
2.4.0	1.16.3	✅
2.4.6	1.17.1 (anaconda 2026.06)	✅

So the specific breaking combination has since been patched upstream and is no longer reproducible. The code is simply fragile to this class of change.

Suggested hardening

Make the calibration robust to the iterate arriving as a 1-element array and ensure scalars are returned, e.g. normalise the input in Vu_error:

def Vu_error(self, Vu, r):
    β, Ve = self.β, self.Ve
    Vu = np.asarray(Vu).item()          # treat the fsolve iterate as a scalar
    a = max(0.0, invp_prime(1 / (β * (Ve - Vu)), r))
    return u(self, 0) - a + β * (p(a, r) * Ve + (1 - p(a, r)) * Vu) - Vu

(Equivalently, pass x0 as a 1-element array and extract results with float(result[0]) / .item().) I verified this variant produces identical results (Vu_aut ≈ 16758.7, p ≈ 0.1) on all four stacks above.

Priority

Low — defensive only. Good first issue. No action is required for current builds; this just prevents a recurrence if the dependency stack moves again.

[mmcky]

Validated — and it turns out this is already fixed by #327, so I'm closing.

I reproduced the failure, traced the exact mechanism, then found the hardening had already landed.

Root cause (numpy ≥ 2.4): numpy 2.4.0 turned converting an ndim > 0 array to a Python scalar (via float()/int()) from a deprecation into a hard TypeError: only 0-dimensional arrays can be converted to Python scalars. Before #327, r_error returned its residual without indexing the fsolve result, so it handed brentq a shape-(1,) array — and brentq calls float() on the residual, which fails on numpy ≥ 2.4.

Minimal reproduction (faithful to the autarky calibration cell):

import numpy as np
import scipy.optimize

class params_instance:
    def __init__(self, r, β=0.999, σ=0.5, w=100):
        self.β, self.σ, self.w, self.r = β, σ, w, r
        self.Ve = (w**(1 - σ) / (1 - σ)) / (1 - β)

def p(a, r):          return 1 - np.exp(-r * a)
def invp_prime(x, r): return -np.log(x / r) / r
def u(self, c):       return c**(1 - self.σ) / (1 - self.σ)

def Vu_error(self, Vu, r):
    a = max(0, invp_prime(1 / (self.β * (self.Ve - Vu)), r))
    return u(self, 0) - a + self.β * (p(a, r) * self.Ve + (1 - p(a, r)) * Vu) - Vu

params = params_instance(r=1e-2)
Vu_error_Λ = lambda Vu, r: Vu_error(params, Vu, r)

def r_error(r):                                                  # pre-#327 form
    Vu_star = scipy.optimize.fsolve(Vu_error_Λ, 15000, args=(r,))   # <-- no [0]; returns a (1,) array
    a_star  = invp_prime(1 / (params.β * (params.Ve - Vu_star)), r)
    return p(a_star, r) - 0.1

scipy.optimize.brentq(r_error, 1e-10, 1 - 1e-10)

Results (the [0] is the fix from #327):

r_error	residual returned	brentq, numpy 2.4.6 / scipy 1.17.1
fsolve(...) (pre-#327)	ndarray, shape (1,)	❌ TypeError: only 0-dimensional arrays can be converted to Python scalars
fsolve(...)[0] (current main)	float64, shape ()	✅ r = 0.000343…

Already resolved by #327 (merged 2026-05-28): indexing fsolve(...)[0] makes r_error return a 0-dimensional scalar, and float() of a 0-d value succeeds on every numpy version (numpy 2.4 only rejects ndim > 0 arrays). I confirmed current main passes on both numpy 2.3.5 / scipy 1.16.3 (Anaconda 2025.12) and numpy 2.4.6 / scipy 1.17.1 (Anaconda 2026.06), so the fix is version-independent and safe for a future Anaconda 2026.06 upgrade.

Closing as resolved. For context, the container-side change that surfaced this downstream is tracked in QuantEcon/actions#28 (lean build container's scientific stack had drifted to numpy 2.4; pinned back to the 2025.12 baseline in QuantEcon/actions#84).