Fixing parameters of a fitting function in Nonlinear Least-Square GSL

Question

I'm working on some code that I'm writing which uses the [GNU Scientific Library (GSL)][1]'s Nonlinear least-squares algorithm for curve fitting.

I have been successful in obtaining a working code that estimate the right parameters from the fitting analysis using a C++ wrapper from https://github.com/Eleobert/gsl-curve-fit/blob/master/example.cpp.

Now, I would like to fix some of the parameters of the function to be fit. And I would like to modify the function in such a way that I can already input the value of the parameter to be fixed.

Any idea on how to do? I'm showing here the full code.

This is the code for performing nonlinear least-squares fitting:

#include <gsl/gsl_vector.h>
#include <gsl/gsl_multifit_nlinear.h>
#include <iostream>
#include <random>
#include <vector>
#include <cassert>
#include <functional>

template <typename F, size_t... Is>
auto gen_tuple_impl(F func, std::index_sequence<Is...> )
{
    return std::make_tuple(func(Is)...);
}

template <size_t N, typename F>
auto gen_tuple(F func)
{
    return gen_tuple_impl(func, std::make_index_sequence<N>{} );
}

auto internal_solve_system(gsl_vector* initial_params, gsl_multifit_nlinear_fdf *fdf,
                           gsl_multifit_nlinear_parameters *params) -> std::vector<double>
{
    // This specifies a trust region method
    const gsl_multifit_nlinear_type *T = gsl_multifit_nlinear_trust;
    const size_t max_iter = 200;
    const double xtol = 1.0e-8;
    const double gtol = 1.0e-8;
    const double ftol = 1.0e-8;

    auto *work = gsl_multifit_nlinear_alloc(T, params, fdf->n, fdf->p);
    int info;

    // initialize solver
    gsl_multifit_nlinear_init(initial_params, fdf, work);
    //iterate until convergence
    gsl_multifit_nlinear_driver(max_iter, xtol, gtol, ftol, nullptr, nullptr, &info, work);

    // result will be stored here
    gsl_vector * y    = gsl_multifit_nlinear_position(work);
    auto result = std::vector<double>(initial_params->size);

    for(int i = 0; i < result.size(); i++)
    {
        result[i] = gsl_vector_get(y, i);
    }

    auto niter = gsl_multifit_nlinear_niter(work);
    auto nfev  = fdf->nevalf;
    auto njev  = fdf->nevaldf;
    auto naev  = fdf->nevalfvv;

    // nfev - number of function evaluations
    // njev - number of Jacobian evaluations
    // naev - number of f_vv evaluations
    //logger::debug("curve fitted after ", niter, " iterations {nfev = ", nfev, "} {njev = ", njev, "} {naev = ", naev, "}");

    gsl_multifit_nlinear_free(work);
    gsl_vector_free(initial_params);
    return result;
}

auto internal_make_gsl_vector_ptr(const std::vector<double>& vec) -> gsl_vector*
{
    auto* result = gsl_vector_alloc(vec.size());
    int i = 0;
    for(const auto e: vec)
    {
        gsl_vector_set(result, i, e);
        i++;
    }
    return result;
}

template<typename C1>
struct fit_data
{
    const std::vector<double>& t;
    const std::vector<double>& y;
    // the actual function to be fitted
    C1 f;
};


template<typename FitData, int n_params>
int internal_f(const gsl_vector* x, void* params, gsl_vector *f)
{
    auto* d  = static_cast<FitData*>(params);
    // Convert the parameter values from gsl_vector (in x) into std::tuple
    auto init_args = [x](int index)
    {
        return gsl_vector_get(x, index);
    };
    auto parameters = gen_tuple<n_params>(init_args);

    // Calculate the error for each...
    for (size_t i = 0; i < d->t.size(); ++i)
    {
        double ti = d->t[i];
        double yi = d->y[i];
        auto func = [ti, &d](auto ...xs)
        {
            // call the actual function to be fitted
            return d->f(ti, xs...);
        };
        auto y = std::apply(func, parameters);
        gsl_vector_set(f, i, yi - y);
    }
    return GSL_SUCCESS;
}

using func_f_type   = int (*) (const gsl_vector*, void*, gsl_vector*);
using func_df_type  = int (*) (const gsl_vector*, void*, gsl_matrix*);
using func_fvv_type = int (*) (const gsl_vector*, const gsl_vector *, void *, gsl_vector *);


auto internal_make_gsl_vector_ptr(const std::vector<double>& vec) -> gsl_vector*;


auto internal_solve_system(gsl_vector* initial_params, gsl_multifit_nlinear_fdf *fdf,
                           gsl_multifit_nlinear_parameters *params) -> std::vector<double>;

