/* 
   migrate.cc
   This file contains methods which will be called by FORTRAN
   subroutines. The purpose of these functions is to migrate
   information provided by FORTRAN users into the Direct
   Search libraries.
   
   Amy Yates, 2001
   The College of William and Mary, Williamsburg, Virginia,
   under advisor Dr. Virginia Torczon   
   
*/

#include <cstdio>           // for printf
#include <iostream>        // for cout
#include <string>          // for strcpy
#include "objective.h"    
#include "cppmat.h"
#include "CompassSearch.h" 
#include "CoordinateSearch.h" 
#include "EdHJSearch.h" 
#include "HJSearch.h" 
#include "NLessSearch.h" 
#include "NMSearch.h" 
#include "SHHSearch.h" 
#include "SMDSearch.h" 
#include "vec.h"
#include <assert.h>

/* The underscore after the names of methods is necessary
   because the FORTRAN appends it to all methods referenced */

// Declare all methods as extern "C" to avoid name-mangling
extern "C" void coord_(int *, double *, double *, double *, int *,
                       int *, double *, int * ,
                       void (*pt2Func)(double *, double *, long *));
extern "C" void compss_(int *, double *, double *, double *, int *,
                        int *, double *, int * ,
                        void (*pt2Func)(double *, double *, long *));
extern "C" void nless_(int *, double *, double *, double *, int *,
                       int *, double *, int * ,
                       void (*pt2Func)(double *, double *, long *));
extern "C" void hj_(int *, double *, double *, double *, int *,
                    int *, double *, int * ,
                    void (*pt2Func)(double *, double *, long *));
extern "C" void edhj_(int *, double *, double *, double *, int *,
                      int *, double *, int * ,
                      void (*pt2Func)(double *, double *, long *));
extern "C" void shhsch_(int *, double *, int *, double *, double *,
                        int *, int *, int *, double *, int * ,
                        void (*pt2Func)(double *, double *, long *));
extern "C" void nmsch_(int *, double *, int *, double *, double *,
                       double *, double *, double *, double *,
                       int *, int *, int *, double *, int * , 
                        void (*pt2Func)(double *, double *, long *));
extern "C" void smdsch_(int *, double *, int *, double *, double *,
                        int *, int *, int *, double *, int * , 
                        void (*pt2Func)(double *, double *, long *));

// Declare a global function pointer
void (*FuncPtr) (double *, double *, long *);

double CallFortranFunction(long n, Vector<double> &x)
{
  double arrayx[n];
  double f=0;
  for (long i=0; i < n; i++)
    {
      arrayx[i] = x[i];
    }
  FuncPtr(arrayx, &f, &n);
  return f;
}

void f_fcn(long vars, Vector<double> &x, double & f, bool & flag, void* nothing)
{
     
  if (vars < 1) {
      cerr << "\nError: Dimension cannot be less than 1!!\n";
      exit(1);
  }
  f = CallFortranFunction(vars, x);
  flag = true;
  return;
}

/* This method sets up a Coordinate search given parameters
   passed in through FORTRAN or C subroutines*/
void coord_(int *dim, double *start, double *step1, double *step2,
             int *max, int *istat, double *aux, int *iaux, 
	     void (*pt2Func)(double *, double *, long *) )
{
  double SMinVal;         // to hold the min value
  long Scalls;            // to hold the number of calls
  
  Vector<double> Sminimum(*dim);   // to store the minimum point later
  
  // initialize an n-entry Vector with values start
  Vector<double>minVec(*dim, start);

  //Check for defaults
  double dstep1, dstep2;
  if (*step1 == -1.0)
    dstep1 = 0.25;
  else
    dstep1 = *step1;
  if (*step2 == -1.0)
    dstep2 = 10E-8;
  else
    dstep2 = *step2;

  // Set the global function pointer
  FuncPtr = pt2Func;
  // Construct a Search object.
  CoordinateSearch Search(*dim, minVec, dstep1, dstep2, f_fcn, NULL);
  
  // Set the maximum number of calls
  Search.SetMaxCallsExact(*max);
    
  // BeginSearch() may go on indefinitely for a poorly defined problem
  Search.BeginSearch();
  
  Search.GetMinPoint(Sminimum);
  Search.GetMinVal(SMinVal);
  Scalls = Search.GetFunctionCalls();

  // Pass back the final step length, minimum function evaluation,
  // and the corresponding x
  aux[0] = Search.GetDelta()*2.0;
  aux[1] = SMinVal;
  for (int i=0; i<*dim;i++)
    aux[i+2] = Sminimum[i];

  // Pass back the number of function evaluations it took
  *iaux = Scalls;
  // If BeginSearch() stopped, assume success
  *istat = 1;
  
  return;
}

