Подгонка нелинейного метода наименьших квадратов GSL не будет сходиться

Я относительно новичок в C++ и особенно в GSL. Тем не менее я предпринял попытку решить проблему подгонки нелинейных функций к некоторым данным с помощью реализаций, предоставляемых GSL. Идея состоит в том, чтобы иметь класс (Clnm), который содержит данные, коэффициенты модели и т. д. Класс Clnm имеет метод Cnlm::fitModel_manu, который принимает начальные значения и указатели функций на минимизируемые функции, которые предоставляются в fit_functions.cpp.

У меня есть минимальный рабочий пример, в котором я сначала генерирую некоторые данные, используя формулу и коэффициенты, как в примере, приведенном в справочном руководстве GSL (Y = A exp(lambda*x)+b,
{ A, лямбда, b}=‹5,0,-0,1,1,0>). Это работает нормально, и я получаю данные, как и ожидалось.

С этими данными создается объект Cnlm. Затем вызывается метод fitModel_manu.

Реализация работает, и код компилируется нормально (с использованием g++ и ISO C++11, g++ -std=c++0x), однако я не могу добиться сходимости, даже если я укажу правильные параметры в качестве первоначального предположения.

iteration: 0    c0 = 1  c1 = 5  c2 = -0.1
        |f(x)| =  1954.18
        status = success
iteration: 1    c0 = -15.513    c1 = 2170.38    c2 = 16.5136
        |f(x)| =  902.806
        status = success
iteration: 2    c0 = -15    c1 = 2170.38    c2 = 16
        |f(x)| =  901.11
status = cannot reach the specified tolerance in F
iteration: 3    c0 = -15    c1 = 2170.38    c2 = 16
        |f(x)| =  901.11 

Я просмотрел довольно много сообщений, но не смог найти решение своей проблемы.

Я в конце моего ума. Я предоставил свой код ниже, если кто-то хочет попробовать его. Может быть, кто-то еще может пролить свет?

Заранее спасибо!

это код

Заголовок, nlls.h:

#ifndef NLLS_H_
#define NLLS_H_

#include<vector>
#include<string>
#include<gsl/gsl_blas.h>
#include<gsl/gsl_randist.h>
#include<gsl/gsl_rng.h>
#include<gsl/gsl_vector.h>
#include<gsl/gsl_multifit_nlin.h>

struct data {
    size_t n;
    double * x;
    double * y;
    double * sigma;
};

struct simulation {
    std::vector<double> x, y, weights, sigma, param;
    size_t p;
    size_t n;
};

// Class definition
class Cnlm {
public:
    Cnlm(simulation xy_data);
    const void printSummary();
    int fitModel_manu(std::vector<double> x_init,
            int (*f)(const gsl_vector *, void *, gsl_vector *),
            int (*df)(const gsl_vector *, void *, gsl_matrix *));
private:
    double m_sumsq, m_DOF, m_c;
    std::vector<double> m_param, m_ERR;
    struct resnorm {
        double chi, chi0;
    } m_resnorm;
    std::string m_solvername;
    unsigned niter, n_fun_eval, n_jacobian_eval;
    int m_info, m_status;
    simulation m_data;
};

void print_state(size_t iter, gsl_multifit_fdfsolver * s);

// function to simulate data
simulation expb_simulate(std::vector<double> x, std::vector<double> param,
        double si);

// fit functions
int expb_f(const gsl_vector * x, void *data, gsl_vector * f);
int expb_df(const gsl_vector * x, void *data, gsl_matrix * J);

#endif /* NLLS_H_ */

main.cpp:

#include<iostream>
#include"nlls.h"

int main() {
    // simulate some data
    std::vector<double> in;
    for(unsigned i = 0; i < 30; i++){
        in.push_back(i+1);
    }
    std::vector<double> exp_par = {1, 5, -0.1};
    simulation exp_sim = expb_simulate(in, exp_par, 0.05);

    // create Clnm object
    Cnlm mynLM(exp_sim);
    // fit with initial parameters 
    std::vector<double> initial_x = {1, 5, -0.1};
    mynLM.fitModel_manu(initial_x, &expb_f, &expb_df);
    mynLM.printSummary();

    return 0;
}

class_functions.cpp:

#include"nlls.h"
#include<iostream>

// constructor
Cnlm::Cnlm(simulation xy_data){
    m_data = xy_data;
}