template<typename C1>
auto curve_fit_impl(func_f_type f, func_df_type df, func_fvv_type fvv, gsl_vector* initial_params, fit_data<C1>& fd) -> std::vector<double>
{
    assert(fd.t.size() == fd.y.size());

    auto fdf = gsl_multifit_nlinear_fdf();
    auto fdf_params = gsl_multifit_nlinear_default_parameters();

    fdf.f   = f;
    fdf.df  = df;
    fdf.fvv = fvv;
    fdf.n   = fd.t.size();
    fdf.p   = initial_params->size;
    fdf.params = &fd;

    // "This selects the Levenberg-Marquardt algorithm with geodesic acceleration."
    fdf_params.trs = gsl_multifit_nlinear_trs_lmaccel;
    return internal_solve_system(initial_params, &fdf, &fdf_params);
}

template<typename Callable>
auto curve_fit(Callable f, const std::vector<double>& initial_params, const std::vector<double>& x, const std::vector<double>& y) -> std::vector<double>
{
    // We can't pass lambdas without convert to std::function.
    constexpr auto n = 3;
    assert(initial_params.size() == n);

    auto params = internal_make_gsl_vector_ptr(initial_params);
    auto fd = fit_data<Callable>{x, y, f};
    return curve_fit_impl(internal_f<decltype(fd), n>, nullptr, nullptr, params,  fd);
}

// linspace from https://github.com/Eleobert/meth/blob/master/interpolators.hpp
template <typename Container>
auto linspace(typename Container::value_type a, typename Container::value_type b, size_t n)
{
    assert(b > a);
    assert(n > 1);

    Container res(n);
    const auto step = (b - a) / (n - 1);
    auto val = a;
    for(auto& e: res)
    {
        e = val;
        val += step;
    }
    return res;
}

This is the function I use for fitting:

double gaussian(double x, double a, double b, double c)
{
    const double z = (x - b) / c;
    return a * std::exp(-0.5 * z * z);
}

And these last lines create a fake dataset of observed data (with some noise which is normally distributed) and test the fitting curve function.

int main()
{
    auto device = std::random_device();
    auto gen    = std::mt19937(device());

    auto xs = linspace<std::vector<double>>(0.0, 1.0, 300);
    auto ys = std::vector<double>(xs.size());

    double a = 5.0, b = 0.4, c = 0.15;

    for(size_t i = 0; i < xs.size(); i++)
    {
        auto y =  gaussian(xs[i], a, b, c);
        auto dist  = std::normal_distribution(0.0, 0.1 * y);
        ys[i] = y + dist(gen);
    }

    auto r = curve_fit(gaussian, {1.0, 0.0, 1.0}, xs, ys);

    std::cout << "result: " << r[0] << ' ' << r[1] << ' ' << r[2] << '\n';
    std::cout << "error : " << r[0] - a << ' ' << r[1] - b << ' ' << r[2] - c << '\n';

}

In this case, I would like to fix one of the a, b, c parameters and estimate the remaining two. For example, fix a and estimate b and c. But I would like to find a solution such that I can input any value to the fixed parameter a, without needing to modify the gaussian function every time.

Answer 1

Ok. Here's the answer based on the code linked in http://github.com/Eleobert/gsl-curve-fit/blob/master/example.cpp. However, this is not the code posted in the question: you should update your question accordingly so that others may take advantage from both the question & answer.

So, basically, the main problem is that GSL is a library written in pure C, whereas you use a high-level wrapper written in C++, published in the aforementioned link . While the wrapper is written pretty well in modern C++, it has one basic problem: it is "stiff" - it can be used only for a subclass of problems it was designed for, and this subclass is a rather narrow subset of the capabilities offered by the original C code.

Let's try to improve it a bit and start from how the wrapper is supposed to be used:

double gaussian(double x, double a, double b, double c)
{
    const double z = (x - b) / c;
    return a * std::exp(-0.5 * z * z);
}

int main()
{
    auto device = std::random_device();
    auto gen = std::mt19937(device());

    auto xs = linspace<std::vector<double>>(0.0, 1.0, 300);
    auto ys = std::vector<double>(xs.size());

    double a = 5.0, b = 0.4, c = 0.15;

    for (size_t i = 0; i < xs.size(); i++)
    {
        auto y = gaussian(xs[i], a, b, c);
        auto dist = std::normal_distribution(0.0, 0.1 * y);
        ys[i] = y + dist(gen);
    }

    auto result = curve_fit(gaussian, {1.0, 0.0, 1.0}, xs, ys);
   // use result
}