/* This method sets up a Compass search given parameters
   passed in through FORTRAN or C subroutines*/
void compss_(int *dim, double *start, double *step1, double *step2,
             int *max, int *istat, double *aux, int *iaux, 
	     void (*pt2Func)(double *, double *, long *) )
{
  double SMinVal;         // to hold the min value
  long Scalls;            // to hold the number of calls
  
  Vector<double> Sminimum(*dim);   // to store the minimum point later
  
  // initialize an n-entry Vector with values start
  Vector<double>minVec(*dim, start);

  //Check for defaults
  double dstep1, dstep2;
  if (*step1 == -1.0)
    dstep1 = 0.25;
  else
    dstep1 = *step1;
  if (*step2 == -1.0)
    dstep2 = 10E-8;
  else
    dstep2 = *step2;

  // Set the global function pointer
  FuncPtr = pt2Func;
  // Construct a Search object.
  CompassSearch Search(*dim, minVec, dstep1, dstep2, f_fcn, NULL);
  
  // Set the maximum number of calls
  Search.SetMaxCallsExact(*max);
    
  // BeginSearch() may go on indefinitely for a poorly defined problem
  Search.BeginSearch();
  
  Search.GetMinPoint(Sminimum);
  Search.GetMinVal(SMinVal);
  Scalls = Search.GetFunctionCalls();

  // Pass back the final step length, minimum function evaluation,
  // and the corresponding x
  aux[0] = Search.GetDelta()*2.0;
  aux[1] = SMinVal;
  for (int i=0; i<*dim;i++)
    aux[i+2] = Sminimum[i];

  // Pass back the number of function evaluations it took
  *iaux = Scalls;
  // If BeginSearch() stopped, assume success
  *istat = 1;
  
  return;
}

/* This method sets up an N-Less search given parameters
   passed in through FORTRAN or C subroutines*/
void nless_(int *dim, double *start, double *step1, double *step2,
             int *max, int *istat, double *aux, int *iaux, 
	     void (*pt2Func)(double *, double *, long *) )
{
  double SMinVal;         // to hold the min value
  long Scalls;            // to hold the number of calls
  
  Vector<double> Sminimum(*dim);   // to store the minimum point later
  
  // initialize an n-entry Vector with values start
  Vector<double>minVec(*dim, start);

  //Check for defaults
  double dstep1, dstep2;
  if (*step1 == -1.0)
    dstep1 = 0.25;
  else
    dstep1 = *step1;
  if (*step2 == -1.0)
    dstep2 = 10E-8;
  else
    dstep2 = *step2;

  // Set the global function pointer
  FuncPtr = pt2Func;
  // Construct a Search object.
  NLessSearch Search(*dim, minVec, dstep1, dstep2, f_fcn, NULL);
  
  // Set the maximum number of calls
  Search.SetMaxCallsExact(*max);
    
  // BeginSearch() may go on indefinitely for a poorly defined problem
  Search.BeginSearch();
  
  Search.GetMinPoint(Sminimum);
  Search.GetMinVal(SMinVal);
  Scalls = Search.GetFunctionCalls();

  // Pass back the final step length, minimum function evaluation,
  // and the corresponding x
  aux[0] = Search.GetDelta()*2.0;
  aux[1] = SMinVal;
  for (int i=0; i<*dim;i++)
    aux[i+2] = Sminimum[i];

  // Pass back the number of function evaluations it took
  *iaux = Scalls;
  // If BeginSearch() stopped, assume success
  *istat = 1;
  
  return;
}

