I am relatively new to C++ and especially the GSL. Nonetheless I have undertaken the endeavour to solve the problem of fitting non-linear functions to some data with the implementations provided by GSL. The idea is to have a Class (Clnm), which holds data and model coefficients etc. The Clnm class has a method Cnlm::fitModel_manu which takes initial values and the function pointers to the functions to be minimized, which are provided in the fit_functions.cpp.
I have a minimal working example where I first generate some data using the formula and coefficients as in the example the GSL reference manual provides ( Y = A exp(lambda*x)+b ,
{A,lambda,b}=<5.0,-0.1,1.0>). This works out fine and I get data as expected.
With this data a Cnlm object is constructed. Then the fitModel_manu method is called.
The implementation works and the code compiles fine (using g++ and ISO C++11, g++ -std=c++0x), however I cannot seem to get a convergence ever, even if I provide the correct parameters as 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 the specified tolerance in F
iteration: 3 c0 = -15 c1 = 2170.38 c2 = 16
|f(x)| = 901.11
I have looked through quite some posts, but I could not find a solution for my problem.
I am at the end of my wits. I have provided my code below if anyone wants to try it out. Maybe someone else can shed some light?
Thanks in advance!
This is the code
The header, nlls.h:
/*
* 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:
/*
* 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:
/*
* 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:
/*
* 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;
}
Aucun commentaire:
Enregistrer un commentaire