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);
}