c++ - GSL non-linear least squares fit won't converge -

i relatively new c++ , gsl. nonetheless have undertaken endeavour solve problem of fitting non-linear functions data implementations provided gsl. idea have class (clnm), holds data , model coefficients etc. clnm class has method cnlm::fitmodel_manu takes initial values , function pointers functions minimized, provided in fit_functions.cpp.

i have minimal working example first generate data using formula , coefficients in example gsl reference manual provides ( y = exp(lambda*x)+b ,
{a,lambda,b}=<5.0,-0.1,1.0>). works out fine , data expected.

with data cnlm object constructed. fitmodel_manu method called.

the implementation works , code compiles fine (using g++ , iso c++11, g++ -std=c++0x), cannot seem convergence ever, if provide correct parameters initial guess.

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 specified tolerance in f iteration: 3    c0 = -15    c1 = 2170.38    c2 = 16         |f(x)| =  901.11  

i have looked through quite posts, not find solution problem.

i @ end of wits. have provided code below if wants try out. maybe else can shed light?

thanks in advance!

---

this code

the header, 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 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 data     std::vector<double> in;     for(unsigned = 0; < 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 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 provide gsl_multifit_function_fdf via function     struct data d = { n, &m_data.y[0], &m_data.x[0], &m_data.sigma[0]}; // 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 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 starting point , 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 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 , 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 = 0; < 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 = 0; < 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 y = vmax * x/(km + x)\n#\n";     std::cout << "# \tvmax \tkm\n# value\t";     for(unsigned = 0; < 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 = 0; < 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 = 0; < s->x->size; i++){         std::cout << "\tc" << << " = " << 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 = gsl_vector_get (x, 0);     double lambda = gsl_vector_get (x, 1);     double b = gsl_vector_get (x, 2);     size_t i;     (i = 0; < n; i++)     {         /* model yi = * exp(-lambda * i) + b */         double t = i;         double yi = * 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 = gsl_vector_get (x, 0);     double lambda = gsl_vector_get (x, 1);     size_t i;     (i = 0; < n; i++)     {         /* jacobian matrix j(i,j) = dfi / dxj, */         /* fi = (yi - yi)/sigma[i],          */         /* yi = * exp(-lambda * i) + b */         /* , xj 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 * * 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 = 0; < out.p; i++){         out.param.push_back(param[i]);     }      // start rng     type = gsl_rng_default;     r = gsl_rng_alloc (type);      // simulate data noise     (unsigned = 0; < 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; } 

i found answer. seems had bug in code. defined struct data as,

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

, upon initialization switched x , y around:

struct data d = { n, &m_data.y[0], &m_data.x[0], &m_data.sigma[0]};

this lead fit trying fit x against y. guess noob mistake :) anyway