/***********************************************************************/
/*
//     This is an example call of MIDACO 6.0
//     -------------------------------------
//
//     MIDACO solves Multi-Objective Mixed-Integer Non-Linear Problems:
//
//
//      Minimize     F_1(X),... F_O(X)  where X(1,...N-NI)   is CONTINUOUS
//                                      and   X(N-NI+1,...N) is DISCRETE
//
//      subject to   G_j(X)  =  0   (j=1,...ME)      equality constraints
//                   G_j(X) >=  0   (j=ME+1,...M)  inequality constraints
//
//      and bounds   XL <= X <= XU
//
//
//     The problem statement of this example is given below. You can use 
//     this example as template to run your own problem. To do so: Replace 
//     the objective functions 'F' (and in case the constraints 'G') given 
//     here with your own problem and follow the below instruction steps.
*/
/***********************************************************************/
#include <omp.h> /* OpenMP-Library */
#include <stdlib.h>
#include <stdio.h>

#include <math.h>
#include <time.h>

double gettime(){time_t second;second=time(NULL);return (double)second;}

/***********************************************************************/
void problem_function( double *F, double *G, double *X, int GTOP_PROBLEM, long int n);
/***********************************************************************/
extern"C"{ int midaco(long int*,long int*,long int*,long int*,long int*,
                      long int*,double*,double*,double*,double*,double*,
                      long int*,long int*,double*,double*,long int*,
                      long int*,long int*,double*,long int*,char*);}
/***********************************************************************/
extern"C"{ int midaco_print(int,long int,long int,long int*,long int*,double*,
                            double*,double*,double*,double*,long int,long int,
                            long int,long int,long int,double*,double*,
                            long int,long int,double*,long int,char*);}
