DM-38616: Add GSL interpolation for Sersic gaussian mixture

include/gsl.h

@@ -0,0 +1,26 @@
+#ifndef GAUSS2DFIT_GSL_H
+#define GAUSS2DFIT_GSL_H
+
+#ifdef GAUSS2D_FIT_HAS_GSL
+
+#include <unordered_map>
+
+#include <gsl/gsl_interp.h>
+
+#include "interpolation.h"
+
+namespace gauss2d::fit {
+
+/**
+ * See GSL docs for 1D interpolation types, or (as of 2.7):
+ * https://www.gnu.org/software/gsl/doc/html/interp.html#d-interpolation-types
+ */
+static const std::unordered_map<InterpType, const gsl_interp_type*> GSLInterpTypes{
+        {InterpType::linear, gsl_interp_linear},
+        {InterpType::polynomial, gsl_interp_polynomial},
+        {InterpType::cspline, gsl_interp_cspline},
+        {InterpType::akima, gsl_interp_akima}};
+}  // namespace gauss2d::fit
+
+#endif  // GAUSS2D_FIT_HAS_GSL
+#endif  // GAUSS2DFIT_GSL_H

include/gslinterpolator.h

@@ -0,0 +1,54 @@
+#ifndef GAUSS2D_FIT_GSLINTERPOLATOR_H
+#define GAUSS2D_FIT_GSLINTERPOLATOR_H
+
+#ifdef GAUSS2D_FIT_HAS_GSL
+
+#include "gsl.h"
+#include "interpolation.h"
+
+#include "gauss2d/object.h"
+
+#include <gsl/gsl_errno.h>
+#include <gsl/gsl_spline.h>
+
+namespace gauss2d::fit {
+
+/**
+ * A wrapper for GSL 1D interpolators.
+ */
+class GSLInterpolator : public Object {
+private:
+    gsl_interp_accel* _acc;
+    gsl_spline* _spline;
+    size_t _n_knots;
+    std::vector<double> _x;
+    std::vector<double> _y;
+
+public:
+    const InterpType interp_type;
+    static constexpr InterpType INTERPTYPE_DEFAULT = InterpType::cspline;
+
+    /// Get the interpolant value for a knot of the given index
+    double get_knot_x(size_t idx) const;
+    /// Get the interpolated function value for a knot of the given index
+    double get_knot_y(size_t idx) const;
+
+    /// Get the interpolated function value for a given interpolant value
+    double eval(double x) const;
+    /// Get the derivative of the interpolated function value for a given interpolant value
+    double eval_deriv(double x) const;
+    /// Get the number of knots
+    size_t size() const;
+
+    std::string repr(bool name_keywords = false) const override;
+    std::string str() const override;
+
+    explicit GSLInterpolator(std::vector<double> x, std::vector<double> y,
+                             InterpType interp_type = INTERPTYPE_DEFAULT);
+    ~GSLInterpolator();
+};
+
+}  // namespace gauss2d::fit
+
+#endif  // GAUSS2D_FIT_HAS_GSL
+#endif  // GAUSS2D_FIT_GSLINTERPOLATOR_H

include/gslsersicmixinterpolator.h

@@ -0,0 +1,57 @@
+#ifndef GAUSS2D_FIT_GSLSERSICMIXINTERPOLATOR_H
+#define GAUSS2D_FIT_GSLSERSICMIXINTERPOLATOR_H
+
+#ifdef GAUSS2D_FIT_HAS_GSL
+
+#include <memory>
+
+#include "gsl.h"
+#include "gslinterpolator.h"
+#include "sersicmix.h"
+
+#include <gsl/gsl_errno.h>
+#include <gsl/gsl_spline.h>
+
+namespace gauss2d::fit {
+
+/**
+ * A SersicMixInterpolator that uses GSL interpolators between knots.
+ */
+class GSLSersicMixInterpolator : public SersicMixInterpolator {
+private:
+    mutable double _final_correction = 1.;
+    const InterpType _interp_type;
+    const unsigned short _order;
+    std::vector<std::pair<std::unique_ptr<GSLInterpolator>, std::unique_ptr<GSLInterpolator>>> _interps;
+
+public:
+    bool correct_final_integral = true;
+    static constexpr InterpType INTERPTYPE_DEFAULT = InterpType::cspline;
+    /// The knot positions and values.
+    const std::vector<SersicMixValues>& knots;
+
+    /// Get the multiplicative factor required to adjust the integral for the final order
+    /// component such that the sum of all integral factors is unity (normalized).
+    double get_final_correction() const;
+
+    std::vector<IntegralSize> get_integralsizes(double sersicindex) const override;
+    std::vector<IntegralSize> get_integralsizes_derivs(double sersicindex) const override;
+
+    InterpType get_interptype() const override;
+    unsigned short get_order() const override;
+
+    const double sersicindex_min;
+    const double sersicindex_max;
+
+    std::string repr(bool name_keywords = false) const override;
+    std::string str() const override;
+
+    explicit GSLSersicMixInterpolator(unsigned short order = SERSICMIX_ORDER_DEFAULT,
+                                      InterpType interp_type = INTERPTYPE_DEFAULT);
+    ~GSLSersicMixInterpolator();
+};
+
+}  // namespace gauss2d::fit
+
+#endif  // GAUSS2D_FIT_HAS_GSL
+#endif  // GAUSS2D_FIT_GSLSERSICMIXINTERPOLATOR_H

