diff --git a/README.md b/README.md index 9e7a0a2..e6b1f05 100644 --- a/README.md +++ b/README.md @@ -225,6 +225,77 @@ The C++ optimizations provide significant speedups over pure Python: - SWIG interpolator and scipy's RegularGridInterpolator behave differently when coordinates aren't linspace - COBYLA in scipy can't set `rhobeg` per dimension individually, requiring problem scaling +### Reparameterized algorithms (4, 5, 6) and grain-bound saturation + +Algorithms 4-6 absorb the simplex constraint +(`f_sca + f_shade + f_bg = 1`, all ≥ 0) into a softmax reparameterization, +removing the need for inequality constraints. They differ in how they handle +the LUT box bounds on dust and grain: + +- **Algorithms 4, 5 (full softmax):** sigmoid reparameterization on dust and + grain, mapping the LUT range onto an unbounded variable in z-space. +- **Algorithm 6 (hybrid, recommended):** softmax for the fractions, but dust + and grain stay in physical units and are *clipped* to the LUT range inside + the objective — turning the bound into a true flat wall. + +The hybrid is the recommended path. On a real 50×50 Sentinel-2 patch (algorithm +benchmark, max_eval=100): + +| Algorithm | Median residual | Grain ≥ 1199 µm | Time | Speedup | +|-------------------------------|-----------------|-----------------|--------|------------| +| 1: COBYLA | 0.1013 | 0 / 2500 | 215 ms | 1.0× | +| 4: NELDERMEAD-softmax (full) | 0.0951 | 4 / 2500 | 94 ms | 2.3× | +| **6: NELDERMEAD-hybrid** | **0.0893** | 13 / 2500 | **84 ms** | **2.6×** | + +#### What grain-bound saturation actually means + +The three algorithms produce different saturation counts for *different reasons*: + +- **COBYLA (0/2500)** is implicitly regularized by the simplex inequality + constraint; its search structure pulls toward the simplex interior, masking + pixels whose true optimum lies at the LUT boundary. +- **Hybrid (13/2500 at max_eval=100)** finds the true boundary optima quickly: + the clip turns the bound into a flat wall, and the simplex contracts against + it in a few iterations and stops. We verified by direct comparison: at every + hybrid-saturated pixel, `grain ≈ 1200` produces a *lower* residual than + COBYLA's interior solution — these are genuine "grain is optimally large" + signals, not optimizer artifacts. +- **Full softmax (4/2500 at max_eval=100, 376/2500 at max_eval=500)** has both + a small set of genuine boundary cases *and* a drift mechanism: the sigmoid's + derivative `d_grain/d_z_grain` vanishes as `z_grain → ∞`, so the optimizer + keeps taking tiny improving steps that push `z_grain` upward without bound. + Raising `max_eval` doesn't approach a fixed point — it accumulates more + drift victims (~360 of them between 100 and 500 evals on this patch). + +The two saturation sets barely overlap: of the 4 softmax-saturated and 13 +hybrid-saturated pixels at max_eval=100, none are common. Different algorithms, +different signals. + +The diagnostic that separates "honest signal" from "drift artifact" is +**stability under max_eval**: + +| Algorithm | Saturation @ max_eval=100 | @ max_eval=500 | Δ | +|----------------------------|---------------------------|----------------------|-------------------| +| 4: NELDERMEAD-softmax | 4 / 2500 (0.16%) | **376 / 2500 (15%)** | **+15% (drift)** | +| 6: NELDERMEAD-hybrid | 13 / 2500 (0.5%) | 17 / 2500 (0.7%) | +0.2% (converged) | + +The hybrid's count is essentially constant; the full softmax's grows linearly +with max_eval. That growth is the failure mode — not the absolute count. + +The hybrid's clip-on-entry turns the LUT bound into a true flat plateau in the +objective: any value of dust or grain outside [min, max] maps to the same model +spectrum as the boundary itself, so further movement contributes nothing to the +residual and the simplex contracts and terminates. The clip introduces a +C^0-but-not-C^1 kink at the bound, which is benign for derivative-free solvers +(Nelder-Mead, COBYLA) but would need care for gradient-based methods. + +**Recommendation:** use algorithm 6 (`NELDERMEAD-hybrid`) for new work, and +treat its saturated pixels as "grain is at the LUT boundary" — flag them +downstream rather than assume they're optimizer noise. If using algorithms 4 +or 5, do not raise `max_eval` above the default of 100 without a downstream +filter, because the saturation count grows with max_eval rather than +converging. + ### Cross-platform numerical reproducibility Inversion results can differ by a few percent between Linux (x86_64, gcc) and macOS diff --git a/spires/invert.py b/spires/invert.py index 4b813c6..f9ae625 100755 --- a/spires/invert.py +++ b/spires/invert.py @@ -44,9 +44,26 @@ def speedy_invert(spectrum_target, spectrum_background, solar_angle, spectrum_sh Default is [0.5, 0.05, 10, 250]. algorithm : int, optional Optimization algorithm to use (default: 2). - 1 = LN_COBYLA (Constrained Optimization BY Linear Approximations), - 2 = LN_NELDERMEAD (Nelder-Mead simplex), - 3 = LD_SLSQP (Sequential Least Squares Programming, not working in C++). + 1 = LN_COBYLA (constrained, derivative-free), + 2 = LN_NELDERMEAD (unconstrained simplex; ignores box bounds), + 3 = LD_SLSQP (gradient-based; degraded — uses NLopt's finite-diff fallback), + 4 = LN_NELDERMEAD on softmax-reparameterized cost (full softmax), + 5 = LN_BOBYQA on softmax-reparameterized cost (quadratic-model variant), + 6 = LN_NELDERMEAD on hybrid: softmax for fractions, clip-on-entry for + dust/grain (recommended for real imagery). + Algorithms 4-6 absorb the simplex (f_sca + f_shade + f_bg = 1, all ≥ 0) + into the parameter transformation, so unconstrained NLopt solvers can + replace COBYLA. On real imagery the hybrid (algorithm 6) is the + recommended default replacement: it beats COBYLA on both fit quality + and speed (~2.6× faster) and stays stable as `max_eval` is raised. + + Note: algorithms 4 and 5 (full softmax) suffer grain-bound saturation + at high `max_eval` — the sigmoid reparameterization on dust/grain is + asymptotically flat near the LUT bounds, letting the optimizer drift + toward the upper bound while still lowering the residual. Algorithm 6 + does not have this problem because dust/grain stay in physical units + with a clip in the objective, turning the bound into a true wall. See + the "Softmax algorithms and grain-bound saturation" note in the README. Returns ------- @@ -134,9 +151,26 @@ def speedy_invert_array1d(spectra_targets, spectra_backgrounds, obs_solar_angles Default is [0.5, 0.05, 10, 250]. algorithm : int, optional Optimization algorithm to use (default: 2). - 1 = LN_COBYLA (Constrained Optimization BY Linear Approximations), - 2 = LN_NELDERMEAD (Nelder-Mead simplex), - 3 = LD_SLSQP (Sequential Least Squares Programming, not working in C++). + 1 = LN_COBYLA (constrained, derivative-free), + 2 = LN_NELDERMEAD (unconstrained simplex; ignores box bounds), + 3 = LD_SLSQP (gradient-based; degraded — uses NLopt's finite-diff fallback), + 4 = LN_NELDERMEAD on softmax-reparameterized cost (full softmax), + 5 = LN_BOBYQA on softmax-reparameterized cost (quadratic-model variant), + 6 = LN_NELDERMEAD on hybrid: softmax for fractions, clip-on-entry for + dust/grain (recommended for real imagery). + Algorithms 4-6 absorb the simplex (f_sca + f_shade + f_bg = 1, all ≥ 0) + into the parameter transformation, so unconstrained NLopt solvers can + replace COBYLA. On real imagery the hybrid (algorithm 6) is the + recommended default replacement: it beats COBYLA on both fit quality + and speed (~2.6× faster) and stays stable as `max_eval` is raised. + + Note: algorithms 4 and 5 (full softmax) suffer grain-bound saturation + at high `max_eval` — the sigmoid reparameterization on dust/grain is + asymptotically flat near the LUT bounds, letting the optimizer drift + toward the upper bound while still lowering the residual. Algorithm 6 + does not have this problem because dust/grain stay in physical units + with a clip in the objective, turning the bound into a true wall. See + the "Softmax algorithms and grain-bound saturation" note in the README. Returns ------- @@ -214,9 +248,26 @@ def speedy_invert_array2d(spectra_targets, spectra_backgrounds, obs_solar_angles Default is [0.5, 0.05, 10, 250]. algorithm : int, optional Optimization algorithm to use (default: 2). - 1 = LN_COBYLA (Constrained Optimization BY Linear Approximations), - 2 = LN_NELDERMEAD (Nelder-Mead simplex), - 3 = LD_SLSQP (Sequential Least Squares Programming, not working in C++). + 1 = LN_COBYLA (constrained, derivative-free), + 2 = LN_NELDERMEAD (unconstrained simplex; ignores box bounds), + 3 = LD_SLSQP (gradient-based; degraded — uses NLopt's finite-diff fallback), + 4 = LN_NELDERMEAD on softmax-reparameterized cost (full softmax), + 5 = LN_BOBYQA on softmax-reparameterized cost (quadratic-model variant), + 6 = LN_NELDERMEAD on hybrid: softmax for fractions, clip-on-entry for + dust/grain (recommended for real imagery). + Algorithms 4-6 absorb the simplex (f_sca + f_shade + f_bg = 1, all ≥ 0) + into the parameter transformation, so unconstrained NLopt solvers can + replace COBYLA. On real imagery the hybrid (algorithm 6) is the + recommended default replacement: it beats COBYLA on both fit quality + and speed (~2.6× faster) and stays stable as `max_eval` is raised. + + Note: algorithms 4 and 5 (full softmax) suffer grain-bound saturation + at high `max_eval` — the sigmoid reparameterization on dust/grain is + asymptotically flat near the LUT bounds, letting the optimizer drift + toward the upper bound while still lowering the residual. Algorithm 6 + does not have this problem because dust/grain stay in physical units + with a clip in the objective, turning the bound into a true wall. See + the "Softmax algorithms and grain-bound saturation" note in the README. bands : numpy.ndarray, optional Band wavelength coordinates of reflectances. Required if interpolator not provided. solar_angles : numpy.ndarray, optional @@ -308,9 +359,26 @@ def speedy_invert_xarray(spectra_targets, spectra_backgrounds, obs_solar_angles, Default is [0.5, 0.05, 10, 250]. algorithm : int, optional Optimization algorithm to use (default: 2). - 1 = LN_COBYLA (Constrained Optimization BY Linear Approximations), - 2 = LN_NELDERMEAD (Nelder-Mead simplex), - 3 = LD_SLSQP (Sequential Least Squares Programming, not working in C++). + 1 = LN_COBYLA (constrained, derivative-free), + 2 = LN_NELDERMEAD (unconstrained simplex; ignores box bounds), + 3 = LD_SLSQP (gradient-based; degraded — uses NLopt's finite-diff fallback), + 4 = LN_NELDERMEAD on softmax-reparameterized cost (full softmax), + 5 = LN_BOBYQA on softmax-reparameterized cost (quadratic-model variant), + 6 = LN_NELDERMEAD on hybrid: softmax for fractions, clip-on-entry for + dust/grain (recommended for real imagery). + Algorithms 4-6 absorb the simplex (f_sca + f_shade + f_bg = 1, all ≥ 0) + into the parameter transformation, so unconstrained NLopt solvers can + replace COBYLA. On real imagery the hybrid (algorithm 6) is the + recommended default replacement: it beats COBYLA on both fit quality + and speed (~2.6× faster) and stays stable as `max_eval` is raised. + + Note: algorithms 4 and 5 (full softmax) suffer grain-bound saturation + at high `max_eval` — the sigmoid reparameterization on dust/grain is + asymptotically flat near the LUT bounds, letting the optimizer drift + toward the upper bound while still lowering the residual. Algorithm 6 + does not have this problem because dust/grain stay in physical units + with a clip in the objective, turning the bound into a true wall. See + the "Softmax algorithms and grain-bound saturation" note in the README. Returns ------- @@ -609,6 +677,8 @@ def speedy_invert_scipy_softmax(interpolator: spires.interpolator.LutInterpolato return res, model_refl + + def speedy_invert_scipy(interpolator: spires.interpolator.LutInterpolator, spectrum_target, spectrum_background, solar_angle, shade=None, scipy_options=None, mode=3, method='SLSQP'): diff --git a/spires/spires.cpp b/spires/spires.cpp index 3d2fbe5..6ecf8bd 100755 --- a/spires/spires.cpp +++ b/spires/spires.cpp @@ -224,6 +224,205 @@ double spectrum_difference(const std::vector& x, } +// ---------------------------------------------------------------------------- +// Softmax reparameterization +// +// Maps an unconstrained vector z ∈ R^4 to physical parameters such that the +// simplex constraints (f_sca, f_shade, f_bg ≥ 0; f_sca + f_shade + f_bg = 1) +// and the box bounds on dust/grain are satisfied by construction. This lets +// unconstrained solvers (Nelder-Mead, BOBYQA) replace COBYLA on the fractional +// sub-problem. +// +// z = [z_snow, z_shade, z_dust, z_grain] +// (f_sca, f_shade, f_bg) = softmax(z_snow, z_shade, 0) // z_bg pinned to 0 +// dust = dust_min + (dust_max - dust_min) * sigmoid(z_dust) +// grain = grain_min + (grain_max - grain_min) * sigmoid(z_grain) +// +// `x_to_z` is used to seed z0 from a user-supplied physical x0 (e.g., the +// default [0.5, 0.05, 10, 250]). `z_to_x` is applied to the optimized z* before +// returning. Both are stable under (eps clamp on) the input range. +// ---------------------------------------------------------------------------- + +static inline void softmax_to_fractions(double z_snow, double z_shade, + double& f_sca, double& f_shade, double& f_bg) { + // Subtract the max for numerical stability — exp(x - max) ∈ (0, 1]. + double m = std::max({z_snow, z_shade, 0.0}); + double e0 = std::exp(z_snow - m); + double e1 = std::exp(z_shade - m); + double e2 = std::exp(-m); + double s = e0 + e1 + e2; + f_sca = e0 / s; + f_shade = e1 / s; + f_bg = e2 / s; +} + + +static inline double sigmoid_to_bounded(double z, double lo, double hi) { + return lo + (hi - lo) / (1.0 + std::exp(-z)); +} + + +std::vector z_to_x(const std::vector& z, + double dust_min, double dust_max, + double grain_min, double grain_max) { + double f_sca, f_shade, f_bg; + softmax_to_fractions(z[0], z[1], f_sca, f_shade, f_bg); + double dust = sigmoid_to_bounded(z[2], dust_min, dust_max); + double grain = sigmoid_to_bounded(z[3], grain_min, grain_max); + return {f_sca, f_shade, dust, grain}; +} + + +std::vector x_to_z(const std::vector& x, + double dust_min, double dust_max, + double grain_min, double grain_max) { + // Inverse of z_to_x. Used to translate a user-supplied physical x0 into + // initial guess in z-space. `eps` keeps logit() / log() finite when the + // user passes 0 or 1 exactly (or x outside the LUT bounds). + constexpr double eps = 1e-6; + auto clamp = [](double v, double lo, double hi) { + return std::min(std::max(v, lo), hi); + }; + + double f_sca = clamp(x[0], eps, 1.0 - eps); + double f_shade = clamp(x[1], eps, 1.0 - eps); + double f_bg = clamp(1.0 - f_sca - f_shade, eps, 1.0 - eps); + // Renormalize after clamping so log-ratios are sensible. + double s = f_sca + f_shade + f_bg; + f_sca /= s; f_shade /= s; f_bg /= s; + double z_snow = std::log(f_sca / f_bg); + double z_shade = std::log(f_shade / f_bg); + + double u_d = clamp((x[2] - dust_min) / (dust_max - dust_min), eps, 1.0 - eps); + double u_g = clamp((x[3] - grain_min) / (grain_max - grain_min), eps, 1.0 - eps); + double z_dust = std::log(u_d / (1.0 - u_d)); + double z_grain = std::log(u_g / (1.0 - u_g)); + + return {z_snow, z_shade, z_dust, z_grain}; +} + + +double spectrum_difference_hybrid(const std::vector& y, + double* spectrum_background, int len_background, + double* spectrum_target, int len_target, + double* spectrum_shade, int len_shade, + double solar_angle, + double* bands, int len_bands, + double* solar_angles, int len_solar_angles, + double* dust_concentrations, int len_dust_concentrations, + double* grain_sizes, int len_grain_sizes, + double* lut, int n_bands, int n_solar_angles, int n_dust_concentrations, int n_grain_sizes) { + /* + Hybrid spectral-difference cost: + - Fractions (f_sca, f_shade, f_bg) via softmax of (y[0], y[1], 0). + - Dust and grain in physical units, CLIPPED to the LUT range on entry. + + The clip turns the LUT-bound into a true wall in the objective: any value + of y[2] / y[3] outside [min, max] is mapped to the same model spectrum as + the boundary itself, so the optimizer sees a flat plateau and stops. This + avoids the asymptotic drift that the full sigmoid reparameterization + suffers (see README "Softmax algorithms and grain-bound saturation"). + + The clip introduces a C^0-but-not-C^1 kink at the bound, which is benign + for derivative-free solvers (Nelder-Mead, COBYLA) but would need care for + gradient-based ones. + + y = [z_snow, z_shade, dust, grain] + */ + double f_sca, f_shade, f_bg; + softmax_to_fractions(y[0], y[1], f_sca, f_shade, f_bg); + + double dust_min = dust_concentrations[0]; + double dust_max = dust_concentrations[len_dust_concentrations - 1]; + double grain_min = grain_sizes[0]; + double grain_max = grain_sizes[len_grain_sizes - 1]; + double dust = std::min(std::max(y[2], dust_min), dust_max); + double grain = std::min(std::max(y[3], grain_min), grain_max); + + double solar_angle_idx = get_idx(solar_angle, solar_angles, len_solar_angles); + double dust_idx = get_idx(dust, dust_concentrations, len_dust_concentrations); + double grain_idx = get_idx(grain, grain_sizes, len_grain_sizes); + + double diff_sq = 0.0; + for (int i = 0; i < len_target; ++i) { + double model_pure = interpolate_idx(lut, n_bands, n_solar_angles, n_dust_concentrations, n_grain_sizes, + i, solar_angle_idx, dust_idx, grain_idx); + double model = model_pure * f_sca + spectrum_shade[i] * f_shade + spectrum_background[i] * f_bg; + double d = spectrum_target[i] - model; + diff_sq += d * d; + } + return std::sqrt(diff_sq); +} + + +// Helpers for translating between user-facing physical x = [f_sca, f_shade, dust, grain] +// and the hybrid's internal y = [z_snow, z_shade, dust, grain]. Dust and grain +// are passed through unchanged; only the fractions are reparameterized. +std::vector y_to_x_hybrid(const std::vector& y) { + double f_sca, f_shade, f_bg; + softmax_to_fractions(y[0], y[1], f_sca, f_shade, f_bg); + return {f_sca, f_shade, y[2], y[3]}; +} + + +std::vector x_to_y_hybrid(const std::vector& x) { + constexpr double eps = 1e-6; + auto clamp = [](double v, double lo, double hi) { + return std::min(std::max(v, lo), hi); + }; + double f_sca = clamp(x[0], eps, 1.0 - eps); + double f_shade = clamp(x[1], eps, 1.0 - eps); + double f_bg = clamp(1.0 - f_sca - f_shade, eps, 1.0 - eps); + double s = f_sca + f_shade + f_bg; + f_sca /= s; f_shade /= s; f_bg /= s; + double z_snow = std::log(f_sca / f_bg); + double z_shade = std::log(f_shade / f_bg); + return {z_snow, z_shade, x[2], x[3]}; +} + + +double spectrum_difference_softmax(const std::vector& z, + double* spectrum_background, int len_background, + double* spectrum_target, int len_target, + double* spectrum_shade, int len_shade, + double solar_angle, + double* bands, int len_bands, + double* solar_angles, int len_solar_angles, + double* dust_concentrations, int len_dust_concentrations, + double* grain_sizes, int len_grain_sizes, + double* lut, int n_bands, int n_solar_angles, int n_dust_concentrations, int n_grain_sizes) { + /* + Spectral-difference cost in unconstrained (softmax-reparameterized) + coordinates. Unconstrained solvers (Nelder-Mead, BOBYQA, BFGS) can be used + directly because the simplex and box constraints are absorbed into the + transformation z -> x. + */ + double f_sca, f_shade, f_bg; + softmax_to_fractions(z[0], z[1], f_sca, f_shade, f_bg); + + double dust_min = dust_concentrations[0]; + double dust_max = dust_concentrations[len_dust_concentrations - 1]; + double grain_min = grain_sizes[0]; + double grain_max = grain_sizes[len_grain_sizes - 1]; + double dust = sigmoid_to_bounded(z[2], dust_min, dust_max); + double grain = sigmoid_to_bounded(z[3], grain_min, grain_max); + + double solar_angle_idx = get_idx(solar_angle, solar_angles, len_solar_angles); + double dust_idx = get_idx(dust, dust_concentrations, len_dust_concentrations); + double grain_idx = get_idx(grain, grain_sizes, len_grain_sizes); + + double diff_sq = 0.0; + for (int i = 0; i < len_target; ++i) { + double model_pure = interpolate_idx(lut, n_bands, n_solar_angles, n_dust_concentrations, n_grain_sizes, + i, solar_angle_idx, dust_idx, grain_idx); + double model = model_pure * f_sca + spectrum_shade[i] * f_shade + spectrum_background[i] * f_bg; + double d = spectrum_target[i] - model; + diff_sq += d * d; + } + return std::sqrt(diff_sq); +} + + double index_to_value(double value, double* coords, int len_coords) { // value in [0, 1] mapped to coords[0] ... coords[len-1]. // Previous version used `idx = value * len_coords` which read past the end @@ -303,6 +502,36 @@ static double spectrum_difference_wrapper(const std::vector& x, std::vec } +static double spectrum_difference_softmax_wrapper(const std::vector& z, std::vector& /*grad*/, void* data) { + ObjectiveData* d = reinterpret_cast(data); + return spectrum_difference_softmax(z, + d->spectrum_background, d->len_background, + d->spectrum_target, d->len_target, + d->spectrum_shade, d->len_shade, + d->solar_angle, + d->bands, d->len_bands, + d->solar_angles, d->len_solar_angles, + d->dust_concentrations, d->len_dust_concentrations, + d->grain_sizes, d->len_grain_sizes, + d->lut, d->n_bands, d->n_solar_angles, d->n_dust_concentrations, d->n_grain_sizes); +} + + +static double spectrum_difference_hybrid_wrapper(const std::vector& y, std::vector& /*grad*/, void* data) { + ObjectiveData* d = reinterpret_cast(data); + return spectrum_difference_hybrid(y, + d->spectrum_background, d->len_background, + d->spectrum_target, d->len_target, + d->spectrum_shade, d->len_shade, + d->solar_angle, + d->bands, d->len_bands, + d->solar_angles, d->len_solar_angles, + d->dust_concentrations, d->len_dust_concentrations, + d->grain_sizes, d->len_grain_sizes, + d->lut, d->n_bands, d->n_solar_angles, d->n_dust_concentrations, d->n_grain_sizes); +} + + static double constraint(const std::vector& x, std::vector& /*grad*/, void* /*data*/) { // f_sca + f_shade <= 1, written as f_sca + f_shade - 1 <= 0 return x[0] + x[1] - 1; @@ -347,7 +576,9 @@ std::vector invert(double* spectrum_background, int len_background, }; nlopt::opt opt; - bool constrained_algorithm; + bool constrained_algorithm = false; + bool softmax_algorithm = false; + bool hybrid_algorithm = false; switch (algorithm) { case 1: { @@ -370,21 +601,99 @@ std::vector invert(double* spectrum_background, int len_background, opt = nlopt::opt(nlopt::LD_SLSQP, 4); constrained_algorithm = true; break; + case 4: + // LN_NELDERMEAD on softmax-reparameterized problem. Unconstrained: + // no bounds, no inequality constraint — simplex and box constraints + // are absorbed into the z -> x transformation. + opt = nlopt::opt(nlopt::LN_NELDERMEAD, 4); + softmax_algorithm = true; + break; + case 5: + // LN_BOBYQA on softmax. Often faster than Nelder-Mead on smooth + // objectives — uses a quadratic model rather than simplex reflections. + opt = nlopt::opt(nlopt::LN_BOBYQA, 4); + softmax_algorithm = true; + break; + case 6: + // LN_NELDERMEAD on hybrid: softmax for fractions, clip-on-entry for + // dust/grain. The clip turns the LUT bounds into a true wall in the + // objective (flat plateau outside), avoiding the asymptotic drift that + // the full sigmoid suffers. Recommended for use on real imagery. + opt = nlopt::opt(nlopt::LN_NELDERMEAD, 4); + hybrid_algorithm = true; + break; default: - throw std::invalid_argument("invert: unknown algorithm code (must be 1=COBYLA, 2=NELDERMEAD, 3=SLSQP)"); + throw std::invalid_argument("invert: unknown algorithm code " + "(1=COBYLA, 2=NELDERMEAD, 3=SLSQP, " + "4=NELDERMEAD-softmax, 5=BOBYQA-softmax, 6=NELDERMEAD-hybrid)"); } + double min_dust_concentration = dust_concentrations[0]; + double max_dust_concentration = dust_concentrations[len_dust_concentrations - 1]; + double min_grain_size = grain_sizes[0]; + double max_grain_size = grain_sizes[len_grain_sizes - 1]; + + if (softmax_algorithm) { + // Unconstrained path: objective is the softmax cost; no bounds; seed + // from x0 by inverting the reparameterization. + // + // Tolerances are set in z-space, which has different scale from x-space: + // - Each unit step in z_snow / z_shade roughly doubles/halves a softmax + // fraction; each unit in z_dust / z_grain moves the bounded variable + // through ~1/4 of its range. So an initial step of ~0.5 is roughly + // equivalent to COBYLA's rhobeg=(0.1, 0.1, 100, 100) in x-space. + // - xtol_abs of 1e-3 in z corresponds to ~0.05% movement in any + // fraction or bounded var. Tighter than the constrained path's + // xtol_rel=1e-2 because BOBYQA otherwise stops prematurely. + opt.set_min_objective(spectrum_difference_softmax_wrapper, &obj_data); + opt.set_maxeval(max_eval); + opt.set_ftol_abs(1e-6); + opt.set_xtol_abs(1e-3); + opt.set_initial_step(std::vector{0.5, 0.5, 0.5, 0.5}); + + std::vector z = x_to_z(x0, + min_dust_concentration, max_dust_concentration, + min_grain_size, max_grain_size); + double minf; + opt.optimize(z, minf); + return z_to_x(z, + min_dust_concentration, max_dust_concentration, + min_grain_size, max_grain_size); + } + + if (hybrid_algorithm) { + // Hybrid path: softmax-reparameterized fractions (z_snow, z_shade) but + // dust/grain stay in physical units, clipped inside the objective. No + // NLopt bounds — the clip turns the LUT bound into a flat plateau + // that the simplex contracts against, terminating cleanly. + // + // Initial step: z-space scale (~0.5) for fractions, physical units + // (~100) for dust/grain — matches COBYLA's rhobeg on the latter two. + opt.set_min_objective(spectrum_difference_hybrid_wrapper, &obj_data); + opt.set_maxeval(max_eval); + opt.set_ftol_abs(1e-4); + opt.set_xtol_rel(1e-2); + opt.set_initial_step(std::vector{0.5, 0.5, 100.0, 100.0}); + + std::vector y = x_to_y_hybrid(x0); + double minf; + opt.optimize(y, minf); + std::vector x = y_to_x_hybrid(y); + // Clamp dust/grain to LUT range on output — the optimizer's working + // value can be outside if the plateau pushed it there; the user- + // facing answer must stay physical. + x[2] = std::min(std::max(x[2], min_dust_concentration), max_dust_concentration); + x[3] = std::min(std::max(x[3], min_grain_size), max_grain_size); + return x; + } + + // Constrained / bounded path (algorithms 1-3). if (constrained_algorithm) { opt.add_inequality_constraint(constraint, &obj_data); } opt.set_min_objective(spectrum_difference_wrapper, &obj_data); opt.set_maxeval(max_eval); - double min_dust_concentration = dust_concentrations[0]; - double max_dust_concentration = dust_concentrations[len_dust_concentrations - 1]; - double min_grain_size = grain_sizes[0]; - double max_grain_size = grain_sizes[len_grain_sizes - 1]; - std::vector lower_bounds = {0.0, 0.0, min_dust_concentration, min_grain_size}; std::vector upper_bounds = {1.0, 1.0, max_dust_concentration, max_grain_size}; opt.set_lower_bounds(lower_bounds); diff --git a/spires/spires.h b/spires/spires.h index 13c9e28..438d81d 100755 --- a/spires/spires.h +++ b/spires/spires.h @@ -62,6 +62,44 @@ double spectrum_difference_scaled(const std::vector &x, double* lut, int n_bands, int n_solar_angles, int n_dust_concentrations, int n_grain_sizes); +double spectrum_difference_softmax(const std::vector& z, + double* spectrum_background, int len_background, + double* spectrum_target, int len_target, + double* spectrum_shade, int len_shade, + double solar_angle, + double* bands, int len_bands, + double* solar_angles, int len_solar_angles, + double* dust_concentrations, int len_dust_concentrations, + double* grain_sizes, int len_grain_sizes, + double* lut, int n_bands, int n_solar_angles, int n_dust_concentrations, int n_grain_sizes); + + +double spectrum_difference_hybrid(const std::vector& y, + double* spectrum_background, int len_background, + double* spectrum_target, int len_target, + double* spectrum_shade, int len_shade, + double solar_angle, + double* bands, int len_bands, + double* solar_angles, int len_solar_angles, + double* dust_concentrations, int len_dust_concentrations, + double* grain_sizes, int len_grain_sizes, + double* lut, int n_bands, int n_solar_angles, int n_dust_concentrations, int n_grain_sizes); + + +std::vector z_to_x(const std::vector& z, + double dust_min, double dust_max, + double grain_min, double grain_max); + + +std::vector x_to_z(const std::vector& x, + double dust_min, double dust_max, + double grain_min, double grain_max); + + +std::vector y_to_x_hybrid(const std::vector& y); +std::vector x_to_y_hybrid(const std::vector& x); + + std::vector invert(double* spectrum_background, int len_background, double* spectrum_target, int len_target, double* spectrum_shade, int len_shade, diff --git a/tests/test_swig.py b/tests/test_swig.py index 4aacdcc..faa493e 100755 --- a/tests/test_swig.py +++ b/tests/test_swig.py @@ -183,3 +183,180 @@ def test_interpolate_all_handles_non_9_band_lut(): dust_concentration=dust_concentration, grain_size=grain_size) assert np.asarray(ret).shape == (n_bands,) + + +def _residual(x, target, background): + return spires.core.spectrum_difference( + x=list(x), + spectrum_background=background, + spectrum_target=target, + spectrum_shade=np.zeros_like(target), + solar_angle=solar_angle, + bands=interpolator.bands, + solar_angles=interpolator.solar_angles, + dust_concentrations=interpolator.dust_concentrations, + grain_sizes=interpolator.grain_sizes, + lut=interpolator.reflectances) + + +@pytest.mark.parametrize("algorithm,name", [(4, "NELDERMEAD-softmax"), + (5, "BOBYQA-softmax"), + (6, "NELDERMEAD-hybrid")]) +def test_invert_softmax(algorithm, name): + """Softmax/hybrid algorithms (4, 5, 6) absorb the simplex constraints into + the parameter transformation, so unconstrained NLopt solvers can be used. + Asserted on physical plausibility and residual fit quality (no pinned + coordinates — they drift across platforms and the softmax surface is + flatter near the bounds than the constrained version).""" + x = spires.core.invert(spectrum_background=spectrum_background, + spectrum_target=spectrum_target, + spectrum_shade=spectrum_shade, + solar_angle=solar_angle, + bands=interpolator.bands, + solar_angles=interpolator.solar_angles, + dust_concentrations=interpolator.dust_concentrations, + grain_sizes=interpolator.grain_sizes, + lut=interpolator.reflectances, + max_eval=500, + x0=x0, + algorithm=algorithm) + x = np.asarray(x) + + # Physical plausibility: simplex + box bounds satisfied by construction. + assert 0 <= x[0] <= 1, f"{name}: f_sca {x[0]} out of [0,1]" + assert 0 <= x[1] <= 1, f"{name}: f_shade {x[1]} out of [0,1]" + assert x[0] + x[1] <= 1 + 1e-6, f"{name}: f_sca + f_shade > 1" + assert (interpolator.dust_concentrations.min() <= x[2] <= + interpolator.dust_concentrations.max()), f"{name}: dust {x[2]} out of LUT range" + assert (interpolator.grain_sizes.min() <= x[3] <= + interpolator.grain_sizes.max()), f"{name}: grain {x[3]} out of LUT range" + + # Residual: softmax variants should match or beat COBYLA's fit on this + # pixel (Python A/B confirmed the constraint was holding COBYLA back). + residual_softmax = _residual(x, spectrum_target, spectrum_background) + x_cobyla = spires.core.invert(spectrum_background=spectrum_background, + spectrum_target=spectrum_target, + spectrum_shade=spectrum_shade, + solar_angle=solar_angle, + bands=interpolator.bands, + solar_angles=interpolator.solar_angles, + dust_concentrations=interpolator.dust_concentrations, + grain_sizes=interpolator.grain_sizes, + lut=interpolator.reflectances, + max_eval=500, x0=x0, algorithm=1) + residual_cobyla = _residual(np.asarray(x_cobyla), spectrum_target, spectrum_background) + assert residual_softmax <= residual_cobyla * 1.05, ( + f"{name}: residual {residual_softmax:.5f} >> COBYLA's {residual_cobyla:.5f}") + + +def _real_imagery_setup(): + """Load the Sentinel-2 LFS subset, or skip the test cleanly. Returns + (R, R0, sza_array, sza_scalar, residual_fn, n_pixels).""" + xr = pytest.importorskip("xarray") + try: + ds = xr.open_dataset('tests/data/sentinel_r_subset.nc') + ds0 = xr.open_dataset('tests/data/sentinel_r0_subset.nc') + except (OSError, ValueError): + pytest.skip("LFS test data not available (run `git lfs pull`)") + + R = np.ascontiguousarray( + ds['reflectance'].isel(time=0).transpose('y', 'x', 'band').values.astype(np.float64)) + R0 = np.ascontiguousarray( + ds0['reflectance'].transpose('y', 'x', 'band').values.astype(np.float64)) + sza = np.full(R.shape[:2], + float(np.nanmean(ds['sun_zenith_grid'].isel(time=0).values))) + sza0 = float(sza[0, 0]) + shade = np.zeros(9) + n = R.shape[0] * R.shape[1] + + def residuals(res): + flat = res.reshape(-1, 4) + flat_t = R.reshape(-1, 9) + flat_b = R0.reshape(-1, 9) + return np.array([ + spires.snow_diff_4(flat[i], flat_t[i], flat_b[i], sza0, interpolator, shade) + for i in range(n) + ]) + + return R, R0, sza, sza0, residuals, n + + +def test_softmax_beats_cobyla_on_real_imagery(): + """End-to-end contract on a real 50x50 Sentinel-2 patch (LFS): the softmax + Nelder-Mead variant (algorithm=4) must produce a median residual that is at + least as good as COBYLA's at the same max_eval, and must not saturate the + grain bound on more than ~2% of pixels at max_eval=100.""" + R, R0, sza, _, residuals, n = _real_imagery_setup() + + res_cobyla = spires.speedy_invert_array2d( + spectra_targets=R, spectra_backgrounds=R0, obs_solar_angles=sza, + interpolator=interpolator, algorithm=1, max_eval=100, x0=np.array(x0)) + res_softmax = spires.speedy_invert_array2d( + spectra_targets=R, spectra_backgrounds=R0, obs_solar_angles=sza, + interpolator=interpolator, algorithm=4, max_eval=100, x0=np.array(x0)) + + r_cobyla = residuals(res_cobyla) + r_softmax = residuals(res_softmax) + + assert np.median(r_softmax) <= np.median(r_cobyla), ( + f"softmax median {np.median(r_softmax):.4f} > COBYLA median {np.median(r_cobyla):.4f}") + + grain_max = interpolator.grain_sizes.max() + grain = res_softmax[..., 3] + n_saturated = int((grain >= grain_max - 1).sum()) + assert n_saturated / n <= 0.02, ( + f"softmax saturated grain on {n_saturated}/{n} pixels (>2%) — " + "the unconstrained problem may be too loose at this max_eval") + + +def test_hybrid_saturation_is_stable_under_max_eval(): + """The hybrid algorithm (6 = softmax for fractions + clip for dust/grain) + should converge to a stable saturation set: pixels whose true optimum lies + at the LUT grain boundary find that boundary in a few iterations and stop, + so the saturation count at max_eval=500 must be close to the count at + max_eval=100. + + The full sigmoid (algorithm 4) lacks this property — its saturation count + grows roughly linearly with max_eval as more pixels drift up the asymptotic + z-ridge. So this test pins the *stability* property that distinguishes + honest boundary signal from optimizer drift. On this patch the hybrid + grows by ~4 pixels (13 → 17) between max_eval 100 and 500; the full + softmax grows by ~370 (4 → 376).""" + R, R0, sza, _, _, n = _real_imagery_setup() + grain_max = interpolator.grain_sizes.max() + + def n_saturated(max_eval): + res = spires.speedy_invert_array2d( + spectra_targets=R, spectra_backgrounds=R0, obs_solar_angles=sza, + interpolator=interpolator, algorithm=6, max_eval=max_eval, + x0=np.array(x0)) + return int((res[..., 3] >= grain_max - 1).sum()) + + sat_100 = n_saturated(100) + sat_500 = n_saturated(500) + growth = (sat_500 - sat_100) / n + + # Stability bar: ≤1% additional pixels saturate when max_eval grows 5×. + # Hybrid currently shows ~0.2% on this patch; full softmax shows ~15%. + # The 1% bar separates the regimes with platform-variance headroom. + assert growth <= 0.01, ( + f"hybrid saturation grew from {sat_100} to {sat_500} of {n} pixels " + f"({growth:.2%}) when max_eval went 100 -> 500. Growth >1% suggests " + "the clip-on-entry mechanism is no longer pinning saturated pixels — " + "the optimizer may be drifting like the full softmax (algorithm 4)") + + +def test_invert_unknown_algorithm_raises(): + """Unknown algorithm codes must raise rather than silently running with an + uninitialized nlopt::opt — guards against the previous fallthrough bug.""" + with pytest.raises(RuntimeError): + spires.core.invert(spectrum_background=spectrum_background, + spectrum_target=spectrum_target, + spectrum_shade=spectrum_shade, + solar_angle=solar_angle, + bands=interpolator.bands, + solar_angles=interpolator.solar_angles, + dust_concentrations=interpolator.dust_concentrations, + grain_sizes=interpolator.grain_sizes, + lut=interpolator.reflectances, + max_eval=100, x0=x0, algorithm=99)