/***********************************************************************/
/************************   MAIN PROGRAM   *****************************/
/***********************************************************************/
int main() 
{ 
      /*
          1) Cassini1
          2) GTOC1
          3) Messenger_Full
          4) Messenger_Reduced
          5) Cassini2
          6) Rosetta
          7) Sagas 
      */
      int GTOP_PROBLEM = 2; /* Select GTOP Benchmark Problem */

      /* Variable and Workspace Declarations */
      long int o,n,ni,m,me,maxeval,maxtime,printeval,save2file,iflag,istop;
      long int liw,lrw,lpf,i,iw[50000],p=1; double rw[50000],pf[50000];
      double   f[10],g[1000],x[1000],xl[1000],xu[1000],param[13];
      char key[] = "************************************************************";

      /*****************************************************************/
      /***  Step 1: Problem definition  ********************************/
      /*****************************************************************/

      /* STEP 1.A: Problem dimensions
      ******************************/
      o  = 1; /* Number of objectives                          */

      /* Adjust dimension n and m to specific GTOP Problem */
      if(GTOP_PROBLEM == 1){n= 6; m=4;}
      if(GTOP_PROBLEM == 2){n= 8; m=6;}
      if(GTOP_PROBLEM == 3){n=26;     }
      if(GTOP_PROBLEM == 4){n=18;     }
      if(GTOP_PROBLEM == 5){n=22;     }
      if(GTOP_PROBLEM == 6){n=22;     }
      if(GTOP_PROBLEM == 7){n=12; m=2;}    
  
      /*  STEP 1.B: Define lower and upper bounds 'xl' & 'xu'  */
      /*  ***************************************************  */
      double pi = 3.14159265359;   
      double fex = 0.0; /* initialize best known value */
      
      if( GTOP_PROBLEM == 1 ){  
        fex = 4.93075;
        double lower[] = {-1000,30,100,30,400,1000};
        double upper[] = {0,400,470,400,2000,6000};
        for( i=0; i<n; i++){ 
        xl[i] = lower[i]; 
        xu[i] = upper[i]; };
      } 
          
      if( GTOP_PROBLEM == 2 ){  
        fex = -1581950.0;     
        double lower[] = {3000, 14, 14, 14, 14, 100, 366, 300};
        double upper[] = {10000, 2000, 2000, 2000, 2000, 9000, 9000, 9000};
        for( i=0; i<n; i++){ 
        xl[i] = lower[i]; 
        xu[i] = upper[i]; };
      }  
      
      if( GTOP_PROBLEM == 3 ){   
        fex = 2.104;    
        double lower[] = {1900,2.5,0,0,100,100,100,100,100,100,0.01,0.01,0.01,0.01,0.01,0.01,1.1,1.1,1.05,1.05,1.05,  -pi,  -pi,  -pi,  -pi,  -pi};
        double upper[] = {2300,4.05,1,1,500,500,500,500,500,600,0.99,0.99,0.99,0.99,0.99,0.99,6,6,6,6,6,   pi,  pi,  pi,  pi,  pi};
        for( i=0; i<n; i++){ 
        xl[i] = lower[i]; 
        xu[i] = upper[i]; };
      }  
      
      if( GTOP_PROBLEM == 4 ){   
        fex = 8.630;    
        double lower[] = {1000,1,0,0,30,30,30,30,0.01,0.01,0.01,0.01,1.1,1.1,1.1,-pi,-pi,-pi};
        double upper[] = {4000,5,1,1,400,400,400,400,0.99,0.99,0.99,0.99,6,6,6,pi,pi,pi};
        for( i=0; i<n; i++){ 
        xl[i] = lower[i]; 
        xu[i] = upper[i]; };
      }  
      
      if( GTOP_PROBLEM == 5 ){   
        fex = 8.383;    
        double lower[] = {-1000,3,0,0,100,100,30,400,800,0.01,0.01,0.01,0.01,0.01,1.05,1.05,1.15,1.7,-pi,-pi,-pi,-pi};
        double upper[] = {0,5,1,1,400,500,300,1600,2200,0.9,0.9,0.9,0.9,0.9,6,6,6.5,291,pi,pi,pi,pi};
        for( i=0; i<n; i++){ 
        xl[i] = lower[i]; 
        xu[i] = upper[i]; };
      }  
      
      if( GTOP_PROBLEM == 6 ){   
        fex = 1.343;    
        double lower[] = {1460,3,0,0,300,150,150,300,700,0.01,0.01,0.01,0.01,0.01,1.06,1.05,1.05,1.05,-pi,-pi,-pi,-pi};
        double upper[] = {1825,5,1,1,500,800,800,800,1850,0.9,0.9,0.9,0.9,0.9,9,9,9,9,pi,pi,pi,pi};
        for( i=0; i<n; i++){ 
        xl[i] = lower[i]; 
        xu[i] = upper[i]; };
      }                     
      
      if( GTOP_PROBLEM == 7 ){   
        fex = 18.194;    
        double lower[] = {7000, 0, 0, 0,   50,  300, 0.01, 0.01, 1.05,   8, -pi, -pi};
        double upper[] = {9100, 7, 1, 1, 2000, 2000,  0.9,  0.9,    7, 500,  pi,  pi};
        for( i=0; i<n; i++){ 
        xl[i] = lower[i]; 
        xu[i] = upper[i]; };
      }    

      /* STEP 1.C: Starting point 'x'  
      ******************************/     
      for( i=0; i<n; i++)
      { 
         x[i] = xl[i]; /* Here for example: starting point = lower bounds */
      } 

      /*****************************************************************/
      /***  Step 2: Choose stopping criteria and printing options   ****/
      /*****************************************************************/

      /* STEP 2.A: Stopping criteria 
      *****************************/
      maxeval = 999999999;    /* Maximum number of function evaluation (e.g. 1000000)  */
      maxtime = 999999999;//60*60*24; /* Maximum time limit in Seconds (e.g. 1 Day = 60*60*24) */

      /* STEP 2.B: Printing options  
      ****************************/
      printeval = 1000000; /* Print-Frequency for current best solution (e.g. 1000) */
      save2file = 1;    /* Save SCREEN and SOLUTION to TXT-files [ 0=NO/ 1=YES]  */
    
      /*****************************************************************/
      /***  Step 3: Choose MIDACO parameters (FOR ADVANCED USERS)    ***/
      /*****************************************************************/
      param[ 0] =  0.0;  /* ACCURACY  */
      param[ 1] =  4.0;  /* SEED      */
      param[ 2] =  fex+abs(fex)/1000.0;  /* FSTOP     */
      param[ 3] =  0.0;  /* ALGOSTOP  */
      param[ 4] =  0.0;  /* EVALSTOP  */
      param[ 5] =  0.0;  /* FOCUS     */
      param[ 6] =  0.0;  /* ANTS      */
      param[ 7] =  0.0;  /* KERNEL    */
      param[ 8] =  0.0;  /* ORACLE    */
      param[ 9] =  0.0;  /* PARETOMAX */
      param[10] =  0.0;  /* EPSILON   */
      param[11] =  0.0;  /* BALANCE   */
      param[12] =  0.0;  /* CHARACTER */

      /*****************************************************************/
      /***  Step 4: Choose Parallelization Factor    *******************/
      /*****************************************************************/  

      p = 4; /* Serial = 1, Parallel = 2,3,4,... */

      /*****************************************************************/
      /*****************************************************************/
      /*****************************************************************/            
      /* Preparations for Parallelization */
      int nthreads,c;   
      double *xxx,*fff,*ggg; 
      #pragma omp parallel private(c) 
      nthreads = omp_get_num_threads(); omp_set_num_threads(p);
      printf("\n [ Available threads on this machine: %i ] \n",nthreads);
      /* Allocate arrays for parallelization */
      fff = (double *) malloc((p*o)*sizeof(double));
      ggg = (double *) malloc((p*m)*sizeof(double));
      xxx = (double *) malloc((p*n)*sizeof(double));     
      /* Store starting point x in xxx array */
      for( c=0; c<p; c++){ for( i=0; i<n; i++){ xxx[c*n+i] = x[i]; }}        
      /*****************************************************************/
      /*   
         Call MIDACO by Reverse Communication
      */
      /*****************************************************************/
      /* Workspace length calculation */
      lrw=sizeof(rw)/sizeof(double); 
      lpf=sizeof(pf)/sizeof(double);   
      liw=sizeof(iw)/sizeof(long int);     
      /* Print midaco headline and basic information */
      midaco_print(1,printeval,save2file,&iflag,&istop,&*f,&*g,&*x,
        &*xl,&*xu,o,n,ni,m,me,&*rw,&*pf,maxeval,maxtime,&*param,p,&*key);
      while(istop==0) /*~~~ Start of the reverse communication loop ~~~*/
      {   
        /* Evaluate block of 'P' iterates in parallel */
        #pragma omp parallel for                             
        for (c=0; c<p; c++)
        { 
          problem_function( &fff[c], &ggg[c*m], &xxx[c*n], GTOP_PROBLEM, n);         }  
                           
          /* Call MIDACO */
          midaco(&p,&o,&n,&ni,&m,&me,&*xxx,&*fff,&*ggg,&*xl,&*xu,&iflag,
                 &istop,&*param,&*rw,&lrw,&*iw,&liw,&*pf,&lpf,&*key);                  
          /* Call MIDACO printing routine */            
          midaco_print(2,printeval,save2file,&iflag,&istop,&*fff,&*ggg,&*xxx,
            &*xl,&*xu,o,n,ni,m,me,&*rw,&*pf,maxeval,maxtime,&*param,p,&*key);   
      } /*~~~End of the reverse communication loop ~~~*/  
      /*****************************************************************/
      // printf("\n Solution f[0] = %f ", fff[0]);
      // printf("\n Solution g[0] = %f ", ggg[0]);
      // printf("\n Solution x[0] = %f ", xxx[0]);     
      // for( i=0; i<n; i++){ x[i] = xxx[i]; } /* Retrieve x from xxx */       
      /*****************************************************************/      
      return 0;
      /*****************************************************************/
      /***********************  END OF MAIN  ***************************/
      /*****************************************************************/ 
}



