/************************************************************************/
/*
opt-amoeba: optimize policy using Numerical Recipes Nelder Mead (Amoeba) algorithm
*/
/************************************************************************/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "nr.h"
#include "nrutil.h"

#include "main.h"
#include "main2.h"

/*****************************************************************************/
/*****************************************************************************/
/*****************************************************************************/

extern SIM sim;

/*****************************************************************************/
/* Numerical Recipes Amoeba stuff */

static float **p; /* ndim+1 x ndim dimensional points in the simplex */
static float p_offsets[MAX_N_PARAMETERS+1];
static float y[MAX_N_PARAMETERS+2]; /* values for points in p */
static float init_parameters[MAX_N_PARAMETERS+1];
static int iter = 0;
#define FTOL 1.0e-6

/*****************************************************************************/
/*****************************************************************************/
/* Main: initialize and then call amoeba optimizer */

int main( int argc, char **argv )
{
  int i, j, k;
  float **matrix();
  double step = 0.01;
  double cost2;
  PARAMETER *read_parameter_file();
  
  init_default_parameters( &sim );

  if ( argc < 2 )
    {
      fprintf( stderr, "Optimize using which parameter file?\n" );
      fprintf( stderr, "%s parameter-file\n", argv[0] );
      exit( -1 );
    }
  sim.params = read_parameter_file( argv[1] );
  sim.n_parameters = process_parameters( sim.params, &sim, 1 );
  if ( sim.n_parameters > MAX_N_PARAMETERS )
    {
      fprintf( stderr, "Too many parameters %d > %d\n",
	       sim.n_parameters, MAX_N_PARAMETERS );
      exit( -1 );
    }
  sprintf( sim.output_file, "%s.new", argv[1] );

  init_sim( &sim );

  sim.n_func_calls_per_eval = 1;
  sim.all_time_low_cost = 1e20;
  sim.debug_criterion = 1;

  parameters_to_vector( sim.params, init_parameters );
  for( i = 0; i < sim.n_parameters; i++ )
    {
      p_offsets[i] = 0.01*init_parameters[i];
      if ( p_offsets[i] == 0.0 )
	p_offsets[i] = 0.01;
    }
  p = matrix( 1, sim.n_parameters+1, 1, sim.n_parameters );
  for( i = 1; i <= sim.n_parameters+1; i++ )
    {
      for ( j = 1; j <= sim.n_parameters; j++ )
	p[i][j] = init_parameters[j-1];
      if ( i >= 2 )
	p[i][i-1] += p_offsets[i-2];
      y[i] = call_many_func( p[i] );
      if ( i >= 2 && sim.status == CRASHED )
	{
	  printf( "Trying other direction: %d.\n", i );
	  p[i][i-1] -= 2*p_offsets[i-2];
	  y[i] = call_many_func( p[i] );
	}
    }

  amoeba( p, y, sim.n_parameters, FTOL, call_many_func, &(sim.iter) );

  cost2 = call_many_func( p[1] );

  printf( "Iterations: %d (%d)\n\n", sim.iter, sim.func_calls );
  printf( "Minimum found at: \n" );
  for ( i = 1; i <= sim.n_parameters; i++ )
    printf( "p[%d] = %18.12f;\n", i-1, p[1][i] );
  printf( "\n\nMinimum function value = %12.6f, should equal %12.6f \n\n",
	  y[1], cost2 );
  printf( "All time low cost: %lg\n", sim.all_time_low_cost );
}

/************************************************************************/
