/* multimin/test_funcs.c
 * 
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Fabrice Rossi
 * 
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 3 of the License, or (at
 * your option) any later version.
 * 
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

#include <config.h>
#include <math.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_multimin.h>

#include "test_funcs.h"

gsl_multimin_function_fdf simpleabs =
{&simpleabs_f,
 &simpleabs_df,
 &simpleabs_fdf,
 2, 0};

gsl_multimin_function simpleabs_fmin =
{&simpleabs_f,
 2, 0};

void simpleabs_initpt (gsl_vector * x)
{
  gsl_vector_set (x, 0, 1.0);
  gsl_vector_set (x, 1, 2.0);
}

void simpleabs_initpt1 (gsl_vector * x)
{
  gsl_vector_set (x, 0, 1.0);
  gsl_vector_set (x, 1, 1.0);
}

double simpleabs_f (const gsl_vector * x, void *params)
{
  double u = gsl_vector_get(x,0);
  double v = gsl_vector_get(x,1);
  double a = u - 1;
  double b = v - 2;
  fcount++;
  return fabs(a) + fabs(b);
}

void simpleabs_df (const gsl_vector * x, void *params, gsl_vector * df)
{
  double u = gsl_vector_get(x,0);
  double v = gsl_vector_get(x,1);
  gcount++;
  gsl_vector_set(df,0, GSL_SIGN(u-1));
  gsl_vector_set(df,1, GSL_SIGN(v-2));  
}

void simpleabs_fdf (const gsl_vector * x, void *params, double * f,
                     gsl_vector * df) 
{
  double u = gsl_vector_get(x,0);
  double v = gsl_vector_get(x,1);
  double a = u - 1;
  double b = v - 2;
  gcount++;
  *f = fabs(a) + fabs(b);
  gsl_vector_set(df,0, GSL_SIGN(u-1));
  gsl_vector_set(df,1, GSL_SIGN(v-2));  
}

gsl_multimin_function_fdf rosenbrock =
{&rosenbrock_f,
 &rosenbrock_df,
 &rosenbrock_fdf,
 2, 0};

gsl_multimin_function rosenbrock_fmin =
{&rosenbrock_f,
 2, 0};

void rosenbrock_initpt (gsl_vector * x)
{
  gsl_vector_set (x, 0, -1.2);
  gsl_vector_set (x, 1, 1.0);
}

void rosenbrock_initpt1 (gsl_vector * x)
{
  gsl_vector_set (x, 0, 1.0);
  gsl_vector_set (x, 1, 1.0);
}

double rosenbrock_f (const gsl_vector * x, void *params)
{
  double u = gsl_vector_get(x,0);
  double v = gsl_vector_get(x,1);
  double a = u - 1;
  double b = u * u - v;
  fcount++;
  return a * a + 10 * b * b;
}

void rosenbrock_df (const gsl_vector * x, void *params, gsl_vector * df)
{
  double u = gsl_vector_get(x,0);
  double v = gsl_vector_get(x,1);
  double b = u * u - v;
  gcount++;
  gsl_vector_set(df,0,2 * (u - 1) + 40 * u * b);
  gsl_vector_set(df,1,-20 * b);  
}

void rosenbrock_fdf (const gsl_vector * x, void *params, double * f,
                     gsl_vector * df) 
{
  double u = gsl_vector_get(x,0);
  double v = gsl_vector_get(x,1);
  double a = u - 1;
  double b = u * u - v;
  gcount++;
  *f = a * a + 10 * b * b;
  gsl_vector_set(df,0,2 * (u - 1) + 40 * u * b);
  gsl_vector_set(df,1,-20 * b);  
}

gsl_multimin_function_fdf roth =
{&roth_f,
 &roth_df,
 &roth_fdf,
 2, 0};

gsl_multimin_function roth_fmin =
{&roth_f,
 2, 0};

void roth_initpt (gsl_vector * x)
{
  gsl_vector_set (x, 0, 4.5);
  gsl_vector_set (x, 1, 3.5);
}

double roth_f (const gsl_vector * x, void *params)
{
  double u = gsl_vector_get(x,0);
  double v = gsl_vector_get(x,1);
  double a = -13.0 + u + ((5.0 - v)*v - 2.0)*v;
  double b = -29.0 + u + ((v + 1.0)*v - 14.0)*v;
  fcount++;
  return a * a + b * b;
}

void roth_df (const gsl_vector * x, void *params, gsl_vector * df)
{
  double u = gsl_vector_get(x,0);
  double v = gsl_vector_get(x,1);
  double a = -13.0 + u + ((5.0 - v)*v - 2.0)*v;
  double b = -29.0 + u + ((v + 1.0)*v - 14.0)*v;
  double c = -2 + v * (10 - 3 * v);
  double d = -14 + v * (2 + 3 * v);
  gcount++;
  gsl_vector_set(df,0,2 * a + 2 * b);
  gsl_vector_set(df,1,2 * a * c + 2 * b * d);
}

void roth_fdf (const gsl_vector * x, void *params, double * f,
               gsl_vector * df) 
{
  *f = roth_f (x,params);
  roth_df(x,params,df);
}

gsl_multimin_function_fdf wood =
{&wood_f,
 &wood_df,
 &wood_fdf,
 4, 0};

gsl_multimin_function wood_fmin =
{&wood_f,
 4, 0};

void
wood_initpt (gsl_vector * x)
{
  gsl_vector_set (x, 0, -3.0);
  gsl_vector_set (x, 1, -1.0);
  gsl_vector_set (x, 2, -3.0);
  gsl_vector_set (x, 3, -1.0);
}

double wood_f (const gsl_vector * x, void *params)
{
  double u1 = gsl_vector_get(x,0);
  double u2 = gsl_vector_get(x,1);
  double u3 = gsl_vector_get(x,2);
  double u4 = gsl_vector_get(x,3);

  double t1 = u1 * u1 - u2;
  double t2 = u3 * u3 - u4;
  fcount++;
  return 100 * t1 * t1 + (1 - u1) * (1 - u1)
    + 90 * t2 * t2 + (1 - u3) * (1 - u3)
    + 10.1 * ( (1 - u2) * (1 - u2) + (1 - u4) * (1 - u4) )
    + 19.8 * (1 - u2) * (1 - u4);
}

void wood_df (const gsl_vector * x, void *params, gsl_vector * df)
{
  double u1 = gsl_vector_get(x,0);
  double u2 = gsl_vector_get(x,1);
  double u3 = gsl_vector_get(x,2);
  double u4 = gsl_vector_get(x,3);

  double t1 = u1 * u1 - u2;
  double t2 = u3 * u3 - u4;
  gcount++;
  gsl_vector_set(df,0, 400 * u1 * t1 - 2 * (1 - u1) );
  gsl_vector_set(df,1, -200 * t1 - 20.2 * (1 - u2) - 19.8 * (1 - u4) );
  gsl_vector_set(df,2, 360 * u3 * t2 - 2 * (1 - u3) );
  gsl_vector_set(df,3, -180 * t2 - 20.2 * (1 - u4) - 19.8 * (1 - u2) );
  
}

void wood_fdf (const gsl_vector * x, void *params, double * f, gsl_vector * df)
{
  wood_df(x,params,df);
  *f=wood_f(x,params);
}

gsl_multimin_function_fdf Nrosenbrock =
{&rosenbrock_f,
 &Nrosenbrock_df,
 &Nrosenbrock_fdf,
 2, 0};

void Nrosenbrock_df (const gsl_vector * x, void *params, gsl_vector * df)
{
  gsl_multimin_function F ;
  F.f = rosenbrock_f;
  F.params = params;
  F.n = x->size;
  gsl_multimin_diff (&F, x, df);
}

void Nrosenbrock_fdf (const gsl_vector * x, void *params, double * f,
                     gsl_vector * df) 
{
  *f = rosenbrock_f (x, params);
  Nrosenbrock_df (x, params, df);
}

gsl_multimin_function_fdf Nroth =
{&roth_f,
 &Nroth_df,
 &Nroth_fdf,
 2, 0};

void Nroth_df (const gsl_vector * x, void *params, gsl_vector * df)
{
  gsl_multimin_function F ;
  F.f = roth_f;
  F.params = params;
  F.n = x->size;
  gsl_multimin_diff (&F, x, df);
}

void Nroth_fdf (const gsl_vector * x, void *params, double * f,
                     gsl_vector * df) 
{
  *f = roth_f (x, params);
  Nroth_df (x, params, df);
}


gsl_multimin_function_fdf Nwood =
{&wood_f,
 &Nwood_df,
 &Nwood_fdf,
 4, 0};

void Nwood_df (const gsl_vector * x, void *params, gsl_vector * df)
{
  gsl_multimin_function F ;
  F.f = wood_f;
  F.params = params;
  F.n = x->size;
  gsl_multimin_diff (&F, x, df);
}

void Nwood_fdf (const gsl_vector * x, void *params, double * f,
                     gsl_vector * df) 
{
  *f = wood_f (x, params);
  Nwood_df (x, params, df);
}


gsl_multimin_function spring_fmin = { &spring_f,
  3, 0
};

void
spring_initpt (gsl_vector * x)
{
  gsl_vector_set (x, 0, 1.0);
  gsl_vector_set (x, 1, 0.0);
  gsl_vector_set (x, 2, 7.0 * M_PI);
}

double
spring_f (const gsl_vector * x, void *params)
{
  double x0 = gsl_vector_get (x, 0);
  double x1 = gsl_vector_get (x, 1);
  double x2 = gsl_vector_get (x, 2);

  double theta = atan2 (x1, x0);
  double r = sqrt (x0 * x0 + x1 * x1);
  double z = x2;
  while (z > M_PI)
    z -= 2.0 * M_PI;
  while (z < -M_PI)
    z += 2.0 * M_PI;
  {
    double tmz = theta - z;
    double rm1 = r - 1.0;
    double ret = 0.1 * (expm1 (tmz * tmz + rm1 * rm1) + fabs (x2 / 10.0));
    return ret;
  }
}