#include <stdio.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_min.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_multimin.h>

int globalcounter;
int globalcounter2;

// f(v0,v1,v2) = (v0-1)^2 + (v1-2)^2 + (v2-5)^2
double
my_f (const gsl_vector *v, void *params)
{
  double a = gsl_vector_get(v,0)-1.0;
  double b = gsl_vector_get(v,1)-2.0;
  double c = gsl_vector_get(v,2)-5.0;

  ++globalcounter;

  return a*a+b*b+c*c;
}

// Partial derivative in each component.  ( d/dv1(f), d/dv2(f), d/dv3(f) )
//  ie d/dv1 (f) = 2(v0-1)^1.1 = 2(v0-1)
void 
my_df (const gsl_vector *v, void *params, 
  gsl_vector *df)
{
  double a = gsl_vector_get(v,0);
  double b = gsl_vector_get(v,1);
  double c = gsl_vector_get(v,2);

  gsl_vector_set(df, 0, 2.0*(a-1.0));
  gsl_vector_set(df, 1, 2.0*(b-2.0));
  gsl_vector_set(df, 2, 2.0*(c-5.0));

  ++globalcounter2;
}


/* Compute both f and df together. */
void 
my_fdf (const gsl_vector *x, void *params,
        double *f, gsl_vector *df)
{
  *f = my_f(x,params);
  my_df(x,params,df);
}

void intro()
{
  printf("f(X) and d/dX f(X) are known functions of X.\n");
  printf("Using library minimizer to solve from a given point X0.\n");

}

int main(int argc, char** argv)
{
  intro();

  size_t iter = 0;
  int status;

  const gsl_multimin_fdfminimizer_type *T;
  gsl_multimin_fdfminimizer *s;

  //double p[3] = {0.0,0.0,0.0};

  gsl_vector *x;
  gsl_multimin_function_fdf my_func;

  my_func.f = &my_f;
  my_func.df = &my_df;
  my_func.fdf = &my_fdf;
  my_func.n = 3;
  //my_func.params = (void *)p;

  /* Starting point */
  x = gsl_vector_alloc (3);
  gsl_vector_set (x, 0, 0.0);
  gsl_vector_set (x, 1, 0.0);
  gsl_vector_set (x, 2, 0.0);

  T = gsl_multimin_fdfminimizer_conjugate_fr;
  s = gsl_multimin_fdfminimizer_alloc (T, 3);
 
  gsl_multimin_fdfminimizer_set (s, &my_func, x, 0.01, 1e-4);

  printf("iteration  x0 x1 x2 f(X)\n");

  globalcounter=0;
  globalcounter2=0;
  do
  {
    iter++;
    status = gsl_multimin_fdfminimizer_iterate (s);

    if (status)
      break;

    status = gsl_multimin_test_gradient (s->gradient, 1e-3);

    if (status==GSL_SUCCESS)
      printf("Minimum found at:\n");

    printf("%5d %.5f %.5f %.5f %10.5f\n", iter,
      gsl_vector_get (s->x, 0), 
      gsl_vector_get (s->x, 1), 
      gsl_vector_get (s->x, 2), 
      s->f);
  }
  while (status == GSL_CONTINUE && iter < 100);

  gsl_multimin_fdfminimizer_free (s);
  gsl_vector_free (x);

  printf("f(X) evaluations=%10d\n", globalcounter);
  printf("d/dX f(X) evaluations=%5d\n", globalcounter2);
  printf("function evaluations=%6d\n",globalcounter+globalcounter2);

  return 0;
}