This code is amazingly simple compared to its original, C-language counterpart. One initializes the x-y pairs of values, here stored as vectors xs and ys and executes a single function that take 4 easy to understand parameters: the function to be fitted to the data, the initial values of the fitting parameters the function depends on, the x values and the corresponding y values of the data to which the function must be fitted.

Your problem is how to keep this high-level interface, but use it for fitting functions where only some parameters are "free", that is, can be changed during the fitting procedure, while the values of others must be fixed. This could be easily achieved using, e.g., global variables that the function has access to, but we hate global variables and never use them without a real cause.

I suggest using a well-known C++ alternative: functors. Look:

struct gaussian_fixed_a
{
    double a;
    gaussian_fixed_a(double a) : a{a} {}
    double operator()(double x, double b, double c) const { return gaussian(x, a, b, c); }
};

This struct/class introduces a new type of function objects. In the constructor, parameter a is passed and stored in an object. Then, there's a function call operator that takes only 3 parameters rather than 4, substituting a from its stored value. This object can pretend to be a Gaussian with a fixed constant and only other 3 arguments, x, b, and c, that can vary.

We would like to use it like this:

    gaussian_fixed_a g(a);
    auto r2 = curve_fit(g, std::array{0.444, 0.11}, xs, ys);

This is almost the same code you'd use for the original wrapper save for 2 differences:

1. You now use an object name (here: g) rather than a function name
2. You have to pass the number of arguments to curve_fit as a compile-time constant, as it is then used internally by it to call a template parametrized by this number. I achieve it by using std::array, for which the compiler can deduce its size at compile time, just as needed. An alternative would ba a nasty template syntax, curve_fit<2>(...

For this to work, you need to change the interface of curve_fit, from

template <typename Callable>
auto curve_fit(Callable f, const std::vector<double>& initial_params, const std::vector<double>& x,
               const std::vector<double>& y) -> std::vector<double>

to

template <typename Callable, auto n>
auto curve_fit(Callable f, const std::array<double, n>& initial_params, const std::vector<double>& x,
               const std::vector<double>& y) -> std::vector<double>

(btw: this -> syntax with well-known type on its right-hand side is not the best one, IMHO, but let it be). The idea is to force the compiler to read the number of fitting parameters at compile time form the size of the array.

Then you need to make a similar adjustment in the argument list of curve_fit_impl - and this is almost it.

Here I spent quite a lot of time trying to figure out why this code does not work. It turned out it had worked all the time, the secret is that if you fit a function to some data, you'd better provide initial values reasonably close to the solution. That's why used this initializer std::array{0.444, 0.11} rather than the original {0, 1}, as the latter does not converge to anything close to the correct answer.

Do we really need to use explicit function objects? Perhaps lambdas will do? Yes, they will - this compiles and runs as expected:

auto r3 = curve_fit([a](double x, double b, double c) { return gaussian(x, a, b, c); }, std::array{0.444, 0.11}, xs, ys);

Here's the full diff between the original and modified code (without lambda):

7a8
> #include <array>
72c73,74
< auto internal_make_gsl_vector_ptr(const std::vector<double>& vec) -> gsl_vector*
---
> template<auto n>
> auto internal_make_gsl_vector_ptr(const std::array<double, n>& vec) -> gsl_vector*
158,159c160,161
< template <typename Callable>
< auto curve_fit(Callable f, const std::vector<double>& initial_params, const std::vector<double>& x,
---
> template <typename Callable, auto n>
> auto curve_fit(Callable f, const std::array<double, n>& initial_params, const std::vector<double>& x,
163,164c165,166
<     constexpr auto n = 3;
<     assert(initial_params.size() == n);
---
>     //    constexpr auto n = 2;
>     //    assert(initial_params.size() == n);
194a197,204
> 
> struct gaussian_fixed_a
> {
>     double a;
>     gaussian_fixed_a(double a) : a{a} {}
>     double operator()(double x, double b, double c) const { return gaussian(x, a, b, c); }
> };
> 
212c222,224
<     auto r = curve_fit(gaussian, {1.0, 0.0, 1.0}, xs, ys);
---
>     auto r = curve_fit(gaussian, std::array{1.0, 0.0, 1.0}, xs, ys);
>     gaussian_fixed_a g(a);
>     auto r2 = curve_fit(g, std::array{0.444, 0.11}, xs, ys);
215a228,230
>     std::cout << "\n";
>     std::cout << "result: " << r2[0] << ' ' << r2[1] << '\n';
>     std::cout << "error : " << r2[0] - b << ' ' << r2[1] - c << '\n';