include/interpolation.h

@@ -0,0 +1,15 @@
+#ifndef GAUSS2DFIT_INTERPOLATION_H
+#define GAUSS2DFIT_INTERPOLATION_H
+
+namespace gauss2d::fit {
+
+enum class InterpType {
+    linear,      ///< Linear interpolation.
+    polynomial,  ///< Polynomial interpolation.
+    cspline,     ///< Cubic spline interpolation.
+    akima,       ///< Akima spline with natural boundary conditions.
+};
+
+}  // namespace gauss2d::fit
+
+#endif  // GAUSS2DFIT_INTERPOLATION_H

include/linearsersicmixinterpolator.h

@@ -19,6 +19,7 @@ class LinearSersicMixInterpolator : public SersicMixInterpolator {
     std::vector<IntegralSize> get_integralsizes(double sersicindex) const override;
     std::vector<IntegralSize> get_integralsizes_derivs(double sersicindex) const override;
 
+    InterpType get_interptype() const override;
     unsigned short get_order() const override;
 
     const double sersicindex_min;

include/meson.build

@@ -11,7 +11,11 @@ headers = [
     'gaussianmodelintegral.h',
     'gaussianparametricellipse.h',
     'gaussianprior.h',
+    'gsl.h',
+    'gslinterpolator.h',
+    'gslsersicmixinterpolator.h',
     'integralmodel.h',
+    'interpolation.h',
     'linearintegralmodel.h',
     'linearsersicmixinterpolator.h',
     'math.h',

include/model.h

@@ -8,10 +8,9 @@
 #include <string>
 
 #include "gauss2d/evaluate.h"
+#include "gauss2d/image.h"
 
 #include "data.h"
-#include "gauss2d/image.h"
-// TODO: Remove when DM-40674 fixed
 #include "parameters.h"
 #include "parametricmodel.h"
 #include "param_filter.h"
@@ -991,10 +990,8 @@ class Model : public ParametricModel {
         for (size_t idx = 0; idx < _size; ++idx) {
             result[idx] = this->_evaluate_observation(idx, print, is_loglike_grad);
         }
-        if(
-            (this->_mode == EvaluatorMode::loglike) || is_loglike_grad
-            || (this->_mode == EvaluatorMode::loglike_image) || (this->_mode == EvaluatorMode::jacobian)
-        ) {
+        if ((this->_mode == EvaluatorMode::loglike) || is_loglike_grad
+            || (this->_mode == EvaluatorMode::loglike_image) || (this->_mode == EvaluatorMode::jacobian)) {
             result[_size] = this->_evaluate_priors(print, normalize_loglike);
         }
 
@@ -1418,8 +1415,8 @@ class Model : public ParametricModel {
                 const auto& grad = grads[idx_param];
 
                 const double value_init = param.get_value();
-                const double value = param.get_value_transformed();
-                double diff = value * findiff_frac;
+                const double value_transformed = param.get_value_transformed();
+                double diff = value_transformed * findiff_frac;
                 if (std::abs(diff) < findiff_add) diff = findiff_add;
                 diff = finite_difference_param(param, diff);
 
@@ -1448,7 +1445,7 @@ class Model : public ParametricModel {
                 const double value_new = param.get_value();
                 if (value_new != value_init) {
                     throw std::logic_error("Could not return param=" + param.str()
-                                           + to_string_float(value_new) + "; diff="
+                                           + " to value_init=" + to_string_float(value_init) + "; diff="
                                            + to_string_float(value_new - value_init) + "); check limits");
                 }
                 if (n_failed > 0) {
@@ -1486,8 +1483,9 @@ class Model : public ParametricModel {
                 }
                 auto& param = (*found).get();
 
-                const double value = param.get_value_transformed();
-                double delta = value * findiff_frac;
+                const double value_init = param.get_value();
+                const double value_transformed = param.get_value_transformed();
+                double delta = value_transformed * findiff_frac;
                 if (std::abs(delta) < findiff_add) delta = findiff_add;
                 delta = finite_difference_param(param, delta);
                 auto result_new = prior->evaluate(true);
@@ -1508,10 +1506,12 @@ class Model : public ParametricModel {
                     }
                     idx_value++;
                 }
-                param.set_value_transformed(value);
-                if (param.get_value_transformed() != value) {
+                param.set_value(value_init);
+                const double value_new = param.get_value();
+                if (value_new != value_init) {
                     throw std::logic_error("Could not return param=" + param.str()
-                                           + " to original value=" + std::to_string(value));
+                                           + " to value_init=" + to_string_float(value_init) + "; diff="
+                                           + to_string_float(value_new - value_init) + "); check limits");
                 }
                 if (n_failed > 0) {
                     std::sort(ratios.begin(), ratios.end());

include/sersicmix.h

@@ -2,14 +2,16 @@
 #define GAUSS2D_FIT_SERSICMIX_H
 
 #include <array>
+#include <memory>
 #include <stdexcept>
 #include <vector>
 
+#include "interpolation.h"
 #include "gauss2d/object.h"
 
 namespace gauss2d::fit {
 
-static const short SERSICMIX_ORDER_DEFAULT = 4;
+static const unsigned short SERSICMIX_ORDER_DEFAULT = 4;
 
 /**
  * A pair of integral - size values for a Gaussian (sub)Component.
@@ -38,9 +40,13 @@ class SersicMixInterpolator : public Object {
     virtual std::vector<IntegralSize> get_integralsizes(double sersicindex) const = 0;
     virtual std::vector<IntegralSize> get_integralsizes_derivs(double sersicindex) const = 0;
 
+    virtual InterpType get_interptype() const = 0;
     virtual unsigned short get_order() const = 0;
 };
 
+const std::shared_ptr<const SersicMixInterpolator> get_sersic_mix_interpolator_default(
+        unsigned short order = SERSICMIX_ORDER_DEFAULT);
+
 /**
  * A vector of IntegralSize values for a given Sersic index.
  */

include/sersicmixcomponent.h

@@ -29,6 +29,8 @@ class SersicMixComponentIndexParameter : public SersicIndexParameter {
     double get_integralratio(unsigned short index) const;
     /// Return the integral ratio derivative for a given Gaussian sub-component index
     double get_integralratio_deriv(unsigned short index) const;
+    static std::shared_ptr<const SersicMixInterpolator> get_interpolator_default(unsigned short order
+                                                                                 = SERSICMIX_ORDER_DEFAULT);
     const parameters::Limits<double>& get_limits_maximal() const override;
     double get_min() const override { return 0.5; }
     double get_max() const override { return 8.0; }
@@ -37,7 +39,8 @@ class SersicMixComponentIndexParameter : public SersicIndexParameter {
     /// Return the size ratio derivative for a given Gaussian sub-component index
     double get_sizeratio_deriv(unsigned short index) const;
 
-    unsigned short order;
+    InterpType get_interptype() const;
+    unsigned short get_order() const;
 
     void set_value(double value) override;
     void set_value_transformed(double value_transformed) override;

meson.build

@@ -14,6 +14,10 @@ endif
 
 gsl_dep = dependency('gsl', required: false)
 
+if gsl_dep.found()
+    add_project_arguments('-DGAUSS2D_FIT_HAS_GSL', language : 'cpp')
+endif
+
 if use_eups
     gauss2d_eupsdir = run_command('eups', 'list', '-d', 'gauss2d', check: true).stdout().strip()
     cpp = meson.get_compiler('cpp')

python/gauss2d/fit/gsl.cc

@@ -0,0 +1,43 @@
+/*
+ * This file is part of gauss2dfit.
+ *
+ * Developed for the LSST Data Management System.
+ * This product includes software developed by the LSST Project
+ * (https://www.lsst.org).
+ * See the COPYRIGHT file at the top-level directory of this distribution
+ * for details of code ownership.
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <https://www.gnu.org/licenses/>.
+ */
+
+#ifdef GAUSS2D_FIT_HAS_GSL
+
+#include <pybind11/attr.h>
+#include <pybind11/pybind11.h>
+#include <pybind11/stl.h>
+
+#include "pybind11.h"
+
+#include "gauss2d/fit/gsl.h"
+
+namespace py = pybind11;
+using namespace pybind11::literals;
+
+namespace g2f = gauss2d::fit;
+
+void bind_gsl(py::module &m) {
+    // Placeholder for now - Python doesn't need to know about gsl_interp_type
+}
+
+#endif  // GAUSS2D_FIT_HAS_GSL