diff --git a/RELEASE.md b/RELEASE.md index a5d3284b..479c3e48 100644 --- a/RELEASE.md +++ b/RELEASE.md @@ -36,7 +36,10 @@ analytic model derivatives and needs far fewer function evaluations. ``Minimizer.LevenbergMarquardt`` is a built-in, dependency-free Levenberg-Marquardt solver with the analytic Jacobian; it reflects steps at - parameter limits and fits data cubes without per-pixel allocations. With + parameter limits and fits data cubes without per-pixel allocations. A fit + that stalls away from a minimum, or whose width collapses below the spacing + of the scan points, is reported as failed (all zeros) rather than returned + as a result, so noisy or signal-free pixels may fail. With ``compute_errors``, Fumili and LevenbergMarquardt report parameter errors from their linearised covariance instead of running Hesse. Pixels for which Fumili does not reach a valid minimum, which Minuit2's implementation cannot @@ -56,7 +59,10 @@ derivatives by the tests). Each trial point solves those exactly and the Levenberg-Marquardt iteration runs over the nonlinear parameters only, so the fit needs fewer evaluations and does not depend on start values of the - linear parameters. Limits on linear parameters, pixels that do not converge + linear parameters. When the basis functions are nearly collinear, as for a + polynomial on a narrow range far from x = 0, one more pass over the data + refines the linear solution to the precision of the other minimizers. + 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 Gaussian, plateau, charge-sharing and S-curve models also @@ -141,6 +147,14 @@ ### Bugfixes: - Fitting a 3D data cube in Python without ``y_err`` now returns ``par_err`` when ``compute_errors`` is set, as fitting a single pixel already did. +- Samples with a zero error are ignored by every minimizer, as documented, + also when they hold NaN. They are now removed before the start values are + estimated; previously a masked NaN at a sample used by the estimate, such + as the first or last point of ``Pol1``, made even Migrad fail. A pixel with + fewer remaining samples than free parameters is reported as failed. +- Fitting a single pixel checks that ``x``, ``y`` and ``y_err`` have the same + size and raises ``RuntimeError`` otherwise. Previously a shorter ``y`` or + ``y_err`` was read past its end. - The Python ``Cluster`` constructors validate that the data array holds exactly one value per pixel and raise ``ValueError`` otherwise. Previously a longer array wrote past the end of the cluster data. diff --git a/docs/src/pyFit.rst b/docs/src/pyFit.rst index 1d5d921e..142b1ca6 100644 --- a/docs/src/pyFit.rst +++ b/docs/src/pyFit.rst @@ -130,7 +130,11 @@ control the minimizer: - ``tolerance`` is the EDM tolerance in Minuit's convention. Migrad, LevenbergMarquardt and VarPro stop when the estimated distance to the minimum is below ``0.002 * tolerance``, Fumili below - ``1e-4 * tolerance``. + ``1e-4 * tolerance``. The EDM is measured in units of ``chi2``, not + relative to the data: without ``y_err`` every sample counts with an error + of 1, so a fit of small values (amplitudes well below 1) meets the + tolerance while still far from the minimum. Pass ``y_err`` or lower + ``tolerance`` for such data. - ``strategy`` is the Minuit2 strategy (0 is fast, 1 is Minuit2's default). LevenbergMarquardt and VarPro ignore it. - ``compute_errors`` adds ``par_err`` to the result. Fixed parameters and diff --git a/include/aare/FitModel.hpp b/include/aare/FitModel.hpp index 6b2e0ceb..f485849e 100644 --- a/include/aare/FitModel.hpp +++ b/include/aare/FitModel.hpp @@ -43,7 +43,9 @@ namespace aare { * report an error of 0 and a failed fit returns zeros. `tolerance` is the * EDM tolerance in Minuit's convention (Migrad, LevenbergMarquardt and * VarPro stop below 0.002 * tolerance, Fumili below - * 1e-4 * tolerance). `max_calls` bounds the work per pixel in the + * 1e-4 * tolerance). The EDM is in units of chi2, so an unweighted fit of + * small values meets it far from the minimum; pass errors or a lower + * tolerance for such data. `max_calls` bounds the work per pixel in the * minimizer's own units (Minuit2 function calls, or model evaluations for * LevenbergMarquardt and VarPro, which spend one evaluation per * iteration and two when a step crosses a limit). `strategy` only affects diff --git a/python/tests/test_Fit.py b/python/tests/test_Fit.py index a5436b60..2e876ed9 100644 --- a/python/tests/test_Fit.py +++ b/python/tests/test_Fit.py @@ -175,6 +175,15 @@ def test_invalid_parameter_access_raises(): model.SetParLimits("mu", 1.0, 1.0) +def test_mismatched_array_sizes_raise(gaussian_data, minimizer): + x, y, y_err = gaussian_data + model = aare.Gaussian(minimizer=minimizer) + with pytest.raises(RuntimeError): + model.fit(x, y[:-1]) + with pytest.raises(RuntimeError): + model.fit(x, y, y_err[:-1]) + + @pytest.mark.parametrize("cls", [aare.RisingScurve, aare.FallingScurve]) def test_scurve_fit(cls, minimizer): truth = np.array([10.0, 0.1, 50.0, 3.0, 1000.0, 2.0]) @@ -239,3 +248,154 @@ def test_all_minimizers_agree_on_3d_data_with_errors(): np.testing.assert_allclose(result["par"], np.broadcast_to(expected, (2, 3, 3)), atol=2e-3) np.testing.assert_allclose(result["chi2"], 0.0, atol=1e-3) np.testing.assert_allclose(result["par_err"], reference["par_err"], rtol=2e-2) + + +# Regressions of the built-in solvers: fits that used to be reported as valid +# although they had stopped away from the minimum or returned garbage. + + +def test_gaussian_fit_to_noise_fails_or_stays_bounded(minimizer): + # A width collapsing between two scan points let the amplitude diverge + # (up to 1e159) in a fit that was still reported as valid. + x = np.linspace(-6.0, 6.0, 61) + y = np.random.default_rng(0).normal(0.0, 1.0, (20, 20, x.size)) + result = aare.Gaussian(minimizer=minimizer).fit(x, y, n_threads=2) + + par = result["par"] + failed = np.all(par == 0.0, axis=-1) + bounded = np.all(np.isfinite(par) & (np.abs(par) < 1e6), axis=-1) + assert np.all(failed | bounded) + assert np.all(np.isfinite(result["chi2"])) + + +def test_step_reflected_onto_start_is_not_convergence(minimizer): + # The step of +2 reflects at the limit 1 back onto the start 0; the + # truncated step must still be tried. + model = aare.Pol1(minimizer=minimizer) + model.SetParameter("p0", 0.0) + model.FixParameter("p1", 0.0) + model.SetParLimits("p0", -1.0, 1.0) + result = model.fit(np.arange(4.0), np.full(4, 2.2)) + + assert result["par"][0] == pytest.approx(1.0, abs=1e-6) + assert result["chi2"][0] == pytest.approx(5.76, rel=1e-6) + + +@pytest.mark.parametrize( + "amplitude, offset", [(1e8, 0.0), (1e10, 0.0), (120.0, 1e9)] +) +def test_large_parameter_does_not_stop_the_fit(minimizer, amplitude, offset): + # One relative step norm over all parameters let a large amplitude or + # position stop the fit at the start estimate of the width. + x = np.linspace(-6.0, 6.0, 61) + y = gaussian(x, amplitude, 0.8, 1.3) + y_err = np.sqrt(y) + 1.0 + result = aare.Gaussian(minimizer=minimizer, max_calls=500).fit( + x + offset, y, y_err + ) + + par = result["par"] + assert par[0] == pytest.approx(amplitude, rel=1e-3) + assert par[1] - offset == pytest.approx(0.8, abs=1e-3) + assert par[2] == pytest.approx(1.3, rel=1e-3) + + +def test_rejected_trial_on_limit_does_not_freeze_interior_parameter(minimizer): + # A rejected trial truncated onto the lower limit froze mu at its start. + x = np.linspace(-0.1, 0.1, 11) + mu = np.sqrt(2.0 * np.log(2.0)) + y = gaussian(x, 10.0, mu, 1.0) + model = aare.Gaussian(minimizer=minimizer) + model.FixParameter("A", 10.0) + model.FixParameter("sigma", 1.0) + model.SetParameter("mu", 3.0) + model.SetParLimits("mu", 0.0, 5.0) + result = model.fit(x, y) + + assert result["par"][1] == pytest.approx(mu, abs=1e-3) + assert result["chi2"][0] < 1e-4 + + +# Samples with a zero error are ignored, also when they hold NaN (a masked +# dead frame): a zero weight does not remove a NaN from a sum. + + +@pytest.mark.parametrize("masked", [4, 0, 9]) +def test_masked_nan_sample_is_ignored(minimizer, masked): + # The start estimate of Pol1 uses the first and the last sample. + x = np.arange(10.0) + y = 2.0 + 0.5 * x + y_err = np.ones_like(y) + y[masked] = np.nan + y_err[masked] = 0.0 + result = aare.Pol1(minimizer=minimizer, compute_errors=True).fit(x, y, y_err) + + np.testing.assert_allclose(result["par"], [2.0, 0.5], atol=1e-5) + assert np.all(result["par_err"] > 0.0) + assert result["chi2"][0] == pytest.approx(0.0, abs=1e-6) + + +def test_too_few_unmasked_samples_fail(minimizer): + x = np.arange(10.0) + y = np.full_like(x, np.nan) + y_err = np.zeros_like(x) + y[0], y_err[0] = 2.0, 1.0 # one sample left for two free parameters + result = aare.Pol1(minimizer=minimizer).fit(x, y, y_err) + assert np.all(result["par"] == 0.0) + assert np.all(result["chi2"] == 0.0) + + +def test_masked_samples_in_cube_match_1d_fits(minimizer): + x = np.linspace(-6.0, 6.0, 61) + y = np.stack([gaussian(x, 80.0 + 10.0 * c, 0.8, 1.2) for c in range(4)]) + y = y.reshape(1, 4, x.size) + y_err = np.ones_like(y) + for col, masked in [(1, 10), (2, 0), (3, 34)]: + y[0, col, masked] = np.nan + y_err[0, col, masked] = 0.0 + model = aare.Gaussian(minimizer=minimizer, compute_errors=True) + cube = model.fit(x, y, y_err, n_threads=1) + + for col in range(4): + single = model.fit(x, y[0, col], y_err[0, col]) + np.testing.assert_allclose(cube["par"][0, col], single["par"], atol=1e-12) + np.testing.assert_allclose(cube["par_err"][0, col], single["par_err"], atol=1e-12) + np.testing.assert_allclose( + cube["par"][0, col], [80.0 + 10.0 * col, 0.8, 1.2], rtol=1e-3 + ) + + +def test_parameter_whose_minimum_lies_on_its_limit_reports_no_error(minimizer): + # At a minimum that lies exactly on the limit the gradient is zero. + x = np.array([-1.0, -1.0, 1.0, 1.0]) + model = aare.Pol1(minimizer=minimizer, compute_errors=True) + model.SetParLimits("p0", 0.0, 1.0) + result = model.fit(x, 1.0 + 2.0 * x) + + assert result["par"][0] == pytest.approx(1.0, abs=1e-6) + assert result["par_err"][0] == 0.0 + assert result["par_err"][1] == pytest.approx(0.5, rel=1e-6) + + +@pytest.mark.parametrize("amplitude", [0.1, 1e-4]) +def test_fit_does_not_stop_at_start_without_a_step(minimizer, amplitude): + # The EDM tolerance is absolute in chi2 units: an unweighted fit of small + # data passes it at the start estimate (sigma 1.3617) already. + x = np.linspace(-6.0, 6.0, 61) + y = gaussian(x, amplitude, 0.8, 1.3) + result = aare.Gaussian(minimizer=minimizer).fit(x, y) + assert abs(result["par"][2] - 1.3) < 0.03 # Migrad ends at 1.3205 + + one_call = aare.Gaussian(minimizer=minimizer, max_calls=1).fit(x, y) + assert np.all(one_call["par"] == 0.0) + + +@pytest.mark.parametrize("lo, span", [(8000, 1000), (20000, 500), (60000, 200)]) +def test_polynomial_far_from_zero_keeps_its_precision(minimizer, lo, span): + # Nearly collinear columns 1, x, x^2: VarPro's normal equations lost six + # digits at x = 20000. + x = np.linspace(lo, lo + span, 40) + c = np.array([1.5 - 2e-3 * lo + 3e-7 * lo**2, 2e-3 - 6e-7 * lo, 3e-7]) + y = c[0] + c[1] * x + c[2] * x**2 + result = aare.Pol2(minimizer=minimizer).fit(x, y) + np.testing.assert_allclose(result["par"], c, rtol=1e-9) diff --git a/src/fit/DampedGaussNewton.hpp b/src/fit/DampedGaussNewton.hpp index 00050142..fdf1670b 100644 --- a/src/fit/DampedGaussNewton.hpp +++ b/src/fit/DampedGaussNewton.hpp @@ -154,15 +154,19 @@ bool gauss_newton_errors(const Mat &H, const std::array &skip, * points outward is frozen for that iteration. * - Convergence in Minuit's convention: the iteration stops when the * estimated distance to the minimum of chi2 (EDM, from the Gauss-Newton - * Hessian) drops below 0.002 * tolerance, or when the step becomes - * negligible. + * Hessian) drops below 0.002 * tolerance. A start point that passes + * already still takes one undamped step, see converged(). When the damped + * step becomes negligible for every parameter before that, the iteration + * has stalled and the fit is valid only if the EDM test passes at that + * point. * - max_calls bounds the number of evaluations: the initial point and every * trial point count, so an iteration costs one evaluation, or two when a * limit is crossed. * - With reject_degenerate, a trial point at which the model has lost its * sensitivity to a free parameter (its Hessian diagonal collapsed by more - * than eight orders of magnitude) is rejected like an inadmissible point, - * so the damping grows and a shorter step is tried instead. + * than eight orders of magnitude relative to the current point or to the + * start point) is rejected like an inadmissible point, so the damping + * grows and a shorter step is tried instead. * - With polish, a converged fit takes one more undamped Gauss-Newton step * and keeps it when it lowers chi2, see polish_step(). * @@ -194,6 +198,7 @@ template class DampedGaussNewton { const Vec &upper, const std::array &fixed, const Options &opt) { nfree_ = 0; + accepted_ = 0; for (int k = 0; k < NP; ++k) { lower_[k] = lower[k]; upper_[k] = upper[k]; @@ -214,6 +219,8 @@ template class DampedGaussNewton { if (!evaluate(eval, q_, cur_)) return res; res.calls = 1; + for (int a = 0; a < nfree_; ++a) + start_diag_[a] = cur_.H[free_[a]][free_[a]]; const double edm_target = 0.002 * opt.tolerance; const int max_calls = opt.max_calls; const bool reject_degenerate = opt.reject_degenerate; @@ -224,7 +231,9 @@ template class DampedGaussNewton { extract_free(); mark_active(); if (gmax_ < 1e-12) { - res.valid = true; + edm(); // the undamped step for converged() + res.valid = + converged(eval, res, max_calls, reject_degenerate, false); break; } // Solve for the damped step first. The solve yields @@ -237,10 +246,9 @@ template class DampedGaussNewton { if (e < edm_target) { // A point already deep inside the tolerance needs no // polishing. - if (opt.polish && e >= 0.01 * edm_target && - res.calls < max_calls) - polish_step(eval, res, reject_degenerate); - res.valid = true; + res.valid = + converged(eval, res, max_calls, reject_degenerate, + opt.polish && e >= 0.01 * edm_target); break; } } @@ -253,8 +261,16 @@ template class DampedGaussNewton { break; continue; } - double step2 = 0.0; - double q2 = 0.0; + // A stalled iteration is not a minimum: repeated rejections + // shrink the damped step until it is negligible anywhere. The + // test looks at the step before it is reflected or truncated at + // a limit, so a step that a limit folds back onto the current + // point still gets its truncated alternative evaluated. + if (negligible_step()) { + res.valid = edm() < edm_target; + break; + } + bool moved = false; bool crossed = false; for (int a = 0; a < nfree_; ++a) { const int k = free_[a]; @@ -271,17 +287,16 @@ template class DampedGaussNewton { } q_new_[a] = reflected; q_alt_[a] = truncated; - const double d = q_new_[a] - q_[a]; - step2 += d * d; - q2 += q_[a] * q_[a]; + moved = moved || reflected != q_[a]; } - if (std::sqrt(step2) <= 1e-8 * (std::sqrt(q2) + 1e-8)) { - res.valid = true; - break; + // A reflection that lands on the current point needs no + // evaluation; the truncated step below is the only candidate. + bool ok = false; + if (moved) { + ++res.calls; + ok = evaluate(eval, q_new_, trial_) && + !(reject_degenerate && degenerate(trial_)); } - ++res.calls; - bool ok = evaluate(eval, q_new_, trial_) && - !(reject_degenerate && degenerate(trial_)); if (crossed && res.calls < max_calls) { // Also try the step truncated at the limit, keep the better. ++res.calls; @@ -300,6 +315,7 @@ template class DampedGaussNewton { std::swap(cur_, trial_); for (int a = 0; a < nfree_; ++a) p_[free_[a]] = q_[a]; + ++accepted_; const double t = 2.0 * rho - 1.0; u *= std::max(1.0 / 3.0, 1.0 - t * t * t); nu = 2.0; @@ -369,31 +385,56 @@ template class DampedGaussNewton { } // A parameter sitting on a limit with the descent direction pointing - // outward is frozen for this iteration. The test deliberately looks at - // the last evaluated point p_, which after a rejected step is the trial - // rather than the current point q_: a rejected trial that came off the - // limit unfreezes the parameter, and the reflected step can then leave - // the limit. Testing q_ instead was measured to leave more fits stuck on - // a limit with a higher chi2 (FallingScurve and GaussianChargeSharing - // cubes), although it helped Pol2. + // outward is frozen for this iteration. Position and gradient are both + // taken at the current point q_, so a parameter stays frozen at the end + // only where the constrained minimum lies on its limit. An earlier + // version tested the last evaluated point p_, a rejected trial, so that a + // trial that came off a limit unfroze the parameter: testing q_ had left + // more FallingScurve and GaussianChargeSharing fits stuck with the width + // on its lower limit and a higher chi2. But pairing a trial position with + // the current gradient could freeze an interior parameter and stop the + // fit there. The width no longer reaches its limit that way, because + // LevenbergMarquardt now rejects degenerate trial points as + // VariableProjection does, see degenerate(). void mark_active() { for (int a = 0; a < nfree_; ++a) { const int k = free_[a]; - active_[a] = (p_[k] <= lower_[k] && g_[a] > 0.0) || - (p_[k] >= upper_[k] && g_[a] < 0.0); + active_[a] = (q_[a] <= lower_[k] && g_[a] > 0.0) || + (q_[a] >= upper_[k] && g_[a] < 0.0); } } + // True when the damped step is negligible for every parameter that is + // not frozen: below 1e-8 of its value plus its own scale 1/sqrt(H_aa), + // the parameter change that raises chi2 by about one. Judging each + // parameter on its own scale keeps a large parameter (an amplitude of + // 1e10, a peak position at 1e9) from hiding the steps of the others. + bool negligible_step() const { + for (int a = 0; a < nfree_; ++a) { + if (active_[a]) + continue; + const double scale = + std::abs(q_[a]) + 1.0 / std::sqrt(std::max(jtj_[a][a], 1e-300)); + if (!(std::abs(delta_[a]) <= 1e-8 * scale)) + return false; + } + return true; + } + // A trial point where the Hessian diagonal of a free parameter collapsed - // relative to the current point: the model no longer responds to that - // parameter there, so the iteration could stop at a spurious stationary - // point (a step function with the width driven to zero, for example). + // relative to the current point or to the start point: the model no + // longer responds to that parameter there, so the iteration could stop + // at a spurious stationary point (a step function with the width driven + // to zero, for example). The start point catches a collapse spread over + // many accepted steps, each of which passes the test against its + // predecessor. bool degenerate(const Normal &trial) const { for (int a = 0; a < nfree_; ++a) { if (active_[a]) continue; const int k = free_[a]; - if (!(trial.H[k][k] >= 1e-8 * cur_.H[k][k])) + const double ref = std::max(cur_.H[k][k], start_diag_[a]); + if (!(trial.H[k][k] >= 1e-8 * ref)) return true; } return false; @@ -425,6 +466,25 @@ template class DampedGaussNewton { return e; } + // The convergence test passed. A point that passes it before any step + // was accepted may merely lie within the tolerance of chi2's scale, as + // the start estimate of an unweighted fit of data of small magnitude + // does, so it still takes the undamped step of polish_step(), as Migrad + // always iterates at least once. Without the budget for that step the + // fit has not converged. With polish, a converged fit takes the same + // step when the budget allows. Returns whether the fit is valid. + template + bool converged(Eval &eval, MinimizerResult &res, int max_calls, + bool reject_degenerate, bool polish) { + const bool first = accepted_ == 0 && nfree_ > 0; + if (!first && !polish) + return true; + if (res.calls >= max_calls) + return !first; + polish_step(eval, res, reject_degenerate); + return true; + } + // One undamped Gauss-Newton step from a point that passed the EDM test, // kept when it lowers F. The test stops the iteration within the // tolerance, but a solver that converges in few iterations stops while @@ -507,6 +567,7 @@ template class DampedGaussNewton { } int nfree_ = 0; + int accepted_ = 0; // trial points accepted in this fit std::array free_{}; std::array free_pos_{}; std::array fixed_{}; @@ -520,6 +581,7 @@ template class DampedGaussNewton { Vec g_{}; Vec delta_{}; Vec gn_step_{}; + Vec start_diag_{}; // Hessian diagonal at the start point Mat jtj_{}; // The two solves that bound the EDM differ by rounding only when the // damping is negligible; the margin keeps the stop decision exact. diff --git a/src/fit/Fit.cpp b/src/fit/Fit.cpp index 37187483..8b34c177 100644 --- a/src/fit/Fit.cpp +++ b/src/fit/Fit.cpp @@ -7,7 +7,9 @@ #include "aare/Models.hpp" #include "aare/utils/par.hpp" #include "aare/utils/task.hpp" +#include #include +#include #include #include #include @@ -141,6 +143,34 @@ template class PixelFitter { bool fit(const FitModel &model, NDView x, NDView y, NDView y_err, double *par_out, double *err_out, double &chi2) { + // Samples with a zero error are ignored. They are dropped here, + // before the start estimate and the minimizer: a masked sample may + // hold NaN, and a zero weight does not remove a NaN from a sum. + if (y_err.size() > 0 && + std::find(y_err.begin(), y_err.end(), 0.0) != y_err.end()) { + keep_x_.clear(); + keep_y_.clear(); + keep_err_.clear(); + for (ssize_t i = 0; i < y.size(); ++i) { + if (y_err[i] == 0.0) + continue; + keep_x_.push_back(x[i]); + keep_y_.push_back(y[i]); + keep_err_.push_back(y_err[i]); + } + const auto n = static_cast(keep_x_.size()); + x = NDView(keep_x_.data(), {n}); + y = NDView(keep_y_.data(), {n}); + y_err = NDView(keep_err_.data(), {n}); + } + // Fewer samples than free parameters leave the fit undetermined, + // and the start estimates read at least two samples. + ssize_t n_free = 0; + for (std::size_t k = 0; k < Model::npar; ++k) + n_free += model.is_user_fixed(static_cast(k)) ? 0 : 1; + if (y.size() < std::max(n_free, 2)) + return fail(par_out, err_out, chi2); + const auto start = start_values(model, x, y); if (model.minimizer() == Minimizer::VarPro) { @@ -150,7 +180,7 @@ template class PixelFitter { if constexpr (model::is_separable::value) { const auto res = vp_.fit(model, x, y, y_err, start, par_out, err_out); - if (res.valid) { + if (res.valid && finite(par_out, err_out, res.chi2)) { chi2 = res.chi2; return true; } @@ -164,15 +194,40 @@ template class PixelFitter { } private: + // The built-in solvers do not guard every operation against NaN or + // infinity, so a result is checked before it is reported as a valid fit. + static bool finite(const double *par, const double *err, double chi2) { + bool ok = std::isfinite(chi2); + for (std::size_t k = 0; k < Model::npar; ++k) + ok = ok && std::isfinite(par[k]) && (!err || std::isfinite(err[k])); + return ok; + } + bool fit_lm(const FitModel &model, NDView x, NDView y, NDView y_err, const std::array &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; + if (res.valid && finite(par_out, err_out, res.chi2)) { + chi2 = res.chi2; + return true; + } + return fail(par_out, err_out, chi2); } + // A failed fit reports zeros. + static bool fail(double *par_out, double *err_out, double &chi2) { + std::fill(par_out, par_out + Model::npar, 0.0); + if (err_out) + std::fill(err_out, err_out + Model::npar, 0.0); + chi2 = 0.0; + return false; + } + + // The samples of the current pixel that have a nonzero error, see fit(). + std::vector keep_x_; + std::vector keep_y_; + std::vector keep_err_; LevenbergMarquardt lm_; std::conditional_t::value, VariableProjection, NotSeparable> @@ -185,6 +240,7 @@ template NDArray fit_pixel(const FitModel &model, NDView x, NDView y, NDView y_err) { constexpr std::size_t npar = Model::npar; + detail::check_fit_pixel_shapes(x, y, y_err); const bool want_errors = model.compute_errors(); const auto result_size = static_cast(want_errors ? (2 * npar + 1) : (npar + 1)); diff --git a/src/fit/Fit.test.cpp b/src/fit/Fit.test.cpp index 3fa39eaa..992500e8 100644 --- a/src/fit/Fit.test.cpp +++ b/src/fit/Fit.test.cpp @@ -9,6 +9,7 @@ #include #include #include +#include #include #include #include @@ -439,6 +440,21 @@ TEST_CASE("FitModel validates parameter names, indices and limits", "[fit]") { CHECK(model.minimizer() == Minimizer::LevenbergMarquardt); } +TEST_CASE("fit_pixel rejects arrays of different sizes", "[fit]") { + const auto minimizer = + GENERATE(Minimizer::Migrad, Minimizer::Fumili, + Minimizer::LevenbergMarquardt, Minimizer::VarPro); + INFO("minimizer " << name(minimizer)); + auto d = gaussian_data(); + const auto model = make_model(minimizer); + NDView short_y(d.y.data(), {n_points - 1}); + NDView short_err(d.y_err.data(), {n_points - 1}); + CHECK_THROWS_AS(aare::fit_pixel(model, d.x.view(), short_y), + std::runtime_error); + CHECK_THROWS_AS(aare::fit_pixel(model, d.x.view(), d.y.view(), short_err), + std::runtime_error); +} + TEST_CASE("Copying a FitModel keeps the minimizer", "[fit]") { const aare::FitModel original( 0, 100, 0.5, false, Minimizer::LevenbergMarquardt); @@ -925,3 +941,368 @@ TEST_CASE("A parameter ending exactly on a limit reports an error of 0", CHECK(res(2) == 0.0); CHECK(res(3) == 0.0); } + +// --------------------------------------------------------------------------- +// Regressions of the built-in solvers: fits that used to be reported as valid +// although they had stopped away from the minimum or returned garbage. +// --------------------------------------------------------------------------- + +namespace { + +// Standard normal noise from std::mt19937, whose output the standard fixes, +// and a Box-Muller transform written out here: std::normal_distribution +// differs between standard libraries. +std::vector normal_noise(std::size_t n, unsigned seed) { + std::mt19937 rng(seed); + const auto uniform = [&rng] { + return (static_cast(rng()) + 0.5) / 4294967296.0; + }; + std::vector out(n); + for (std::size_t i = 0; i < n; i += 2) { + const double r = std::sqrt(-2.0 * std::log(uniform())); + const double phi = 2.0 * std::acos(-1.0) * uniform(); + out[i] = r * std::cos(phi); + if (i + 1 < n) + out[i + 1] = r * std::sin(phi); + } + return out; +} + +} // namespace + +TEST_CASE("A Gaussian fit to pure noise fails or stays finite and bounded", + "[fit]") { + // Fitting noise, a narrow width can collapse between two scan points; + // the amplitude that fits a single point then diverges. Such a fit must + // be reported as failed, not as a valid result with |A| up to 1e159. + const auto minimizer = + GENERATE(Minimizer::Migrad, Minimizer::Fumili, + Minimizer::LevenbergMarquardt, Minimizer::VarPro); + INFO("minimizer " << name(minimizer)); + constexpr ssize_t rows = 20; + constexpr ssize_t cols = 20; + NDArray x({n_points}); + for (ssize_t i = 0; i < n_points; ++i) + x(i) = -6.0 + 0.2 * static_cast(i); + NDArray y({rows, cols, n_points}); + const auto noise = normal_noise(static_cast(y.size()), 1); + std::copy(noise.begin(), noise.end(), y.begin()); + + NDArray par({rows, cols, 3}); + NDArray chi2({rows, cols}); + const auto model = make_model(minimizer); + aare::fit_3d(model, x.view(), y.view(), NDView{}, par.view(), + NDView{}, chi2.view(), 1); + + int absurd = 0; + for (ssize_t row = 0; row < rows; ++row) { + for (ssize_t col = 0; col < cols; ++col) { + bool ok = std::isfinite(chi2(row, col)); + for (ssize_t k = 0; k < 3; ++k) + ok = ok && std::isfinite(par(row, col, k)) && + std::abs(par(row, col, k)) < 1e6; + absurd += ok ? 0 : 1; + } + } + CHECK(absurd == 0); +} + +TEST_CASE("A step reflected back onto the start point is not convergence", + "[fit]") { + // p0 starts at 0 with limits [-1, 1]; the damped step of +2 reflects at + // the upper limit back onto 0. The step truncated at the limit must still + // be tried: the constrained minimum is p0 = 1. + const auto minimizer = + GENERATE(Minimizer::Migrad, Minimizer::Fumili, + Minimizer::LevenbergMarquardt, Minimizer::VarPro); + INFO("minimizer " << name(minimizer)); + NDArray x({4}); + NDArray y({4}, 2.2); + for (ssize_t i = 0; i < 4; ++i) + x(i) = static_cast(i); + + auto model = make_model(minimizer); + model.SetParameter("p0", 0.0); + model.FixParameter("p1", 0.0); + model.SetParLimits("p0", -1.0, 1.0); + const auto res = aare::fit_pixel(model, x.view(), y.view()); + CHECK(res(0) == Approx(1.0).margin(limit_margin(minimizer))); + CHECK(res(2) == Approx(5.76).epsilon(1e-6)); +} + +TEST_CASE("A large parameter does not hide the steps of the others", "[fit]") { + // The negligible-step test used one relative norm over all parameters, + // so an amplitude of 1e10 or a position of 1e9 stopped the fit at the + // start estimate of the width. + const auto minimizer = + GENERATE(Minimizer::Migrad, Minimizer::Fumili, + Minimizer::LevenbergMarquardt, Minimizer::VarPro); + INFO("minimizer " << name(minimizer)); + const auto [amplitude, offset] = GENERATE( + std::pair{1e8, 0.0}, std::pair{1e10, 0.0}, std::pair{120.0, 1e9}); + INFO("amplitude " << amplitude << ", offset " << offset); + + Data d{NDArray({n_points}), NDArray({n_points}), + NDArray({n_points})}; + fill_gaussian(d.x, d.y, amplitude, 0.8, 1.3); + for (ssize_t i = 0; i < n_points; ++i) { + d.x(i) += offset; + d.y_err(i) = std::sqrt(d.y(i)) + 1.0; + } + const auto model = make_model(minimizer, false, 500); + const auto res = + aare::fit_pixel(model, d.x.view(), d.y.view(), d.y_err.view()); + // The EDM tolerance leaves a few 1e-4 relative, far below the errors; + // the stuck fits were 0.7 % (A = 1e8) and 5 % (start estimate) off. + CHECK(res(0) == Approx(amplitude).epsilon(1e-3)); + CHECK(res(1) - offset == Approx(0.8).margin(1e-3)); + CHECK(res(2) == Approx(1.3).epsilon(1e-3)); +} + +TEST_CASE("A rejected trial on a limit does not freeze an interior parameter", + "[fit]") { + // mu starts at 3 inside [0, 5]. A trial truncated onto the lower limit is + // rejected; testing that trial's position against the gradient of the + // accepted point froze mu, and the fit stopped at the start value. + const auto minimizer = + GENERATE(Minimizer::Migrad, Minimizer::Fumili, + Minimizer::LevenbergMarquardt, Minimizer::VarPro); + INFO("minimizer " << name(minimizer)); + constexpr ssize_t n = 11; + const double mu = std::sqrt(2.0 * std::log(2.0)); + NDArray x({n}); + NDArray y({n}); + for (ssize_t i = 0; i < n; ++i) { + x(i) = -0.1 + 0.02 * static_cast(i); + y(i) = 10.0 * std::exp(-0.5 * (x(i) - mu) * (x(i) - mu)); + } + auto model = make_model(minimizer); + model.FixParameter("A", 10.0); + model.FixParameter("sigma", 1.0); + model.SetParameter("mu", 3.0); + model.SetParLimits("mu", 0.0, 5.0); + const auto res = aare::fit_pixel(model, x.view(), y.view()); + CHECK(res(1) == Approx(mu).margin(1e-3)); + CHECK(res(3) < 1e-4); +} + +// --------------------------------------------------------------------------- +// Samples with a zero error are ignored. They may hold NaN (a masked dead +// frame), which a zero weight does not remove from a sum, so they must not +// reach the start estimate or the minimizer at all. +// --------------------------------------------------------------------------- + +TEST_CASE("Samples with a zero error are ignored even when NaN", "[fit]") { + const auto minimizer = + GENERATE(Minimizer::Migrad, Minimizer::Fumili, + Minimizer::LevenbergMarquardt, Minimizer::VarPro); + INFO("minimizer " << name(minimizer)); + const double nan = std::numeric_limits::quiet_NaN(); + + constexpr ssize_t n = 10; + NDArray x({n}); + NDArray y({n}); + NDArray y_err({n}, 1.0); + for (ssize_t i = 0; i < n; ++i) { + x(i) = static_cast(i); + y(i) = 2.0 + 0.5 * x(i); + } + const auto line = make_model(minimizer, true); + + // The start estimate of Pol1 uses the first and the last sample. + const auto masked = GENERATE(as{}, 4, 0, n - 1); + SECTION("masked NaN on a line") { + INFO("masked sample " << masked); + y(masked) = nan; + y_err(masked) = 0.0; + const auto res = + aare::fit_pixel(line, x.view(), y.view(), y_err.view()); + CHECK(res(0) == Approx(2.0).epsilon(1e-5).margin(1e-6)); + CHECK(res(1) == Approx(0.5).epsilon(1e-5).margin(1e-6)); + CHECK(res(2) > 0.0); + CHECK(res(4) == Approx(0.0).margin(1e-6)); + } + + SECTION("masked NaN at the peak of a Gaussian") { + auto d = gaussian_data(); + constexpr ssize_t peak = 34; // x = 0.8 + d.y(peak) = nan; + d.y_err(peak) = 0.0; + const auto model = make_model(minimizer); + const auto res = + aare::fit_pixel(model, d.x.view(), d.y.view(), d.y_err.view()); + const double tol = truth_tolerance(minimizer); + CHECK(res(0) == Approx(120.0).epsilon(tol)); + CHECK(res(1) == Approx(0.8).epsilon(tol)); + CHECK(res(2) == Approx(1.3).epsilon(tol)); + } + + SECTION("too few samples left fail cleanly") { + // One sample left for two free parameters, then none at all. + for (ssize_t i = 1; i < n; ++i) { + y(i) = nan; + y_err(i) = 0.0; + } + auto res = aare::fit_pixel(line, x.view(), y.view(), y_err.view()); + for (ssize_t k = 0; k < res.size(); ++k) + CHECK(res(k) == 0.0); + y_err(0) = 0.0; + res = aare::fit_pixel(line, x.view(), y.view(), y_err.view()); + for (ssize_t k = 0; k < res.size(); ++k) + CHECK(res(k) == 0.0); + } +} + +TEST_CASE("Masked samples in a data cube match per-pixel fits", "[fit]") { + // The kept samples of each pixel go through reused buffers, which must + // not carry anything over from one pixel to the next. + const auto minimizer = + GENERATE(Minimizer::Migrad, Minimizer::Fumili, + Minimizer::LevenbergMarquardt, Minimizer::VarPro); + INFO("minimizer " << name(minimizer)); + constexpr ssize_t rows = 2; + constexpr ssize_t cols = 3; + // masked sample per pixel, -1 for none + constexpr ssize_t masked[rows][cols] = {{-1, 10, 0}, + {34, -1, n_points - 1}}; + NDArray x({n_points}); + NDArray tmp({n_points}); + NDArray y({rows, cols, n_points}); + NDArray y_err({rows, cols, n_points}, 1.0); + for (ssize_t row = 0; row < rows; ++row) { + for (ssize_t col = 0; col < cols; ++col) { + fill_gaussian(x, tmp, 80.0 + 10.0 * static_cast(col), 0.8, + 1.0 + 0.2 * static_cast(row)); + for (ssize_t i = 0; i < n_points; ++i) + y(row, col, i) = tmp(i); + if (masked[row][col] >= 0) { + y(row, col, masked[row][col]) = + std::numeric_limits::quiet_NaN(); + y_err(row, col, masked[row][col]) = 0.0; + } + } + } + + NDArray par({rows, cols, 3}); + NDArray err({rows, cols, 3}); + NDArray chi2({rows, cols}); + const auto model = make_model(minimizer, true); + aare::fit_3d(model, x.view(), y.view(), y_err.view(), par.view(), + err.view(), chi2.view(), 1); + + for (ssize_t row = 0; row < rows; ++row) { + for (ssize_t col = 0; col < cols; ++col) { + INFO("pixel " << row << ", " << col); + NDView yv(&y(row, col, 0), {n_points}); + NDView ev(&y_err(row, col, 0), {n_points}); + const auto single = aare::fit_pixel(model, x.view(), yv, ev); + for (ssize_t k = 0; k < 3; ++k) { + CHECK(par(row, col, k) == Approx(single(k)).margin(1e-12)); + CHECK(err(row, col, k) == Approx(single(3 + k)).margin(1e-12)); + } + CHECK(chi2(row, col) == Approx(single(6)).margin(1e-12)); + // Close to the truth within the EDM tolerance; the exact check + // is the agreement with the per-pixel fit above. + const double tol = 1e-3; + CHECK(par(row, col, 0) == + Approx(80.0 + 10.0 * static_cast(col)).epsilon(tol)); + CHECK(par(row, col, 1) == Approx(0.8).epsilon(tol)); + CHECK(par(row, col, 2) == + Approx(1.0 + 0.2 * static_cast(row)).epsilon(tol)); + } + } +} + +TEST_CASE("A parameter whose minimum lies on its limit reports no error", + "[fit]") { + // At a minimum that lies exactly on a limit the gradient is zero, so a + // test for an outward gradient does not see the parameter as bound. + const auto minimizer = + GENERATE(Minimizer::Migrad, Minimizer::Fumili, + Minimizer::LevenbergMarquardt, Minimizer::VarPro); + INFO("minimizer " << name(minimizer)); + + SECTION("linear parameter") { + NDArray x({4}); + NDArray y({4}); + const double xs[4] = {-1.0, -1.0, 1.0, 1.0}; + for (ssize_t i = 0; i < 4; ++i) { + x(i) = xs[i]; + y(i) = 1.0 + 2.0 * xs[i]; + } + auto model = make_model(minimizer, true); + model.SetParLimits("p0", 0.0, 1.0); + const auto res = aare::fit_pixel(model, x.view(), y.view()); + CHECK(res(0) == Approx(1.0).margin(limit_margin(minimizer))); + CHECK(res(2) == 0.0); + CHECK(res(3) == Approx(0.5).epsilon(1e-6)); // 1 / sqrt(sum x^2) + } + + SECTION("nonlinear parameter") { + auto d = gaussian_data(); + auto model = make_model(minimizer, true); + model.SetParLimits("mu", -5.0, 0.8); + const auto res = + aare::fit_pixel(model, d.x.view(), d.y.view(), d.y_err.view()); + CHECK(res(1) == Approx(0.8).margin(1e-6)); + CHECK(res(4) == 0.0); + CHECK(res(3) > 0.0); + CHECK(res(5) > 0.0); + } +} + +TEST_CASE("A fit does not stop at its start point without taking a step", + "[fit]") { + // The EDM tolerance is absolute in chi2 units. An unweighted fit of data + // of small magnitude passes it at the start estimate already; the fit + // must still take a step, as Migrad always iterates at least once. + const auto minimizer = + GENERATE(Minimizer::Migrad, Minimizer::Fumili, + Minimizer::LevenbergMarquardt, Minimizer::VarPro); + INFO("minimizer " << name(minimizer)); + const double amplitude = GENERATE(0.1, 1e-4); + INFO("amplitude " << amplitude); + Data d{NDArray({n_points}), NDArray({n_points}), + NDArray({n_points}, 1.0)}; + fill_gaussian(d.x, d.y, amplitude, 0.8, 1.3); + + // The start estimate of sigma is 1.3617; Migrad ends at 1.3205. + const auto model = make_model(minimizer); + const auto res = aare::fit_pixel(model, d.x.view(), d.y.view()); + CHECK(std::abs(res(2) - 1.3) < 0.03); + + // Without the budget for a single step the fit has not converged. + const auto one_call = + make_model(minimizer, false, 1); + const auto none = aare::fit_pixel(one_call, d.x.view(), d.y.view()); + for (ssize_t k = 0; k < none.size(); ++k) + CHECK(none(k) == 0.0); +} + +TEST_CASE("A polynomial far from x = 0 keeps its precision", "[fit]") { + // On a narrow range far from 0 the columns 1, x, x^2 are nearly + // collinear. VarPro solves them from the normal equations and lost six + // digits at x = 20000; from about x = 60000 it falls back to LM. + const auto minimizer = + GENERATE(Minimizer::Migrad, Minimizer::Fumili, + Minimizer::LevenbergMarquardt, Minimizer::VarPro); + INFO("minimizer " << name(minimizer)); + const auto [lo, span] = + GENERATE(std::pair{8000.0, 1000.0}, std::pair{20000.0, 500.0}, + std::pair{60000.0, 200.0}); + INFO("x in [" << lo << ", " << lo + span << "]"); + // 1.5 + 2e-3 (x - lo) + 3e-7 (x - lo)^2 expanded in powers of x. + const double c[3] = {1.5 - 2e-3 * lo + 3e-7 * lo * lo, 2e-3 - 6e-7 * lo, + 3e-7}; + constexpr ssize_t n = 40; + NDArray x({n}); + NDArray y({n}); + for (ssize_t i = 0; i < n; ++i) { + x(i) = lo + span * static_cast(i) / static_cast(n - 1); + y(i) = c[0] + c[1] * x(i) + c[2] * x(i) * x(i); + } + const auto model = make_model(minimizer); + const auto res = aare::fit_pixel(model, x.view(), y.view()); + for (ssize_t k = 0; k < 3; ++k) + CHECK(res(k) == Approx(c[k]).epsilon(1e-9)); +} diff --git a/src/fit/FitHelpers.hpp b/src/fit/FitHelpers.hpp index 062fa640..a82599cf 100644 --- a/src/fit/FitHelpers.hpp +++ b/src/fit/FitHelpers.hpp @@ -35,6 +35,16 @@ std::array start_values(const FitModel &model, return start; } +/** @brief Validate the array sizes passed to fit_pixel. */ +inline void check_fit_pixel_shapes(NDView x, NDView y, + NDView y_err) { + if (x.size() != y.size()) + throw std::runtime_error("fit_pixel: x.size() must match y.size()."); + if (y_err.size() > 0 && y_err.size() != y.size()) + throw std::runtime_error( + "fit_pixel: y and y_err must have identical size."); +} + /** @brief Validate the array shapes passed to fit_3d. */ template void check_fit_3d_shapes(NDView x, NDView y, diff --git a/src/fit/LevenbergMarquardt.hpp b/src/fit/LevenbergMarquardt.hpp index c862707e..8a6a7213 100644 --- a/src/fit/LevenbergMarquardt.hpp +++ b/src/fit/LevenbergMarquardt.hpp @@ -26,8 +26,9 @@ namespace aare::detail { * equations g = J^T r and J^T J, which DampedGaussNewton iterates on: it * handles fixed parameters, limits, damping and the convergence test, see * there. Rejected trials are rare with its damping rule, and a rejected one - * only wastes a single pass. As in Minuit2, max_calls = 0 selects the - * default budget 200 + 100 * npar + 5 * npar^2. + * only wastes a single pass. Trial points at which the model loses its + * sensitivity to a free parameter are rejected. As in Minuit2, max_calls = 0 + * selects the default budget 200 + 100 * npar + 5 * npar^2. * * Errors are Gauss-Newton estimates, sqrt(diag((J^T J)^-1)), over the free * parameters that do not sit on a limit. Fixed and limit-bound parameters @@ -95,6 +96,11 @@ template class LevenbergMarquardt { // Start cautiously: the estimated start values of the linear // parameters can be far off, and the first steps are then wild. opt.damping0 = 0.1; + // A width driven onto its lower limit turns a model into a step + // function whose chi2 no longer responds to the width and jumps + // with the position; the iteration then stalls there. Such trial + // points are rejected, as in VariableProjection. + opt.reject_degenerate = true; const Result res = driver_.fit(eval, start, lower, upper, fixed, opt); if (!res.valid) { std::fill(par_out, par_out + npar, 0.0); diff --git a/src/fit/VariableProjection.hpp b/src/fit/VariableProjection.hpp index f0059ca4..0ea0cdbc 100644 --- a/src/fit/VariableProjection.hpp +++ b/src/fit/VariableProjection.hpp @@ -101,6 +101,7 @@ template class VariableProjection { } pvec_.assign(start.begin(), start.end()); have_last_ = false; + have_ref_ = false; nfree_lin_ = 0; for (int j = 0; j < nlin; ++j) { const auto k = static_cast(Traits::linear_index(j)); @@ -195,8 +196,9 @@ template class VariableProjection { private: // F = chi2 / 2, gradient and Kaufman Hessian of the projected residual // at the nonlinear point alpha, and the linear solution in out.aux. - // Returns false when the model rejects the point or the basis functions - // are linearly dependent there. + // Returns false when the model rejects the point, or when the basis + // functions are numerically linearly dependent there or one of them has + // vanished on the scan points, see reduce(). bool evaluate(const Vec &alpha, Normal &out) { for (int k = 0; k < nnl; ++k) pvec_[Traits::nonlinear_par[k]] = alpha[k]; @@ -293,13 +295,64 @@ template class VariableProjection { for (int b = 0; b < nfree_lin_; ++b) A[a][b] = M[ja][free_lin_[b]]; } + // A free basis function that has all but vanished on the scan points + // (a Gaussian narrower than the point spacing, centred between two + // points) leaves its coefficient undetermined, and the solve returns + // an absurd one that fits a single point. Its squared norm is compared + // with that at the pixel's first evaluation, the start point. + if (!have_ref_) { + for (int j = 0; j < nlin; ++j) + ref_norm2_[j] = M[j][j]; + have_ref_ = true; + } + for (int a = 0; a < nfree_lin_; ++a) { + const int j = free_lin_[a]; + if (!(M[j][j] >= vanished_basis * ref_norm2_[j])) + return false; + } Mat L{}; if (!cholesky_factor(nfree_lin_, A, L)) return false; + // Nearly collinear basis functions: a pivot of the factorisation far + // below its diagonal element means that the linear solution loses + // about log10(diagonal / pivot) digits, beyond which it is noise. + bool refine = false; + for (int a = 0; a < nfree_lin_; ++a) { + if (!(L[a][a] * L[a][a] >= collinear_basis * A[a][a])) + return false; + refine = refine || L[a][a] * L[a][a] < refine_basis * A[a][a]; + } Vec beta_free{}; cholesky_substitute(nfree_lin_, L, v, beta_free); for (int a = 0; a < nfree_lin_; ++a) pvec_[Traits::linear_index(free_lin_[a])] = beta_free[a]; + // Short of that, the digits lost to the normal equations are + // recovered by one step of iterative refinement on the data (the + // corrected semi-normal equations): the correction solves the same + // equations for Phi^T W r of the residual of the first solution. A + // polynomial fitted far from x = 0 needs it; x in [20000, 20500] + // otherwise costs a Pol2 six digits. + if (refine) { + Vec beta; + for (int j = 0; j < nlin; ++j) + beta[j] = pvec_[Traits::linear_index(j)]; + Vec s{}; + for (ssize_t i = 0; i < n; ++i) { + double r = yd[i]; + for (int j = 0; j < nlin; ++j) + r -= beta[j] * pc[j][i]; + const double rw = Weighted ? r * w2[i] : r; + for (int j = 0; j < nlin; ++j) + s[j] += pc[j][i] * rw; + } + Vec sf{}; + for (int a = 0; a < nlin; ++a) + sf[a] = a < nfree_lin_ ? s[free_lin_[a]] : 0.0; + Vec delta{}; + cholesky_substitute(nfree_lin_, L, sf, delta); + for (int a = 0; a < nfree_lin_; ++a) + pvec_[Traits::linear_index(free_lin_[a])] += delta[a]; + } Vec beta; for (int j = 0; j < nlin; ++j) { beta[j] = pvec_[Traits::linear_index(j)]; @@ -438,6 +491,14 @@ template class VariableProjection { Mat last_M_{}; std::array, ndrv> last_Q_{}; Mat last_N_{}; + // Squared norms of the basis functions at the start point, see reduce(). + bool have_ref_ = false; + Vec ref_norm2_{}; + static constexpr double vanished_basis = 1e-8; + static constexpr double collinear_basis = 1e-10; + // Pivot ratio below which the linear solution is refined, about four + // lost digits. + static constexpr double refine_basis = 1e-4; int nfree_lin_ = 0; std::array free_lin_{}; std::array lin_fixed_{};