/* This method sets up a Hooke and Jeeves search given parameters
   passed in through FORTRAN or C subroutines*/
void hj_(int *dim, double *start, double *step1, double *step2,
             int *max, int *istat, double *aux, int *iaux, 
	     void (*pt2Func)(double *, double *, long *) )
{
  double SMinVal;         // to hold the min value
  long Scalls;            // to hold the number of calls
  
  Vector<double> Sminimum(*dim);   // to store the minimum point later
  
  // initialize an n-entry Vector with values start
  Vector<double>minVec(*dim, start);

  //Check for defaults
  double dstep1, dstep2;
  if (*step1 == -1.0)
    dstep1 = 0.25;
  else
    dstep1 = *step1;
  if (*step2 == -1.0)
    dstep2 = 10E-8;
  else
    dstep2 = *step2;

  // Set the global function pointer
  FuncPtr = pt2Func;
  // Construct a Search object.
  HJSearch Search(*dim, minVec, dstep1, dstep2, f_fcn, NULL);
  
  // Set the maximum number of calls
  Search.SetMaxCallsExact(*max);
    
  // BeginSearch() may go on indefinitely for a poorly defined problem
  Search.BeginSearch();
  
  Search.GetMinPoint(Sminimum);
  Search.GetMinVal(SMinVal);
  Scalls = Search.GetFunctionCalls();

  // Pass back the final step length, minimum function evaluation,
  // and the corresponding x
  aux[0] = Search.GetDelta()*2.0;
  aux[1] = SMinVal;
  for (int i=0; i<*dim;i++)
    aux[i+2] = Sminimum[i];

  // Pass back the number of function evaluations it took
  *iaux = Scalls;
  // If BeginSearch() stopped, assume success
  *istat = 1;
  
  return;
}

/* This method sets up an edited Hooke and Jeeves search given 
   parameters passed in through FORTRAN or C subroutines*/
void edhj_(int *dim, double *start, double *step1, double *step2,
           int *max, int *istat, double *aux, int *iaux, 
           void (*pt2Func)(double *, double *, long *) )
{
  double SMinVal;         // to hold the min value
  long Scalls;            // to hold the number of calls
  
  Vector<double> Sminimum(*dim);   // to store the minimum point later
  
  // initialize an n-entry Vector with values start
  Vector<double>minVec(*dim, start);

  //Check for defaults
  double dstep1, dstep2;
  if (*step1 == -1.0)
    dstep1 = 0.25;
  else
    dstep1 = *step1;
  if (*step2 == -1.0)
    dstep2 = 10E-8;
  else
    dstep2 = *step2;

  // Set the global function pointer
  FuncPtr = pt2Func;
  // Construct a Search object.
  EdHJSearch Search(*dim, minVec, dstep1, dstep2, f_fcn, NULL);
  
  // Set the maximum number of calls
  Search.SetMaxCallsExact(*max);
    
  // BeginSearch() may go on indefinitely for a poorly defined problem
  Search.BeginSearch();
  
  Search.GetMinPoint(Sminimum);
  Search.GetMinVal(SMinVal);
  Scalls = Search.GetFunctionCalls();

  // Pass back the final step length, minimum function evaluation,
  // and the corresponding x
  aux[0] = Search.GetDelta()*2.0;
  aux[1] = SMinVal;
  for (int i=0; i<*dim;i++)
    aux[i+2] = Sminimum[i];

  // Pass back the number of function evaluations it took
  *iaux = Scalls;
  // If BeginSearch() stopped, assume success
  *istat = 1;
  
  return;
}

/* This method sets up a simplex search based on the method
   of Spendley, Hext, and Himsworth given parameters passed 
   in through FORTRAN or C subroutines*/
