Re: [Help-gsl] Nonlinear fitting

[Ken O'Brien]

Hi,

Here are the relevant functions in the attachment.

Thanks in advance

Ken


> On Fri, 2014-07-04 at 08:11 +0200, Marco Maggi wrote: 
> 
> > k...@kenobrien.org wrote:
> > 
> > > Hi,
> > > I am trying to fit the nonlinear equation:
> > > P = C0 + C1*X + C2*X^r
> > > where C0, C1, C2 and r are variables to be fitted, P and X are vectors
> > > of known observations.
> > > I have read the docs and I'm trying to modify the example program. I'm
> > > struggling to provide the expb_f, expb_df, expb_fdf functions.
> > > I would like advice on how to proceed.
> > 
> > Can  you  post a  short  version  of  the program  you  have
> > composed so far?
> 
> 
> 
> 

-- 
PhD student,
Simulation Science and Extreme Events Cluster / Heterogeneous Computing
Laboratory
University College Dublin,
Ireland

struct data {
        size_t n;
        double * y;
        double * cpu;
        double * power;
};


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 *cpu = ((struct data *) data)->cpu;
        double *power = ((struct data *)data)->power;

        double C0 = gsl_vector_get (x, 0);
        double C1 = gsl_vector_get (x, 1);
        double C2 = gsl_vector_get (x, 2);
        double r  = gsl_vector_get (x, 3);
        size_t i;

        for (i = 0; i < n; i++) {
      /* Model Yi = C0 * C1*cpu + C2*cpu^r */
      double t = i;
      double Yi = C0 + C1 * cpu[i] + C2*pow(cpu[i], r);
      gsl_vector_set (f, i, (Yi - y[i]));
    }

        return GSL_SUCCESS;
}

int expb_df(const gsl_vector *x, void *data, gsl_matrix *J) {
        return GSL_SUCCESS;
}

int expb_fdf(const gsl_vector *x, void *data, gsl_vector *f, gsl_matrix *J) {
        return GSL_SUCCESS;
}


/*
        Implement P=Co+C1*Ucpu+C2*ucpu^r
        @param observations Number of observations made.
        @param cpu_util CPU utilization as percentage of total system.
        @param power_util Watt readings during this CPU utilization.
*/

void
print_state (size_t iter, gsl_multifit_fdfsolver * s)
{
  printf ("iter: %3u x = % 15.8f % 15.8f % 15.8f "
          "|f(x)| = %g\n",
          (unsigned int) iter,
          gsl_vector_get (s->x, 0), 
          gsl_vector_get (s->x, 1),
          gsl_vector_get (s->x, 2), 
          gsl_blas_dnrm2 (s->f));
}

void
fitting_cpu_model_1(int observations, double *cpu_util, double *power_util,
                    char *dir)
{
        const gsl_multifit_fdfsolver_type *T;
        gsl_multifit_fdfsolver *s;
        const int p = 4;
        int status;
        gsl_matrix *covar = gsl_matrix_alloc (p, p);    
        double y[observations], sigma[observations];
        struct data d = { observations, y, cpu_util, power_util};
        unsigned int iter = 0;
        gsl_multifit_function_fdf f;
        double x_init[4] = { 1.0, 0.0, 0.0, 0.0 };
        gsl_vector_view x = gsl_vector_view_array (x_init, p);  

        f.f = &expb_f;
        f.df = &expb_df;
        f.fdf = &expb_fdf;
        f.n = observations;
        f.p = p;
        f.params = &d;
        
        T = gsl_multifit_fdfsolver_lmsder;
        s = gsl_multifit_fdfsolver_alloc(T, observations, p);
        
        gsl_multifit_fdfsolver_set(s, &f, &x.vector);


        do {    
                iter++;
                status = gsl_multifit_fdfsolver_iterate(s);
                printf ("status = %s\n", gsl_strerror (status));
                print_state(iter, s);
                if(status)
                        break;
                status = gsl_multifit_test_delta(s->dx, s->x, 1e-4, 1e-4);
        }
        while(status == GSL_CONTINUE && iter < 500);

        gsl_multifit_covar(s->J, 0.0, covar);
        
        #define FIT(i) gsl_vector_get(s->x, i)
        #define ERR(i) sqrt(gsl_matrix_get(covar,i,i))

        double chi = gsl_blas_dnrm2(s->f);
    double dof = observations - p;
    double c = GSL_MAX_DBL(1, chi / sqrt(dof)); 

    printf("chisq/dof = %g\n",  pow(chi, 2.0) / dof);

    printf ("C0      = %.5f +/- %.5f\n", FIT(0), c*ERR(0));
    printf ("C1 = %.5f +/- %.5f\n", FIT(1), c*ERR(1));
    printf ("C2 = %.5f +/- %.5f\n", FIT(2), c*ERR(2));
    printf ("r     = %.5f +/- %.5f\n", FIT(3), c*ERR(3));
        
        gsl_multifit_fdfsolver_free (s);
        gsl_matrix_free (covar);
}

signature.asc
Description: This is a digitally signed message part