Added scurve fitting (#168)

- added scurve fitting with two different signs (scurve, scurve2)
- at the moment no option to set initial parameters

---------

Co-authored-by: JulianHeymes <julian.heymes@psi.ch>

include/aare/Fit.hpp

@@ -15,6 +15,12 @@ NDArray<double, 1> gaus(NDView<double, 1> x, NDView<double, 1> par);
 double pol1(const double x, const double *par);
 NDArray<double, 1> pol1(NDView<double, 1> x, NDView<double, 1> par);
 
+double scurve(const double x, const double *par);
+NDArray<double, 1> scurve(NDView<double, 1> x, NDView<double, 1> par);
+
+double scurve2(const double x, const double *par);
+NDArray<double, 1> scurve2(NDView<double, 1> x, NDView<double, 1> par);
+
 } // namespace func
 
 
@@ -25,6 +31,9 @@ std::array<double, 3> gaus_init_par(const NDView<double, 1> x, const NDView<doub
 
 std::array<double, 2> pol1_init_par(const NDView<double, 1> x, const NDView<double, 1> y);
 
+std::array<double, 6> scurve_init_par(const NDView<double, 1> x, const NDView<double, 1> y); 
+std::array<double, 6> scurve2_init_par(const NDView<double, 1> x, const NDView<double, 1> y);
+
 static constexpr int DEFAULT_NUM_THREADS = 4;
 
 /**
@@ -38,7 +47,7 @@ NDArray<double, 1> fit_gaus(NDView<double, 1> x, NDView<double, 1> y);
 /**
  * @brief Fit a 1D Gaussian to each pixel. Data layout [row, col, values]
  * @param x x values
- * @param y y vales, layout [row, col, values]
+ * @param y y values, layout [row, col, values]
  * @param n_threads number of threads to use
  */
 
@@ -51,7 +60,7 @@ NDArray<double, 3> fit_gaus(NDView<double, 1> x, NDView<double, 3> y,
 /**
  * @brief Fit a 1D Gaussian with error estimates
  * @param x x values
- * @param y y vales, layout [row, col, values]
+ * @param y y values, layout [row, col, values]
  * @param y_err error in y, layout [row, col, values]
  * @param par_out output parameters
  * @param par_err_out output error parameters
@@ -64,7 +73,7 @@ void fit_gaus(NDView<double, 1> x, NDView<double, 1> y, NDView<double, 1> y_err,
  * @brief Fit a 1D Gaussian to each pixel with error estimates. Data layout
  * [row, col, values]
  * @param x x values
- * @param y y vales, layout [row, col, values]
+ * @param y y values, layout [row, col, values]
  * @param y_err error in y, layout [row, col, values]
  * @param par_out output parameters, layout [row, col, values]
  * @param par_err_out output parameter errors, layout [row, col, values]
@@ -88,5 +97,19 @@ void fit_pol1(NDView<double, 1> x, NDView<double, 3> y, NDView<double, 3> y_err,
               NDView<double, 3> par_out, NDView<double, 3> par_err_out,NDView<double, 2> chi2_out,
               int n_threads = DEFAULT_NUM_THREADS);
 
+NDArray<double, 1> fit_scurve(NDView<double, 1> x, NDView<double, 1> y); 
+NDArray<double, 3> fit_scurve(NDView<double, 1> x, NDView<double, 3> y, int n_threads); 
+void fit_scurve(NDView<double, 1> x, NDView<double, 1> y, NDView<double, 1> y_err,
+              NDView<double, 1> par_out, NDView<double, 1> par_err_out, double& chi2); 
+void fit_scurve(NDView<double, 1> x, NDView<double, 3> y, NDView<double, 3> y_err,
+              NDView<double, 3> par_out, NDView<double, 3> par_err_out, NDView<double, 2> chi2_out,
+              int n_threads); 
 
+NDArray<double, 1> fit_scurve2(NDView<double, 1> x, NDView<double, 1> y); 
+NDArray<double, 3> fit_scurve2(NDView<double, 1> x, NDView<double, 3> y, int n_threads); 
+void fit_scurve2(NDView<double, 1> x, NDView<double, 1> y, NDView<double, 1> y_err,
+              NDView<double, 1> par_out, NDView<double, 1> par_err_out, double& chi2); 
+void fit_scurve2(NDView<double, 1> x, NDView<double, 3> y, NDView<double, 3> y_err,
+              NDView<double, 3> par_out, NDView<double, 3> par_err_out, NDView<double, 2> chi2_out,
+              int n_threads); 
 } // namespace aare

python/aare/__init__.py

@@ -15,7 +15,7 @@ from .ClusterFinder import ClusterFinder, ClusterCollector, ClusterFinderMT, Clu
 from .ClusterVector import ClusterVector
 
 
-from ._aare import fit_gaus, fit_pol1
+from ._aare import fit_gaus, fit_pol1, fit_scurve, fit_scurve2
 from ._aare import Interpolator
 from ._aare import calculate_eta2
 

python/aare/func.py

@@ -1 +1 @@
-from ._aare import gaus, pol1
+from ._aare import gaus, pol1, scurve, scurve2

python/src/fit.hpp

@@ -55,6 +55,47 @@ void define_fit_bindings(py::module &m) {
         )",
         py::arg("x"), py::arg("par"));
 
+    m.def(
+        "scurve",
+        [](py::array_t<double, py::array::c_style | py::array::forcecast> x,
+           py::array_t<double, py::array::c_style | py::array::forcecast> par) {
+            auto x_view = make_view_1d(x);
+            auto par_view = make_view_1d(par);
+            auto y = new NDArray<double, 1>{aare::func::scurve(x_view, par_view)};
+            return return_image_data(y);
+        },
+        R"(
+        Evaluate a 1D scurve function for all points in x using parameters par. 
+        
+        Parameters
+        ----------
+        x : array_like
+            The points at which to evaluate the scurve function.
+        par : array_like
+            The parameters of the scurve function. The first element is the background slope, the second element is the background intercept, the third element is the mean, the fourth element is the standard deviation, the fifth element is inflexion point count number, and the sixth element is C.     
+        )",
+        py::arg("x"), py::arg("par"));
+
+    m.def(
+        "scurve2",
+        [](py::array_t<double, py::array::c_style | py::array::forcecast> x,
+           py::array_t<double, py::array::c_style | py::array::forcecast> par) {
+            auto x_view = make_view_1d(x);
+            auto par_view = make_view_1d(par);
+            auto y = new NDArray<double, 1>{aare::func::scurve2(x_view, par_view)};
+            return return_image_data(y);
+        },
+        R"(
+        Evaluate a 1D scurve2 function for all points in x using parameters par. 
+        
+        Parameters
+        ----------
+        x : array_like
+            The points at which to evaluate the scurve function.
+        par : array_like
+            The parameters of the scurve2 function. The first element is the background slope, the second element is the background intercept, the third element is the mean, the fourth element is the standard deviation, the fifth element is inflexion point count number, and the sixth element is C.     
+        )",
+        py::arg("x"), py::arg("par"));
 
     m.def(
         "fit_gaus",
@@ -235,6 +276,180 @@ n_threads : int, optional
         R"(
 Fit a 1D polynomial to data with error estimates.
 
+Parameters
+----------
+x : array_like
+    The x values.
+y : array_like
+    The y values.
+y_err : array_like
+    The error in the y values.
+n_threads : int, optional
+    The number of threads to use. Default is 4.
+)",
+        py::arg("x"), py::arg("y"), py::arg("y_err"), py::arg("n_threads") = 4);
+
+//=========
+    m.def(
+        "fit_scurve",
+        [](py::array_t<double, py::array::c_style | py::array::forcecast> x,
+           py::array_t<double, py::array::c_style | py::array::forcecast> y,
+           int n_threads) {
+            if (y.ndim() == 3) {
+                auto par = new NDArray<double, 3>{};
+
+                auto x_view = make_view_1d(x);
+                auto y_view = make_view_3d(y);
+                *par = aare::fit_scurve(x_view, y_view, n_threads);
+                return return_image_data(par);
+            } else if (y.ndim() == 1) {
+                auto par = new NDArray<double, 1>{};
+                auto x_view = make_view_1d(x);
+                auto y_view = make_view_1d(y);
+                *par = aare::fit_scurve(x_view, y_view);
+                return return_image_data(par);
+            } else {
+                throw std::runtime_error("Data must be 1D or 3D");
+            }
+        },
+        py::arg("x"), py::arg("y"), py::arg("n_threads") = 4);
+
+    m.def(
+        "fit_scurve",
+        [](py::array_t<double, py::array::c_style | py::array::forcecast> x,
+           py::array_t<double, py::array::c_style | py::array::forcecast> y,
+           py::array_t<double, py::array::c_style | py::array::forcecast> y_err,
+           int n_threads) {
+            if (y.ndim() == 3) {
+                auto par = new NDArray<double, 3>({y.shape(0), y.shape(1), 6});
+
+                auto par_err =
+                    new NDArray<double, 3>({y.shape(0), y.shape(1), 6});
+
+                auto y_view = make_view_3d(y);
+                auto y_view_err = make_view_3d(y_err);
+                auto x_view = make_view_1d(x);
+
+                auto chi2 = new NDArray<double, 2>({y.shape(0), y.shape(1)});
+
+                aare::fit_scurve(x_view, y_view, y_view_err, par->view(),
+                               par_err->view(), chi2->view(), n_threads);
+                return py::dict("par"_a = return_image_data(par),
+                                "par_err"_a = return_image_data(par_err), 
+                                "chi2"_a = return_image_data(chi2),
+                                "Ndf"_a = y.shape(2) - 2);
+
+
+            } else if (y.ndim() == 1) {
+                auto par = new NDArray<double, 1>({2});
+                auto par_err = new NDArray<double, 1>({2});
+
+                auto y_view = make_view_1d(y);
+                auto y_view_err = make_view_1d(y_err);
+                auto x_view = make_view_1d(x);
+
+                double chi2 = 0;
+
+                aare::fit_scurve(x_view, y_view, y_view_err, par->view(),
+                               par_err->view(), chi2);
+                return py::dict("par"_a = return_image_data(par),
+                                "par_err"_a = return_image_data(par_err),
+                                "chi2"_a = chi2, "Ndf"_a = y.size() - 2);
+
+            } else {
+                throw std::runtime_error("Data must be 1D or 3D");
+            }
+        },
+        R"(
+Fit a 1D polynomial to data with error estimates.
+
+Parameters
+----------
+x : array_like
+    The x values.
+y : array_like
+    The y values.
+y_err : array_like
+    The error in the y values.
+n_threads : int, optional
+    The number of threads to use. Default is 4.
+)",
+        py::arg("x"), py::arg("y"), py::arg("y_err"), py::arg("n_threads") = 4);
+    
+
+    m.def(
+        "fit_scurve2",
+        [](py::array_t<double, py::array::c_style | py::array::forcecast> x,
+           py::array_t<double, py::array::c_style | py::array::forcecast> y,
+           int n_threads) {
+            if (y.ndim() == 3) {
+                auto par = new NDArray<double, 3>{};
+
+                auto x_view = make_view_1d(x);
+                auto y_view = make_view_3d(y);
+                *par = aare::fit_scurve2(x_view, y_view, n_threads);
+                return return_image_data(par);
+            } else if (y.ndim() == 1) {
+                auto par = new NDArray<double, 1>{};
+                auto x_view = make_view_1d(x);
+                auto y_view = make_view_1d(y);
+                *par = aare::fit_scurve2(x_view, y_view);
+                return return_image_data(par);
+            } else {
+                throw std::runtime_error("Data must be 1D or 3D");
+            }
+        },
+        py::arg("x"), py::arg("y"), py::arg("n_threads") = 4);
+
+    m.def(
+        "fit_scurve2",
+        [](py::array_t<double, py::array::c_style | py::array::forcecast> x,
+           py::array_t<double, py::array::c_style | py::array::forcecast> y,
+           py::array_t<double, py::array::c_style | py::array::forcecast> y_err,
+           int n_threads) {
+            if (y.ndim() == 3) {
+                auto par = new NDArray<double, 3>({y.shape(0), y.shape(1), 6});
+
+                auto par_err =
+                    new NDArray<double, 3>({y.shape(0), y.shape(1), 6});
+
+                auto y_view = make_view_3d(y);
+                auto y_view_err = make_view_3d(y_err);
+                auto x_view = make_view_1d(x);
+
+                auto chi2 = new NDArray<double, 2>({y.shape(0), y.shape(1)});
+
+                aare::fit_scurve2(x_view, y_view, y_view_err, par->view(),
+                               par_err->view(), chi2->view(), n_threads);
+                return py::dict("par"_a = return_image_data(par),
+                                "par_err"_a = return_image_data(par_err), 
+                                "chi2"_a = return_image_data(chi2),
+                                "Ndf"_a = y.shape(2) - 2);
+
+
+            } else if (y.ndim() == 1) {
+                auto par = new NDArray<double, 1>({6});
+                auto par_err = new NDArray<double, 1>({6});
+
+                auto y_view = make_view_1d(y);
+                auto y_view_err = make_view_1d(y_err);
+                auto x_view = make_view_1d(x);
+
+                double chi2 = 0;
+
+                aare::fit_scurve2(x_view, y_view, y_view_err, par->view(),
+                               par_err->view(), chi2);
+                return py::dict("par"_a = return_image_data(par),
+                                "par_err"_a = return_image_data(par_err),
+                                "chi2"_a = chi2, "Ndf"_a = y.size() - 2);
+
+            } else {
+                throw std::runtime_error("Data must be 1D or 3D");
+            }
+        },
+        R"(
+Fit a 1D polynomial to data with error estimates.
+
 Parameters
 ----------
 x : array_like

src/Fit.cpp

@@ -34,6 +34,30 @@ NDArray<double, 1> pol1(NDView<double, 1> x, NDView<double, 1> par) {
     return y;
 }
 
+double scurve(const double x, const double * par) {
+    return (par[0] + par[1] * x) + 0.5 * (1 + erf((x - par[2]) / (sqrt(2) * par[3]))) * (par[4] + par[5] * (x - par[2]));
+}
+
+NDArray<double, 1> scurve(NDView<double, 1> x, NDView<double, 1> par) {
+    NDArray<double, 1> y({x.shape()}, 0);
+    for (ssize_t i = 0; i < x.size(); i++) {
+        y(i) = scurve(x(i), par.data());
+    }
+    return y;
+}
+
+double scurve2(const double x, const double * par) {
+    return (par[0] + par[1] * x) + 0.5 * (1 - erf((x - par[2]) / (sqrt(2) * par[3]))) * (par[4] + par[5] * (x - par[2]));
+}
+
+NDArray<double, 1> scurve2(NDView<double, 1> x, NDView<double, 1> par) {
+    NDArray<double, 1> y({x.shape()}, 0);
+    for (ssize_t i = 0; i < x.size(); i++) {
+        y(i) = scurve2(x(i), par.data());
+    }
+    return y;
+}
+
 } // namespace func
 
 NDArray<double, 1> fit_gaus(NDView<double, 1> x, NDView<double, 1> y) {
@@ -273,4 +297,229 @@ NDArray<double, 3> fit_pol1(NDView<double, 1> x, NDView<double, 3> y,
     return result;
 }
 
+// ~~ S-CURVES ~~
+
+// SCURVE --
+std::array<double, 6> scurve_init_par(const NDView<double, 1> x, const NDView<double, 1> y){
+        // Estimate the initial parameters for the fit
+        std::array<double, 6> start_par{0, 0, 0, 0, 0, 0};
+
+        auto ymax = std::max_element(y.begin(), y.end());
+        auto ymin = std::min_element(y.begin(), y.end());
+        start_par[4] = *ymin + (*ymax - *ymin) / 2;
+        
+        // Find the first x where the corresponding y value is above the threshold (start_par[4])
+        for (ssize_t i = 0; i < y.size(); ++i) {
+            if (y[i] >= start_par[4]) {
+                start_par[2] = x[i];
+                break; // Exit the loop after finding the first valid x
+            }
+        }
+
+        start_par[3] = 2 * sqrt(start_par[2]); 
+        start_par[0] = 100; 
+        start_par[1] = 0.25;
+        start_par[5] = 1;
+        return start_par;
+}
+
+// - No error
+NDArray<double, 1> fit_scurve(NDView<double, 1> x, NDView<double, 1> y) {
+    NDArray<double, 1> result = scurve_init_par(x, y);
+    lm_status_struct status;
+
+    lmcurve(result.size(), result.data(), x.size(), x.data(), y.data(),
+            aare::func::scurve, &lm_control_double, &status);
+
+    return result;
+}
+
+NDArray<double, 3> fit_scurve(NDView<double, 1> x, NDView<double, 3> y, int n_threads) {
+    NDArray<double, 3> result({y.shape(0), y.shape(1), 6}, 0);
+
+    auto process = [&x, &y, &result](ssize_t first_row, ssize_t last_row) {
+        for (ssize_t row = first_row; row < last_row; row++) {
+            for (ssize_t col = 0; col < y.shape(1); col++) {
+                NDView<double, 1> values(&y(row, col, 0), {y.shape(2)});
+                auto res = fit_scurve(x, values);
+                result(row, col, 0) = res(0);
+                result(row, col, 1) = res(1);
+                result(row, col, 2) = res(2);
+                result(row, col, 3) = res(3);
+                result(row, col, 4) = res(4);
+                result(row, col, 5) = res(5);
+            }
+        }
+    };
+
+    auto tasks = split_task(0, y.shape(0), n_threads);
+    RunInParallel(process, tasks);
+    return result;
+}
+
+// - Error
+void fit_scurve(NDView<double, 1> x, NDView<double, 1> y, NDView<double, 1> y_err,
+              NDView<double, 1> par_out, NDView<double, 1> par_err_out, double& chi2) {
+
+    // Check that we have the correct sizes
+    if (y.size() != x.size() || y.size() != y_err.size() ||
+        par_out.size() != 6 || par_err_out.size() != 6) {
+        throw std::runtime_error("Data, x, data_err must have the same size "
+                                 "and par_out, par_err_out must have size 6");
+    }
+
+    lm_status_struct status;
+    par_out = scurve_init_par(x, y);
+    std::array<double, 36> cov = {0}; // size 6x6
+    // std::array<double, 4> cov{0, 0, 0, 0};
+
+    lmcurve2(par_out.size(), par_out.data(), par_err_out.data(), cov.data(),
+             x.size(), x.data(), y.data(), y_err.data(), aare::func::scurve,
+             &lm_control_double, &status);
+
+    // Calculate chi2
+    chi2 = 0;
+    for (ssize_t i = 0; i < y.size(); i++) {
+        chi2 += std::pow((y(i) - func::pol1(x(i), par_out.data())) / y_err(i), 2);
+    }
+}
+
+void fit_scurve(NDView<double, 1> x, NDView<double, 3> y, NDView<double, 3> y_err,
+              NDView<double, 3> par_out, NDView<double, 3> par_err_out, NDView<double, 2> chi2_out,
+              int n_threads) {
+
+    auto process = [&](ssize_t first_row, ssize_t last_row) {
+        for (ssize_t row = first_row; row < last_row; row++) {
+            for (ssize_t col = 0; col < y.shape(1); col++) {
+                NDView<double, 1> y_view(&y(row, col, 0), {y.shape(2)});
+                NDView<double, 1> y_err_view(&y_err(row, col, 0),
+                                             {y_err.shape(2)});
+                NDView<double, 1> par_out_view(&par_out(row, col, 0),
+                                               {par_out.shape(2)});
+                NDView<double, 1> par_err_out_view(&par_err_out(row, col, 0),
+                                                   {par_err_out.shape(2)});
+
+                fit_scurve(x, y_view, y_err_view, par_out_view, par_err_out_view, chi2_out(row, col));
+
+            }
+        }
+    };
+
+    auto tasks = split_task(0, y.shape(0), n_threads);
+    RunInParallel(process, tasks);
+
+}
+
+// SCURVE2 ---
+
+std::array<double, 6> scurve2_init_par(const NDView<double, 1> x, const NDView<double, 1> y){
+        // Estimate the initial parameters for the fit
+        std::array<double, 6> start_par{0, 0, 0, 0, 0, 0};
+
+        auto ymax = std::max_element(y.begin(), y.end());
+        auto ymin = std::min_element(y.begin(), y.end());
+        start_par[4] = *ymin + (*ymax - *ymin) / 2;
+        
+        // Find the first x where the corresponding y value is above the threshold (start_par[4])
+        for (ssize_t i = 0; i < y.size(); ++i) {
+            if (y[i] <= start_par[4]) {
+                start_par[2] = x[i];
+                break; // Exit the loop after finding the first valid x
+            }
+        }
+
+        start_par[3] = 2 * sqrt(start_par[2]); 
+        start_par[0] = 100; 
+        start_par[1] = 0.25;
+        start_par[5] = -1;
+        return start_par;
+}
+
+// - No error
+NDArray<double, 1> fit_scurve2(NDView<double, 1> x, NDView<double, 1> y) {
+    NDArray<double, 1> result = scurve2_init_par(x, y);
+    lm_status_struct status;
+
+    lmcurve(result.size(), result.data(), x.size(), x.data(), y.data(),
+            aare::func::scurve2, &lm_control_double, &status);
+
+    return result;
+}
+
+NDArray<double, 3> fit_scurve2(NDView<double, 1> x, NDView<double, 3> y, int n_threads) {
+    NDArray<double, 3> result({y.shape(0), y.shape(1), 6}, 0);
+
+    auto process = [&x, &y, &result](ssize_t first_row, ssize_t last_row) {
+        for (ssize_t row = first_row; row < last_row; row++) {
+            for (ssize_t col = 0; col < y.shape(1); col++) {
+                NDView<double, 1> values(&y(row, col, 0), {y.shape(2)});
+                auto res = fit_scurve2(x, values);
+                result(row, col, 0) = res(0);
+                result(row, col, 1) = res(1);
+                result(row, col, 2) = res(2);
+                result(row, col, 3) = res(3);
+                result(row, col, 4) = res(4);
+                result(row, col, 5) = res(5);
+            }
+        }
+    };
+
+    auto tasks = split_task(0, y.shape(0), n_threads);
+    RunInParallel(process, tasks);
+    return result;
+}
+
+// - Error
+void fit_scurve2(NDView<double, 1> x, NDView<double, 1> y, NDView<double, 1> y_err,
+              NDView<double, 1> par_out, NDView<double, 1> par_err_out, double& chi2) {
+
+    // Check that we have the correct sizes
+    if (y.size() != x.size() || y.size() != y_err.size() ||
+        par_out.size() != 6 || par_err_out.size() != 6) {
+        throw std::runtime_error("Data, x, data_err must have the same size "
+                                 "and par_out, par_err_out must have size 6");
+    }
+
+    lm_status_struct status;
+    par_out = scurve2_init_par(x, y);
+    std::array<double, 36> cov = {0}; // size 6x6
+    // std::array<double, 4> cov{0, 0, 0, 0};
+
+    lmcurve2(par_out.size(), par_out.data(), par_err_out.data(), cov.data(),
+             x.size(), x.data(), y.data(), y_err.data(), aare::func::scurve2,
+             &lm_control_double, &status);
+
+    // Calculate chi2
+    chi2 = 0;
+    for (ssize_t i = 0; i < y.size(); i++) {
+        chi2 += std::pow((y(i) - func::pol1(x(i), par_out.data())) / y_err(i), 2);
+    }
+}
+
+void fit_scurve2(NDView<double, 1> x, NDView<double, 3> y, NDView<double, 3> y_err,
+              NDView<double, 3> par_out, NDView<double, 3> par_err_out, NDView<double, 2> chi2_out,
+              int n_threads) {
+
+    auto process = [&](ssize_t first_row, ssize_t last_row) {
+        for (ssize_t row = first_row; row < last_row; row++) {
+            for (ssize_t col = 0; col < y.shape(1); col++) {
+                NDView<double, 1> y_view(&y(row, col, 0), {y.shape(2)});
+                NDView<double, 1> y_err_view(&y_err(row, col, 0),
+                                             {y_err.shape(2)});
+                NDView<double, 1> par_out_view(&par_out(row, col, 0),
+                                               {par_out.shape(2)});
+                NDView<double, 1> par_err_out_view(&par_err_out(row, col, 0),
+                                                   {par_err_out.shape(2)});
+
+                fit_scurve2(x, y_view, y_err_view, par_out_view, par_err_out_view, chi2_out(row, col));
+
+            }
+        }
+    };
+
+    auto tasks = split_task(0, y.shape(0), n_threads);
+    RunInParallel(process, tasks);
+
+}
+
 } // namespace aare