/*
	A simple program to test the global optimizer.
*/

#include "global.h"
#include "tools.h"
#include "testfun.h"

#define STRLEN_MAX 80

int main() {
  bool AVfail, AVflag;
  int testfnc, dim, axis, i;
  double AVbest;
  Pdom  Dom;
  Pgrad Grad;
  Pobj  Obj;

  cout << "Global Optimizer v0.2" << endl;
  cout << "1.  Rosenbrock (min=0)" << endl;
  cout << "2.  McCormick (min=-1.91)" << endl;
  cout << "3.  Box & Betts (min=0)" << endl;
  cout << "4.  Paviani (min=-45.7)" << endl;
  cout << "5.  Generalized Rosenbrock" << endl;
  cout << "6.  Goldstein and Price" << endl;
  cout << "7.  Shekel 5" << endl;
  cout << "8.  Shekel 7" << endl;
  cout << "9.  Shekel 10" << endl;
  cout << "10. Levy4 (min=-21.502)" << endl;
  cout << "11. Levy5 (min=-11.504)" << endl;
  cout << "12. Levy6 (min=-11.504)" << endl;
  cout << "13. Levy7 (min=-11.504)" << endl;
  cout << "14. Griewank (min=0)" << endl;
  cout << "15. Six-hump Camel (min=-1.03) " << endl;
  cout << "16. Branin" << endl;
  cout << "17. Shubert (min=-24.06)" << endl;
  cout << "18. Hansen (min=-176.54)" << endl;
  cout << "19. Cola" << endl;
  cout << "20. Fitting problem" << endl;
  cout << endl << "Select test function :";
  cin >> testfnc;

  switch (testfnc) {
  case 1:
    dim=2;
    Dom=Domain_Rosenbrock;
    Obj=Objective_Rosenbrock;
    Grad=Gradient_Rosenbrock;
    break;
  case 2:
    dim=2;
    Dom=Domain_McCormick;
    Obj=Objective_McCormick;
    Grad=Gradient_McCormick;
    break;
  case 3:
    dim=3;
    Dom=Domain_BoxBetts; 
    Obj=Objective_BoxBetts;
    Grad=Gradient_BoxBetts; 
    break;
  case 4:
    dim=10;
    Dom=Domain_Paviani;
    Obj=Objective_Paviani;
    Grad=Gradient_Paviani;
    break;
  case 5:
    cout << "Enter problem dimension ";
    int ros_dim;   
    cin >> ros_dim;
    dim=ros_dim;
    Dom=Domain_Generalized_Rosenbrock;   
    Obj=Objective_Generalized_Rosenbrock;
    Grad=Gradient_Generalized_Rosenbrock;
    break;
  case 6:
    dim=2;
    Dom=Domain_Goldstein_Price;
    Obj=Objective_Goldstein_Price;
    Grad=Gradient_Goldstein_Price;
    break;
  case 7:
    dim=4;
    Dom=Domain_Shekel;
    Obj=Objective_Shekel5;
    Grad=Gradient_Shekel5;
    break;
  case 8:
    dim=4;
    Dom=Domain_Shekel;
    Obj=Objective_Shekel7;
    Grad=Gradient_Shekel7;
    break;
  case 9:
    dim=4;
    Dom=Domain_Shekel;
    Obj=Objective_Shekel;
    Grad=Gradient_Shekel;
    break;
  case 10:
    dim=4;
    Dom=Domain_Levy;
    Obj=Objective_Levy;
    Grad=Gradient_Levy;  
    break;
  case 11:
    dim=5;
    Dom=Domain_Levy; 
    Obj=Objective_Levy; 
    Grad=Gradient_Levy; 
    break;
  case 12:
    dim=6;
    Dom=Domain_Levy;
    Obj=Objective_Levy;
    Grad=Gradient_Levy;
    break;
  case 13:
    dim=7;
    Dom=Domain_Levy;
    Obj=Objective_Levy;
    Grad=Gradient_Levy;
    break;
  case 14:
    dim=10;
    Dom=Domain_Griewank;
    Obj=Objective_Griewank;
    Grad=Gradient_Griewank;
    break;
  case 15:
    dim=2;
    Dom=Domain_Camel6;
    Obj=Objective_Camel6;
    Grad=Gradient_Camel6;
    break;
  case 16:
    dim=2;
    Dom=Domain_Branin;
    Obj=Objective_Branin;
    Grad=Gradient_Branin;
    break;
  case 17:
    dim=2;
    Dom=Domain_Shubert;
    Obj=Objective_Shubert;
    Grad=Gradient_Shubert;
    break;
  case 18:
    dim=2;
    Dom=Domain_Hansen;
    Obj=Objective_Hansen;
    Grad=Gradient_Hansen;
    break;
  case 19:
    dim=17;
    Dom=Domain_Cola;
    Obj=Objective_Cola;
    Grad=Gradient_Cola;
    break;
  case 20:
    dim=12;
    Dom=Domain_myMinHeight;
    Obj=Objective_myMinHeight;
    Grad=Gradient_myMinHeight;
    break;
  default:
    cout << "Error : Function not defined" << endl;
    exit(1);
  }

  cout << "Dimension=" <<dim << endl;
  TBox D(dim);
  Dom(D);
  cout << "Domain=";
  for (i=0; i<dim; i++)
    cout << "[" << D.lb(i) << "," << D.ub(i) << "]";
  cout << endl << endl;

  GlobalParams params;
  cout << "Enter time limit (seconds) ";
  cin >> params.maxtime;
  if (params.maxtime<1) {   
    cout << "Warning: time limit set to 1 second\n";
  }
//    cout << "Number of deterministic points ";
//    cin >> params.det_pnts;

  cout << "Use factory settings (y/n) ";
  char str[STRLEN_MAX]; cin >> str;
  if (str[0]=='y') {
    params.det_pnts=2*dim+1; 
    params.rnd_pnts=0;
    params.eps_cl=0.1; params.rshift=0.3;
    params.mu=1.0E-4; AVflag=FALSE;
  }
  else {
    cout << "Number of deterministic points ";
    cin >> params.det_pnts;
    cout << "Numer of stochastic points ";
    cin >> params.rnd_pnts;
    cout << "Radius of attraction ";
    cin >> params.eps_cl;
    cout << "Parameter rshift ";
    cin >> params.rshift;
    cout << "Parameter mu ";
    cin >> params.mu;
    cout << "Use the AV initialization (y/n) ";
    cin >> str;
    if (str[0]=='y') AVflag=TRUE; else AVflag=FALSE;
  }

  Global Problem(D,Obj, Grad, params);
  RVector x_av(dim);
  if (AVflag==TRUE) {
    cout << "Enter time limit for each coordinate direction (seconds) ";
    cin >> params.maxtime;
    if (params.maxtime<1) {
      cout << "Warning: time limit set to 1 second\n";
    }
    params.det_pnts=3;
    TBox I(1);
    Global AV(I, Obj, Grad, params);

    x_av=0.0; AVfail=FALSE;
    for (axis=0; axis<Problem.dim; axis++) {
      cout << "### axis=" << axis << " ###" << endl;
      I.lb=(D.lb)(axis); I.ub=(D.ub)(axis);
      AV.SetDomain(I);
      AV.ClearSolSet();
      AV.Search(axis, x_av);

      if (AV.NoMinimizers()) {
	cout << "AV failed with axis=" << axis << endl;
	AVfail=TRUE; break;
      }
    }

    if (AVfail==FALSE) {
      AVbest=AV.GetMinValue();
      cout << "### AV Located x=" << x_av << " fbound=" << AVbest << endl;
      RVector AVx(Problem.dim);
      AVx=x_av;

      // Use result from AV for new fbound
      Problem.SetMinValue(AVbest);

      // Add the best point found to the initial box (domain)
      Problem.AddPoint(x_av, AVbest);      
    }
  }

  // Perform the main search
  cout << "### Starting main search ###" << endl;
  Problem.Search(-1, x_av);

  cout << "Optimization terminated. Current set of minimizers is" << endl;
  if (Problem.NoMinimizers() && AVflag==FALSE)
    cout << "### No improvement found ###" << endl;
  else
    Problem.DispMinimizers();
}