Add Minimizer::VarPro, a variable projection minimizer

Every bundled model now declares the parameters it is linear in
(linear_par, basis_and_grad) and, except the polynomials, a vectorisable
bulk basis pass (basis_columns) built on the new model::fast_exp. The
Levenberg-Marquardt iteration is factored into a shared DampedGaussNewton
driver with unchanged results; VarPro solves the linear parameters exactly
at every trial point and iterates over the nonlinear ones, rejecting
degenerate trial points and polishing the converged point with one
undamped step. Limits on linear parameters, pixels that do not converge
and models without the separable structure fall back to
LevenbergMarquardt. Parameter errors are the same Gauss-Newton estimates.
The [fit] tests run VarPro through every configuration and check the
separable structure of each model against its derivatives.

Measured on an M1 Max, single thread: FallingScurve population 5.5 us per
fit against 11.7 us with LevenbergMarquardt, and 2 wrong minima in 4096
pixels against 30; Gaussian 3.3 vs 4.3 us; GaussianChargeSharingKb 33 vs
49 us; Pol2 0.34 vs 1.52 us. Migrad stays the default.

Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>

src/Fit.cpp

@@ -3,6 +3,7 @@
 #include "FitHelpers.hpp"
 #include "FitMinuit2.hpp"
 #include "LevenbergMarquardt.hpp"
+#include "VariableProjection.hpp"
 #include "aare/Models.hpp"
 #include "aare/utils/par.hpp"
 #include "aare/utils/task.hpp"
@@ -10,6 +11,7 @@
 #include <cstddef>
 #include <stdexcept>
 #include <string>
+#include <type_traits>
 #include <vector>
 
 namespace aare {
@@ -117,11 +119,13 @@ std::vector<std::string> FitModel<Model>::GetParNames() const {
 
 namespace detail {
 
+struct NotSeparable {};
+
 /**
  * @brief Fits pixels with the minimizer selected by the model.
  *
- * Keep one instance per thread: the Levenberg-Marquardt buffers are reused
- * between pixels so that a data cube is fitted without per-pixel
+ * Keep one instance per thread: the buffers of the built-in minimizers are
+ * reused between pixels so that a data cube is fitted without per-pixel
  * allocations. The Minuit2 minimizers create their state per pixel.
  */
 template <typename Model> class PixelFitter {
@@ -139,18 +143,40 @@ template <typename Model> class PixelFitter {
              double *err_out, double &chi2) {
         const auto start = start_values(model, x, y);
 
-        if (model.minimizer() == Minimizer::LevenbergMarquardt) {
-            const auto res =
-                lm_.fit(model, x, y, y_err, start, par_out, err_out);
-            chi2 = res.valid ? res.chi2 : 0.0;
-            return res.valid;
+        if (model.minimizer() == Minimizer::VarPro) {
+            // A model without the separable structure, a pixel that does
+            // not converge and a pixel whose linear solution violates a
+            // limit fall back to the full solver.
+            if constexpr (model::is_separable<Model>::value) {
+                const auto res =
+                    vp_.fit(model, x, y, y_err, start, par_out, err_out);
+                if (res.valid) {
+                    chi2 = res.chi2;
+                    return true;
+                }
+            }
+            return fit_lm(model, x, y, y_err, start, par_out, err_out, chi2);
         }
+        if (model.minimizer() == Minimizer::LevenbergMarquardt)
+            return fit_lm(model, x, y, y_err, start, par_out, err_out, chi2);
         return fit_pixel_minuit2(model, x, y, y_err, start, par_out, err_out,
                                  chi2);
     }
 
   private:
+    bool fit_lm(const FitModel<Model> &model, NDView<double, 1> x,
+                NDView<double, 1> y, NDView<double, 1> y_err,
+                const std::array<double, Model::npar> &start, double *par_out,
+                double *err_out, double &chi2) {
+        const auto res = lm_.fit(model, x, y, y_err, start, par_out, err_out);
+        chi2 = res.valid ? res.chi2 : 0.0;
+        return res.valid;
+    }
+
     LevenbergMarquardt<Model> lm_;
+    std::conditional_t<model::is_separable<Model>::value,
+                       VariableProjection<Model>, NotSeparable>
+        vp_;
 };
 
 } // namespace detail