#ifndef _GLOBALGSLFUNCS_H_
#define _GLOBALGSLFUNCS_H_

residFunc *residFuncPtr; // global pointer to a fitter class for gsl callback

int fGsl(const gsl_vector *x, void *data, gsl_vector *f) {
  unsigned int n = residFuncPtr->getn(), p = residFuncPtr->getp(), i;
  vector<double> xVector(p, 0.0), fVector(n, 0.0);
  for (i = 0; i < p; i++) xVector[i] = gsl_vector_get(x, (size_t)i);
  residFuncPtr->resid(&xVector, data, &fVector);
  for (i = 0; i < n; i++) gsl_vector_set(f, (size_t)i, fVector[i]);
  return GSL_SUCCESS;
}

int dfGsl(const gsl_vector *x, void *data, gsl_matrix *J){
  int n = residFuncPtr->getn(), p = residFuncPtr->getp();
  int i, j;
  vector<double> xVector(p);
  vector< vector<double> > JVector(n,p);
  for (i = 0; i < p; i++) xVector[i] = gsl_vector_get(x, (size_t)i);
  residFuncPtr->deriv(&xVector, data, &JVector);
  for (i = 0; i < n; i++) {
    for (j = 0; j < p; j++) {
      gsl_matrix_set(J, (size_t)i, (size_t)j, JVector[i][j]);
    }
  }
  return GSL_SUCCESS;
}

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

#endif // _GLOBALGSLFUNCS_H_