help-gsl
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: how to debug multifit_nlinear


From: Patrick Alken
Subject: Re: how to debug multifit_nlinear
Date: Wed, 19 Jun 2024 10:03:41 -0400
User-agent: Mozilla Thunderbird

Hello, it looks like a problem with the finite-difference approximation to the Jacobian matrix. Can you provide a MWE we can test?

On 6/19/24 03:23, PICCA Frederic-Emmanuel wrote:
[External email - use caution]


Hello,

I am using

#include <gsl/gsl_multifit_nlinear.h>

in order to do a least square fitting of this function

int model_f (const gsl_vector *x, void *data, gsl_vector *f)
{
         size_t n = ((struct data *)data)->n;
         float *t = ((struct data *)data)->t;
         float *y = ((struct data *)data)->y;
         int* lista = ((struct data *)data)->lista;
         size_t n_lista = ((struct data *)data)->n_lista;
         float* parameters = ((struct data *)data)->parameters;
         size_t n_parameters = ((struct data *)data)->n_parameters;
         FitFunction funcs = ((struct data *)data)->funcs;

         float dyda;

         size_t i;

         for (i=0; i<n_lista; ++i){
                 parameters[lista[i] - 1] = gsl_vector_get(x, i);
         }

         for(i=0; i<n; ++i){
                 float ymod;

                 (*funcs)(t[i], parameters, &ymod, &dyda, n_parameters, 1);
                 ymod = t[i];
                 fprintf(stdout, "f(i=%d, t=%f) = %f (%f)\n", i, t[i], ymod, 
y[i]);
                 gsl_vector_set (f, i, ymod - y[i]);
         }
         fprintf(stdout, "\n");

         return GSL_SUCCESS;
}

with

         const gsl_multifit_nlinear_type *T = gsl_multifit_nlinear_trust;
         gsl_multifit_nlinear_parameters fdf_params = 
gsl_multifit_nlinear_default_parameters();
         struct data d = {
                 .n = n,
                 .t = xx,
                 .y = yy,
                 .wgt = wgt,
                 .parameters = a,
                 .n_parameters = ma,
                 .lista = lista,
                 .n_lista = p,
                 .funcs = funcs,
         };

         gsl_rng_env_setup();
         r = gsl_rng_alloc(gsl_rng_default);

         /* define the function to be minimized */
         fdf.f = model_f;
         fdf.df = NULL;   /* set to NULL for finite-difference Jacobian */
         fdf.fvv = NULL;     /* not using geodesic acceleration */
         fdf.n = n;
         fdf.p = p;
         fdf.params = &d;

         /* allocate workspace with default parameters */
         w = gsl_multifit_nlinear_alloc (T, &fdf_params, n, p);

         /* initialize solver with starting point and weights */
         gsl_multifit_nlinear_winit (&x.vector, &wts.vector, &fdf, w);

         /* compute initial cost function */
         /* f = gsl_multifit_nlinear_residual(w); */
         /* gsl_blas_ddot(f, f, &chisq0); */

         /* solve the system with a maximum of 100 iterations */
         status = gsl_multifit_nlinear_driver(10000, xtol, gtol, ftol,
                                              callback, NULL, &info, w);


but when I fit the function (There is only 2 parameters and around 133
data points points), I have this error

iter  0: param0 = 0.2731, param1 = 0.000, cond(J) =   4.0000, |f(x)| = 566.7044
gsl: qrpt.c:345: ERROR: rank must have 0 < rank <= N
Default GSL error handler invoked.

Is it possible to explain the possible cause of this error. And  is
there a degug mode in order to observe the  internal of
gsl_multifit_nlinéar_driver during the fit.

thanks for your help

Frederic




reply via email to

[Prev in Thread] Current Thread [Next in Thread]