mpfit_nan_guard: sanitize non-finite numerical derivatives in mp_fdjac2

mp_fdjac2's finite-difference Jacobian writes were only protected by
assert(isfinite(fjac[ij])), which -DNDEBUG (the Pi's Release build) strips
entirely. A degenerate residual (e.g. right after a Reinit-triggered scene
reset) can make wa[i]/fvec[i] non-finite, writing NaN/Inf straight into
fjac. That poisons fnorm/ratio in mpfit()'s outer Levenberg-Marquardt loop;
since every comparison against a NaN ratio is false, the iteration counter
(gated on the ratio>=p0001 success branch) never advances and the maxiter
check never fires. Confirmed via reproduction: ~465k spin iterations/sec
with zero progress, matching the gdb captures of agent threads stuck at
95-100% CPU for minutes after Reinit.

Sanitize fjac to 0 on a non-finite derivative in both the one-sided and
two-sided (which had no assert at all) paths, so the column is treated the
same way mp_qrfac already handles a zero-norm column.

Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>

Author: rocketmark

redist/mpfit/mpfit.c

@@ -1228,6 +1228,23 @@ static int mp_fdjac2(mp_func funct, int m, int n, int *ifree, int npar, FLT *x,
 					/* Non-debug path for speed */
 					for (i = 0; i < m; i++, ij++) {
 						fjac[ij] = (wa[i] - fvec[i]) / h; /* fjac[i+m*j] */
+						/* Stagehand patch: the assert below is the only protection
+						 * against a non-finite derivative here, and asserts are
+						 * compiled out under -DNDEBUG (the Pi's Release build). A
+						 * degenerate user function (e.g. a residual that goes Inf
+						 * right after a Reinit-triggered scene reset) then writes
+						 * NaN/Inf straight into fjac, which poisons ratio/fnorm in
+						 * mpfit()'s outer loop. Every comparison against a NaN
+						 * ratio is false, so the iteration counter -- which only
+						 * advances inside the ratio>=p0001 success branch -- never
+						 * increments and the maxiter check never fires: an
+						 * unbounded CPU-spin hang instead of a clean failure
+						 * return. Sanitize to zero so the column is treated as
+						 * having no measurable gradient here, same outcome
+						 * mp_qrfac already handles safely via its ajnorm==0
+						 * guard. */
+						if (!isfinite(fjac[ij]))
+							fjac[ij] = 0;
 						assert(isfinite(fjac[ij]));
 					}
 				} else {
@@ -1264,6 +1281,11 @@ static int mp_fdjac2(mp_func funct, int m, int n, int *ifree, int npar, FLT *x,
 					/* Non-debug path for speed */
 					for (i = 0; i < m; i++, ij++) {
 						fjac[ij] = (fjac[ij] - wa[i]) / (2 * h); /* fjac[i+m*j] */
+						/* Stagehand patch: same non-finite-derivative hang as the
+						 * one-sided path above -- this path had no assert at
+						 * all, debug or release. */
+						if (!isfinite(fjac[ij]))
+							fjac[ij] = 0;
 					}
 				} else {
 					/* Debug path for correctness */

src/test_cases/CMakeLists.txt

@@ -7,7 +7,7 @@ SET(SURVIVE_TESTS
         reproject_residual_props event_queue_props
         residual_cascade_props
         variance_gate_props
-        normal_filter_props)
+        normal_filter_props mpfit_props)
 
 set(barycentric_svd_ADDITIONAL_SRCS ../barycentric_svd/barycentric_svd.c)
 set(reproject_residual_props_ADDITIONAL_SRCS ../barycentric_svd/barycentric_svd.c)

src/test_cases/mpfit_props.c

@@ -0,0 +1,98 @@
+// Regression test for the mpfit numerical-derivative NaN hang (Stagehand patch).
+//
+// redist/mpfit/mpfit.c's mp_fdjac2() computes a finite-difference Jacobian
+// when no analytic derivative is supplied. Its only protection against a
+// degenerate (non-finite) residual poisoning fjac was assert(isfinite(...)),
+// which -DNDEBUG (the Pi's Release build) strips entirely. A NaN/Inf in fjac
+// poisons fnorm/ratio in mpfit()'s outer Levenberg-Marquardt loop; since the
+// iteration counter only advances inside the `ratio >= p0001` success branch,
+// and every comparison against a NaN ratio is false, the loop never
+// increments `iter`, the maxiter check never fires, and mpfit() spins at
+// ~100% CPU forever instead of returning a clean failure status.
+//
+// This test drives mpfit() through the numerical-derivative path with a user
+// function that returns a non-finite residual once perturbed by the
+// finite-difference step, and checks mpfit() returns promptly with a finite
+// result. It also confirms a well-conditioned fit on the same code path still
+// converges normally (no regression on the happy path).
+
+#include "test_case.h"
+#include <math.h>
+#include "mpfit/mpfit.h"
+#include <stdio.h>
+
+// y = p0 * x, with no analytic derivative supplied (dvec ignored) so mpfit
+// falls back to mp_fdjac2's numerical derivative. Once p0 is perturbed past
+// a threshold the residual goes non-finite, simulating a degenerate
+// real-world residual (e.g. right after a Reinit-triggered scene reset).
+struct LinearData {
+	FLT *x;
+	FLT *y;
+	int n;
+};
+
+static int linear_degenerate(int m, int n, FLT *p, FLT *dy, FLT **dvec, void *priv) {
+	struct LinearData *d = (struct LinearData *)priv;
+	for (int i = 0; i < m; i++) {
+		if (p[0] > 1e10) {
+			dy[i] = INFINITY;
+		} else {
+			dy[i] = d->y[i] - p[0] * d->x[i];
+		}
+	}
+	return 0;
+}
+
+static int linear_well_conditioned(int m, int n, FLT *p, FLT *dy, FLT **dvec, void *priv) {
+	struct LinearData *d = (struct LinearData *)priv;
+	for (int i = 0; i < m; i++) {
+		dy[i] = d->y[i] - (p[0] * d->x[i] + p[1]);
+	}
+	return 0;
+}
+
+TEST(MpfitProps, NonFiniteResidualDoesNotHang) {
+	FLT x[5] = {0, 1, 2, 3, 4};
+	FLT y[5] = {0, 2, 4, 6, 8};
+	struct LinearData d = {x, y, 5};
+
+	// Start absurdly large so the finite-difference probe step pushes p[0]
+	// past the 1e10 threshold and the user function returns Inf.
+	FLT p[1] = {1e15};
+	mp_result result = {0};
+	mp_config cfg = {0};
+	cfg.maxiter = 50;
+
+	int status = mpfit(linear_degenerate, 5, 1, p, 0, &cfg, &d, &result);
+
+	// The point of the patch is that this returns at all (no hang) and
+	// doesn't leave NaN/Inf in the solver's outputs.
+	if (!isfinite(result.bestnorm) || !isfinite(p[0])) {
+		fprintf(stderr, "NonFiniteResidualDoesNotHang FAILED: status=%d bestnorm=%f p0=%f\n", status, result.bestnorm,
+				p[0]);
+		return -1;
+	}
+	return 0;
+}
+
+TEST(MpfitProps, WellConditionedFitStillConverges) {
+	// y = 2x + 1, same numerical-derivative code path as the test above.
+	FLT x[6] = {0, 1, 2, 3, 4, 5};
+	FLT y[6] = {1, 3, 5, 7, 9, 11};
+	struct LinearData d = {x, y, 6};
+
+	FLT p[2] = {0, 0};
+	mp_result result = {0};
+
+	int status = mpfit(linear_well_conditioned, 6, 2, p, 0, 0, &d, &result);
+
+	if (status <= 0) {
+		fprintf(stderr, "WellConditionedFitStillConverges FAILED: mpfit status=%d\n", status);
+		return -1;
+	}
+	if (fabs(p[0] - 2.0) > 1e-4 || fabs(p[1] - 1.0) > 1e-4) {
+		fprintf(stderr, "WellConditionedFitStillConverges FAILED: p=[%.6f,%.6f], want [2,1]\n", p[0], p[1]);
+		return -1;
+	}
+	return 0;
+}