/***********************************************************************/
/*********************   OPTIMIZATION PROBLEM   ************************/
/***********************************************************************/
#include "trajobjfuns.h"
#include <vector>
#define VENUS 2
#define EARTH 3
#define MARS 4
#define JUPITER 5
#define SATURN 6
/***********************************************************************/
/***********************************************************************/
/***********************************************************************/
void problem_function(double *f, double *g, double *x, int GTOP_PROBLEM, long int n)
{
  double DVtot, DVonboard;
  std::vector<double> solution, rp;   
  for (int i=0; i<n; ++i){ solution.push_back(x[i]); } 

  if(GTOP_PROBLEM == 1)
  { 
/* Objective function value F(X) (denoted here as f[0]) */
   f[0] = cassini1(solution,rp);   
/* Constraints */
   g[0] = rp[0] - 6351.8; 
   g[1] = rp[1] - 6351.8;
   g[2] = rp[2] - 6778.1;
   g[3] = rp[3] - 671492.0;
  }
  
 
  if(GTOP_PROBLEM == 2)
  { 
/* Objective function value F(X) (denoted here as f[0]) */
   f[0] = gtoc1(solution,rp) - 2000000.0; 
/* Constraints */
   g[0] = rp[0] - 6351.8;
   g[1] = rp[1] - 6778.1;
   g[2] = rp[2] - 6351.8;
   g[3] = rp[3] - 6778.1;
   g[4] = rp[4] - 600000.0;
   g[5] = rp[5] - 70000.0;
  }
    
  
  if(GTOP_PROBLEM == 3){ f[0] = messengerfull(solution); }
  if(GTOP_PROBLEM == 4){ f[0] = messenger(solution); } 
  if(GTOP_PROBLEM == 5){ f[0] = cassini2(solution); } 
  if(GTOP_PROBLEM == 6){ f[0] = rosetta(solution); }        
 
 
  if(GTOP_PROBLEM == 7)
  { 
/* Objective function value F(X) (denoted here as f[0]) */
   f[0] = sagas(solution,DVtot,DVonboard);   
/* Constraints */
   g[0] = 6.782 - DVtot;
   g[1] = 1.782 - DVonboard;
  }       
}
/***********************************************************************/
/*************************** END OF FILE *******************************/
/***********************************************************************/