int Cnlm::fitModel_manu(std::vector<double> x_init,  int (*f)(const gsl_vector * , void *, gsl_vector *), int (*df)(const gsl_vector * , void *, gsl_matrix * )){

    const size_t n = m_data.x.size();
    const size_t p = m_data.p;
    m_info = 0;
    unsigned iter = 0;

    //construct data struct to provide to gsl_multifit_function_fdf via function
    struct data d = { n, &m_data.y[0], &m_data.x[0], &m_data.sigma[0]}; // from simulated data
    gsl_vector_view w = gsl_vector_view_array(&m_data.weights[0], n);

    const gsl_multifit_fdfsolver_type * T = gsl_multifit_fdfsolver_lmsder;
    gsl_multifit_fdfsolver *s;

    gsl_matrix *J = gsl_matrix_alloc(n, p);
    gsl_matrix *covar = gsl_matrix_alloc (p, p);

    gsl_multifit_function_fdf FUN;
    gsl_vector *res_f;

    FUN.f = f;
    FUN.df = df;
    //FUN.df = NULL;
    FUN.n = n;
    FUN.p = p;
    FUN.params = &d;

    // initial values for determination of parameters
    gsl_vector_view x = gsl_vector_view_array (&x_init[0], p);

    // create new solver
    s = gsl_multifit_fdfsolver_alloc (T, n, p);

    /* initialize solver with starting point and weights */
    gsl_multifit_fdfsolver_wset (s, &FUN, &x.vector, &w.vector);

    /* compute initial residual norm */
    res_f = gsl_multifit_fdfsolver_residual(s);
    m_resnorm.chi0 = gsl_blas_dnrm2(res_f);

    print_state(iter, s);

    // solve by iteration
    do{
        iter++;
        m_status = gsl_multifit_fdfsolver_iterate (s);
        std::cout << "status = " << gsl_strerror (m_status) << std::endl;

        print_state (iter, s);

        if(m_status) break;

        m_status = gsl_multifit_test_delta (s->dx, s->x, 1e-4, 1e-4);

    }while (m_status == GSL_CONTINUE && iter < 500);

    // compute covariance matrix and best fit parameters
    gsl_multifit_fdfsolver_jac(s, J);
    gsl_multifit_covar (J, 0.0, covar);

    /* compute final residual norm */
    m_resnorm.chi = gsl_blas_dnrm2(res_f);

    m_solvername = gsl_multifit_fdfsolver_name(s);
    m_info = 0;
    niter = gsl_multifit_fdfsolver_niter(s);
    n_fun_eval = FUN.nevalf;
    n_jacobian_eval = FUN.nevaldf;
    m_DOF = n-p; //compute degrees of freedom
    m_c = GSL_MAX_DBL(1, m_resnorm.chi / sqrt(m_DOF));
    for(unsigned i = 0; i < p; i++){
        m_param.push_back(gsl_vector_get(s->x, i));
        m_ERR.push_back(gsl_matrix_get(covar,i,i));


    }
    gsl_multifit_fdfsolver_free (s);
    gsl_matrix_free (covar);
    gsl_matrix_free (J);

    return 0;
}

const void Cnlm::printSummary(){
    std::cout << "# Data\n#\tx\ty\tweight\tsigma\n";
    for(unsigned i = 0; i < m_data.x.size(); i++){
        std::cout << "\t" << m_data.x[i] << "\t" << m_data.y[i]  << "\t" << m_data.weights[i]  << "\t" << m_data.sigma[i] << std::endl;
    }
    std::cout << "#\n# Summary\n# -------------------------------------------\n";
    std::cout << "# least squares fit estimates for Y = Vmax * X/(Km + X)\n#\n";
    std::cout << "# \tVmax \tKm\n# value\t";
    for(unsigned i = 0; i < m_param.size(); i++){
        if(i == 0)  std::cout << "\t" << m_param[i];
        else std::cout << "\t" << m_param[i];
    }
    std::cout << std::endl;
    std::cout << "# err\t";
    for(unsigned i = 0; i < m_ERR.size(); i++){
        if(i == 0)  std::cout << "\t" << m_ERR[i];
        else std::cout << "\t" << m_ERR[i];
    }
    std::cout << std::endl;
    std::cout << "#\n# Chisq/DOF: " << pow(m_resnorm.chi, 2.0)/m_DOF << std::endl;
}


// General functions
void print_state (size_t iter, gsl_multifit_fdfsolver * s){
    std::cout << "iteration: " << iter;
    for(unsigned i = 0; i < s->x->size; i++){
        std::cout << "\tc" << i << " = " << gsl_vector_get (s->x, i);
    }
    std::cout << "\n\t\t|f(x)| =  " << gsl_blas_dnrm2 (s->f) << std::endl;
}

fit_functions.cpp:

#include"nlls.h"


int expb_f (const gsl_vector * x, void *data, gsl_vector * f)
{
    size_t n = ((struct data *)data)->n;
    double *y = ((struct data *)data)->y;
    double A = gsl_vector_get (x, 0);
    double lambda = gsl_vector_get (x, 1);
    double b = gsl_vector_get (x, 2);
    size_t i;
    for (i = 0; i < n; i++)
    {
        /* Model Yi = A * exp(-lambda * i) + b */
        double t = i;
        double Yi = A * exp (-lambda * t) + b;
        gsl_vector_set (f, i, Yi - y[i]);
    }
    return GSL_SUCCESS;
}

int expb_df (const gsl_vector * x, void *data,  gsl_matrix * J)
{
    size_t n = ((struct data *)data)->n;
    double A = gsl_vector_get (x, 0);
    double lambda = gsl_vector_get (x, 1);
    size_t i;
    for (i = 0; i < n; i++)
    {
        /* Jacobian matrix J(i,j) = dfi / dxj, */
        /* where fi = (Yi - yi)/sigma[i],
         */
        /*
Yi = A * exp(-lambda * i) + b */
        /* and the xj are the parameters (A,lambda,b) */
        double t = i;
        double e = exp(-lambda * t);
        gsl_matrix_set (J, i, 0, e);
        gsl_matrix_set (J, i, 1, -t * A * e);
        gsl_matrix_set (J, i, 2, 1.0);
    }
    return GSL_SUCCESS;
}

simulation expb_simulate(std::vector<double> x, std::vector<double> param, double si){

    gsl_rng * r;
    const gsl_rng_type * type;
    gsl_rng_env_setup();
    unsigned n = x.size();

    simulation out;
    out.p = param.size();
    out.n = n;
    out.x.assign(x.begin(), x.end());
    out.sigma.assign(n, si);

    for(unsigned i = 0; i < out.p; i++){
        out.param.push_back(param[i]);
    }

    // start RNG
    type = gsl_rng_default;
    r = gsl_rng_alloc (type);

    // simulate data with some noise
    for (unsigned i = 0; i < n; i++)
    {
        double yi = out.param[0]+ out.param[1] * exp (out.param[2] * x[i]);
        double s = si;
        double dy = gsl_ran_gaussian(r, s);

        out.y.push_back(yi + dy);
        out.weights.push_back(1.0/(s*s));
    }

    gsl_rng_free (r);
    return out;
}