diff --git a/src/fit/DampedGaussNewton.hpp b/src/fit/DampedGaussNewton.hpp index 0d45a7b6..00050142 100644 --- a/src/fit/DampedGaussNewton.hpp +++ b/src/fit/DampedGaussNewton.hpp @@ -282,7 +282,7 @@ template class DampedGaussNewton { ++res.calls; bool ok = evaluate(eval, q_new_, trial_) && !(reject_degenerate && degenerate(trial_)); - if (crossed) { + if (crossed && res.calls < max_calls) { // Also try the step truncated at the limit, keep the better. ++res.calls; const bool ok_alt = evaluate(eval, q_alt_, alt_) && @@ -326,9 +326,15 @@ template class DampedGaussNewton { /** @brief F, gradient and Hessian at the final point. */ const Normal &normal() const { return cur_; } - /** @brief True for a fixed parameter or one frozen on a limit. */ + /** + * @brief True for a parameter without an error: fixed, frozen on a limit, + * or ending exactly on one. The frozen flag alone misses the last case, + * because it needs the gradient to point outward and that gradient is + * zero when the unconstrained minimum coincides with the limit. + */ bool skipped(int k) const { - return fixed_[k] || (free_pos_[k] >= 0 && active_[free_pos_[k]]); + return fixed_[k] || (free_pos_[k] >= 0 && active_[free_pos_[k]]) || + p_[k] <= lower_[k] || p_[k] >= upper_[k]; } /** @brief Gauss-Newton errors at the final point, see diff --git a/src/fit/Fit.test.cpp b/src/fit/Fit.test.cpp index 75b9cf8d..3fa39eaa 100644 --- a/src/fit/Fit.test.cpp +++ b/src/fit/Fit.test.cpp @@ -873,3 +873,55 @@ TEST_CASE("fast_exp matches std::exp", "[fit]") { CHECK(aare::model::fast_exp(0.0) == 1.0); CHECK(aare::model::fast_exp(-1000.0) == aare::model::fast_exp(-700.0)); } + +TEST_CASE("Fit with a one-sided limit far from the start value", "[fit]") { + const auto minimizer = + GENERATE(Minimizer::Migrad, Minimizer::Fumili, + Minimizer::LevenbergMarquardt, Minimizer::VarPro); + INFO("minimizer " << name(minimizer)); + constexpr ssize_t n = 20; + NDArray x({n}); + NDArray y({n}); + for (ssize_t i = 0; i < n; ++i) { + x(i) = static_cast(i); + y(i) = 2.0 + 0.5 * x(i); + } + const double inf = std::numeric_limits::infinity(); + const bool upper = GENERATE(true, false); + auto model = make_model(minimizer); + if (upper) + model.SetParLimits(0, -inf, 1e9); + else + model.SetParLimits(0, -1e9, inf); + + const auto res = aare::fit_pixel(model, x.view(), y.view()); + CHECK(res(0) == Approx(2.0).epsilon(1e-6)); + CHECK(res(1) == Approx(0.5).epsilon(1e-6)); +} + +TEST_CASE("A parameter ending exactly on a limit reports an error of 0", + "[fit]") { + const auto minimizer = + GENERATE(Minimizer::Migrad, Minimizer::Fumili, + Minimizer::LevenbergMarquardt, Minimizer::VarPro); + INFO("minimizer " << name(minimizer)); + constexpr ssize_t n = 3; + NDArray x({n}); + NDArray y({n}, 2.0); + for (ssize_t i = 0; i < n; ++i) + x(i) = static_cast(i); + // The unconstrained minimum p0 = 2 coincides with the limit. + const bool lower = GENERATE(true, false); + auto model = make_model(minimizer, true); + model.FixParameter(1, 0.0); + if (lower) + model.SetParLimits(0, 2.0, 3.0); + else + model.SetParLimits(0, 1.0, 2.0); + + const auto res = aare::fit_pixel(model, x.view(), y.view()); + REQUIRE(res.size() == 5); + CHECK(res(0) == Approx(2.0).epsilon(1e-6)); + CHECK(res(2) == 0.0); + CHECK(res(3) == 0.0); +} diff --git a/src/fit/FitMinuit2.cpp b/src/fit/FitMinuit2.cpp index 4c07040a..326c4b11 100644 --- a/src/fit/FitMinuit2.cpp +++ b/src/fit/FitMinuit2.cpp @@ -48,11 +48,13 @@ user_parameters(const FitModel &model, NDView x, } else if (std::isfinite(lo) || std::isfinite(hi)) { // Minuit2's transformation for one-sided limits converges markedly // slower than the two-sided one, so an open side gets a wide bound. + // It is placed relative to the start value, which then always + // lies inside the range however far it is from the finite limit. const double span = std::max(1e6, 1e3 * std::abs(start[i])); if (std::isfinite(lo)) - upar.SetLimits(idx, lo, lo + span); + upar.SetLimits(idx, lo, std::max(start[i], lo) + span); else - upar.SetLimits(idx, hi - span, hi); + upar.SetLimits(idx, std::min(start[i], hi) - span, hi); } if (model.is_user_fixed(idx)) upar.Fix(idx); diff --git a/src/fit/VariableProjection.hpp b/src/fit/VariableProjection.hpp index f94dd9d5..f0059ca4 100644 --- a/src/fit/VariableProjection.hpp +++ b/src/fit/VariableProjection.hpp @@ -169,14 +169,22 @@ template class VariableProjection { driver_.normal().aux[static_cast(j)]; } - // The linear solve knows nothing about limits. + // The linear solve knows nothing about limits. A solution beyond a + // limit fails; one within rounding of a limit is moved onto it, so + // that it counts as on the limit whichever way the rounding went. for (int j = 0; j < nlin; ++j) { if (lin_fixed_[j]) continue; const auto k = static_cast(Traits::linear_index(j)); - if (pvec_[k] < model.lower_limit(k) || - pvec_[k] > model.upper_limit(k)) + const double lo = model.lower_limit(k); + const double hi = model.upper_limit(k); + const double tol = 1e-12 * std::max(1.0, std::abs(pvec_[k])); + if (pvec_[k] < lo - tol || pvec_[k] > hi + tol) return fail(); + if (std::abs(pvec_[k] - lo) <= tol) + pvec_[k] = lo; + else if (std::abs(pvec_[k] - hi) <= tol) + pvec_[k] = hi; } std::copy(pvec_.begin(), pvec_.end(), par_out); if (err_out) @@ -403,9 +411,16 @@ template class VariableProjection { H[a][c] += g[a] * g[c] * wi; } } + // Fixed parameters and parameters that end on a limit have no error. + // The driver knows the nonlinear ones; the linear solve can land a + // linear parameter exactly on a limit without violating it. std::array skip; - for (int k = 0; k < npar; ++k) - skip[k] = model.is_user_fixed(static_cast(k)); + for (int k = 0; k < npar; ++k) { + const auto idx = static_cast(k); + const double v = pvec_[static_cast(k)]; + skip[k] = model.is_user_fixed(idx) || v <= model.lower_limit(idx) || + v >= model.upper_limit(idx); + } if constexpr (nnl > 0) { for (int k = 0; k < nnl; ++k) skip[Traits::nonlinear_par[k]] =