void shhsch_(int *dim, double *start, int * right, double *step1,
             double *step2, int *max, int *stop, int *istat, double *aux,
             int *iaux, void (*pt2Func)(double *, double *, long *) )
{
  double SMinVal;         // to hold the min value
  long Scalls;            // to hold the number of calls
  
  Vector<double> Sminimum(*dim);   // to store the minimum point later
  
  // initialize an n-entry Vector with values start
  Vector<double>minVec(*dim, start);
    
  //Check for defaults
  double* dstep1 = new double[*dim];
  double dstep2;
  if (step1[0] == -1.0)
  {
    for (int j = 0; j < *dim; j++)
      dstep1[j] = 2.0;
  }
  else
  {
    for (int j = 0; j < *dim; j++)
      dstep1[j] = step1[j];
  }
  if (*step2 == -1.0)
    dstep2 = 10E-8;
  else
    dstep2 = *step2;

  // Set the global function pointer
  FuncPtr = pt2Func;
  // Construct a Search object.
  SHHSearch Search(*dim, minVec, dstep1[0], dstep2, f_fcn, NULL);
  
  // Choose simplex (Default is right)
  if (!*right)
  {
    Search.ChooseRegularSimplex();
  }
  else
  {
    // Make sure edge lengths are set correctly
    // (Variable-length right simplex is not automatically constructed)
    Vector<double>temp(*dim, dstep1);
    Search.SetStartingEdgeLengths(temp);
  }

  // Set the maximum number of calls
  Search.SetMaxCallsExact(*max);

  // If user wants to stop on standard deviation
  if (*stop)
    Search.Set_Stop_on_std();
    
  // BeginSearch() may go on indefinitely for a poorly defined problem
  Search.BeginSearch();
  
  Search.GetMinPoint(Sminimum);
  Search.GetMinVal(SMinVal);
  Scalls = Search.GetFunctionCalls();

  // Pass back the final step length, minimum function evaluation,
  // and the corresponding x
  aux[0] = Search.GetDelta();
  aux[1] = SMinVal;
  for (int i=0; i<*dim;i++)
    aux[i+2] = Sminimum[i];

  // Pass back the number of function evaluations it took
  *iaux = Scalls;
  // If BeginSearch() stopped, assume success
  *istat = 1;
  
  return;
}

/* This method sets up a simplex search based on the method
   of Nelder and Meade given parameters passed in through 
   FORTRAN or C subroutines*/
void nmsch_(int *dim, double *start, int *right, double *sigma,
            double *alpha, double *beta, double*gamma, double *step1,
            double *step2, int *max, int *stop, int *istat, double *aux,
            int *iaux, void (*pt2Func)(double *, double *, long *) )
{
  double SMinVal;         // to hold the min value
  long Scalls;            // to hold the number of calls
  
  Vector<double> Sminimum(*dim);   // to store the minimum point later
  
  // initialize an n-entry Vector with values start
  Vector<double>minVec(*dim, start);

  //Check for defaults
  double dsigma, dalpha, dbeta, dgamma;
  if (*sigma == -1.0)
    dsigma = 0.5;
  else
    dsigma = *sigma;
  if (*alpha == -1.0)
    dalpha = 1.0;
  else
    dalpha = *alpha;
  if (*beta == -1.0)
    dbeta = 0.5;
  else
    dbeta = *beta;
  if (*gamma == -1.0)
    dgamma = 2.0;
  else
    dgamma = *gamma;
  
  double* dstep1 = new double[*dim];
  double dstep2;
  if (step1[0] == -1.0)
  {
    for (int j = 0; j < *dim; j++)
      dstep1[j] = 2.0;
  }
  else
  {
    for (int j = 0; j < *dim; j++)
      dstep1[j] = step1[j];
  }
  if (*step2 == -1.0)
    dstep2 = 10E-8;
  else
    dstep2 = *step2;

  // Set the global function pointer
  FuncPtr = pt2Func;
  // Construct a Search object.
  NMSearch Search(*dim, minVec, dsigma, dalpha, dbeta, dgamma, dstep1[0],
                  dstep2, f_fcn, NULL);
  
  // Choose simplex (Default is right)
  if (!*right)
  {
    Search.ChooseRegularSimplex();
  }
  else
  {
    // Make sure edge lengths are set correctly
    // (Variable-length right simplex is not automatically constructed)
    Vector<double>temp(*dim, dstep1);
    Search.SetStartingEdgeLengths(temp);
  }

  // Set the maximum number of calls
  Search.SetMaxCallsExact(*max);

  // If user wants to stop on standard deviation
  if (*stop)
    Search.Set_Stop_on_std();
    
  // BeginSearch() may go on indefinitely for a poorly defined problem
  Search.BeginSearch();
  
  Search.GetMinPoint(Sminimum);
  Search.GetMinVal(SMinVal);
  Scalls = Search.GetFunctionCalls();

  // Pass back the final step length, minimum function evaluation,
  // and the corresponding x
  aux[0] = Search.GetDelta();
  aux[1] = SMinVal;
  for (int i=0; i<*dim;i++)
    aux[i+2] = Sminimum[i];

  // Pass back the number of function evaluations it took
  *iaux = Scalls;
  // If BeginSearch() stopped, assume success
  *istat = 1;
  
  return;
}

/* This method sets up a sequential multidirectional search
   given parameters passed in through FORTRAN or C subroutines*/
void smdsch_(int *dim, double *start, int * right, double *step1,
             double *step2, int *max, int *stop, int *istat, double *aux,
             int *iaux, void (*pt2Func)(double *, double *, long *) )
{
  double SMinVal;         // to hold the min value
  long Scalls;            // to hold the number of calls
  
  Vector<double> Sminimum(*dim);   // to store the minimum point later
  
  // initialize an n-entry Vector with values start
  Vector<double>minVec(*dim, start);

  //Check for defaults
  double* dstep1 = new double[*dim];
  double dstep2;
  if (step1[0] == -1.0)
  {
    for (int j = 0; j < *dim; j++)
      dstep1[j] = 2.0;
  }
  else
  {
    for (int j = 0; j < *dim; j++)
      dstep1[j] = step1[j];
  }
  if (*step2 == -1.0)
    dstep2 = 10E-8;
  else
    dstep2 = *step2;

  // Set the global function pointer
  FuncPtr = pt2Func;
  // Construct a Search object.
  SMDSearch Search(*dim, minVec, dstep1[0], dstep2, f_fcn, NULL);
  
  // Choose simplex (Default is right)
  if (!*right)
  {
    Search.ChooseRegularSimplex();
  }
  else
  {
    // Make sure edge lengths are set correctly
    // (Variable-length right simplex is not automatically constructed)
    Vector<double>temp(*dim, dstep1);
    Search.SetStartingEdgeLengths(temp);
  }

  // Set the maximum number of calls
  Search.SetMaxCallsExact(*max);

  // If user wants to stop on standard deviation
  if (*stop)
    Search.Set_Stop_on_std();
    
  // BeginSearch() may go on indefinitely for a poorly defined problem
  Search.BeginSearch();
  
  Search.GetMinPoint(Sminimum);
  Search.GetMinVal(SMinVal);
  Scalls = Search.GetFunctionCalls();

  // Pass back the final step length, minimum function evaluation,
  // and the corresponding x
  aux[0] = Search.GetDelta();
  aux[1] = SMinVal;
  for (int i=0; i<*dim;i++)
    aux[i+2] = Sminimum[i];

  // Pass back the number of function evaluations it took
  *iaux = Scalls;
  // If BeginSearch() stopped, assume success
  *istat = 1;
  
  return;
}


/* Need a small wrapper for the fortran main program, since fortran
   code should not be used as the main program. For more info, see:
   http://www.chiralcomp.com/support/mixing_f77_c_cpp/c++interface.html
*/
extern "C" int mymain_();
int main()
{
  return mymain_();
}
