/**
  * This program conducts a genetic algorithm on "chromosomes" of
  * 6 real parameters a1, b1, c1, a2, b2, c2 maximizing whatever
  * function is supplied at the method "objective".  Here that function
  * is supplied by the java class amerputExerBoundary; a program that
  * is outlined in the text as Algorithm 14 in Section 3.5.6.
  *
  * Created on 9/3/2010
  *
  * (c) Copyright Ronald W. Shonkwiler
  * All rights reserved.
  *
  * RSS style GA to find the value of an American put option
  * This will be the maximum option value over all possible
  * exercise boundary curves.
  *  Each curve is the sum of two parameterized log functions
  *
  * Adjust the parameters dfS0, dfX, ... just below using an editor
  * and recompile.
  */
import java.util.*;
import static java.lang.System.err;

public class
amerputGAweb
{
  static int nIterations= 20; //120; //GA generations
  static int nTrials = 3;

  static double dfS0 = 100; //default S0
  static double dfX = 100;  //default X
  static int dfnDays = 40;  //default T in days
  static double dfrfr = 0.26; //risk-free rate
  static double dfvol = 0.40; //volatility
  static double dtINdays = 5; //1 for daily, could be fractional

  static double daysInYr = 365;
  static String TinYears =
    new String(String.format("%.2f",dfnDays/daysInYr));
  static String yearVSyears = new String(
    (TinYears.equals("1.0")?"year":"years"));

  static String[] helptext =
  {
"Implementation of the RSS style GA for maximizing an American.",
"put option value over parameterized exercise boundary curves.",
   "Using the following parameters (change in source and recompile)",
   "  -S<stock price>     default "+dfS0+"",
   "  -X<strike price>    default "+dfX+"",
   "  -T<time to expiry>  (days) default "+dfnDays+
     " ("+TinYears+" "+yearVSyears+")",
   "  -v<volatility>      default "+dfvol+" (annualized)",
   "  -r<riskfree rate>   default= "+dfrfr+" (annualized)",
   "  -d<dtInDays>        time increment(days), default "+dtINdays+"",
   "",
  };

//=== problem specifics
  static double aBot = 0, aTop = 0.1; //range of a parameter
  static double bBot = 0.1, bTop = 4.0;
  static double cBot = 0.05, cTop = 1.5;

    //pertains to all instances
  static final int NCEXIT=60; //no change exit condition
  static Random rng;
  static long seed;
    //command line variable
  static int vfreq=5; //vfreq=0;

  static ArrayList minPerformance = new ArrayList();
  static ArrayList maxPerformance = new ArrayList();

  static final int BUGINVALID=-1; //mark popspace free
  static final double tempSTART=0.88; //living prob. 'temperature'

  static final int STDZ=12; //population size
  static final int HALFSTDZ=STDZ/2;  //half std pop size
  static final int THREEQSTD=3*STDZ/4;  // 3 quarters std pop size
  static final int POPSPACE=3*STDZ;  //space for critters
  static final int MAXMATE=THREEQSTD;  //max number of matings
  static final int MAXMUTE=HALFSTDZ; // max number of mutations

//=== genetic specifics
    //a chromosome here is a list of 6 real values between 
    //given limits, see Bot and Top
  static double mutLimit = 0.005; //mutate up to this much
  int popSize = STDZ;

//ggg===
  critterAP[] c0 = new critterAP[POPSPACE]; //ptrs to critters
  critterAP bestCrit = new critterAP(this);
  critterAP rerunCrit = new critterAP(this);
  int[] role = new int[POPSPACE];

  long fevals=0L; //keep track of fcn evaluations
  boolean logProg = false;
  double logBase = 0; //only log fitnesses exceeding this

  double S0, X, vol, rfr;
  int nDays;
  amerputExerBoundary ameb = null;

  double besta1=0, bestb1=0, bestc1=0;
  double besta2=0, bestb2=0, bestc2=0;
//ccc===
  public
  amerputGAweb(double s0, double r, double v, double x, int ndays)
  {
    S0 = s0; rfr = r; vol = v; X = x; nDays = ndays;

    seed = System.currentTimeMillis();
//      seed = 1174741248177L; // uncomment for testing
    rng = new Random(seed);
    System.out.println("\nseed: "+seed);
    System.out.println("amerputGAweb@cons dtINdays= "+dtINdays);

      //set up the option parameters
    ameb = new amerputExerBoundary(S0,rfr,vol,X,nDays,rng);
    ameb.setdtINdays(dtINdays);

    double bestOverAll = 0;
    double runOut;
    ArrayList[] runLogs = new ArrayList[nTrials];

      //initialize the min/max record
    for(int i=1; i <= nIterations; ++i)
    {
      if( vfreq!=0 && i%vfreq == 0 )
      {
        minPerformance.add(Double.MAX_VALUE);
        maxPerformance.add((Double)0.);
      }
    }

    double[] avgNvar = new double[2];
    for( int i=0; i<nTrials; ++i)
    {
      if( i>0 )
      {
        logProg = true;
        logBase = bestOverAll;
      }

      runLogs[i] = new ArrayList();
      runOut = runGA(runLogs[i]);

      if( runOut > bestOverAll )
      {
         bestOverAll = runOut;
         besta1 = bestCrit.a1;
         bestb1 = bestCrit.b1;
         bestc1 = bestCrit.c1;
         besta2 = bestCrit.a2;
         bestb2 = bestCrit.b2;
         bestc2 = bestCrit.c2;
      }
err.format("trial %d, returnValue = %f\n",i,runOut);
err.format("  a1= %.4f\n",bestCrit.a1);
err.format("  b1= %.4f\n",bestCrit.b1);
err.format("  c1= %.4f\n",bestCrit.c1);
err.format("  a2= %.4f\n",bestCrit.a2);
err.format("  b2= %.4f\n",bestCrit.b2);
err.format("  c2= %.4f\n",bestCrit.c2);
err.format("Rerunning the objective with these parameters:");
  avgNvar = objective(
bestCrit.a1,bestCrit.b1,bestCrit.c1,
bestCrit.a2,bestCrit.b2,bestCrit.c2);
err.format("  result= %.3f, var= %.4f\n",avgNvar[0],avgNvar[1]);
    }

/*
      for( int i=0; i<nTrials; ++i)
      {
        System.out.println("run "+i+" record");
        for( int t=0; t<runLogs[i].size(); ++t)
        {
          System.out.format("t= %3d, best= %.3f\n",
            t,((Double)runLogs[i].get(t)).doubleValue());
        }
      }

      System.out.println(
        "min/max record every "+vfreq+" over "+nTrials+" runs");
      for( int j=0; j< minPerformance.size(); ++j)
      {
        System.out.format(
          "j= %4d, min %.3f, max %.3f\n",(j+1)*vfreq,
          ((Double)minPerformance.get(j)).doubleValue(),
          ((Double)maxPerformance.get(j)).doubleValue());
      }
*/

    System.out.println(
      "in "+nTrials+" trials, best overall= "+bestOverAll+" parameters:");
    System.out.format(" a1= %.4f\n",besta1);
    System.out.format(" b1= %.4f\n",bestb1);
    System.out.format(" c1= %.4f\n",bestc1);
    System.out.format(" a2= %.4f\n",besta2);
    System.out.format(" b2= %.4f\n",bestb2);
    System.out.format(" c2= %.4f\n",bestc2);
  }
//===
    /**
    * need to control the variance in the individual
    * option price calculation
    */
  double[]
  objective( 
       double a1, double b1, double c1,
       double a2, double b2, double c2
    )
  {
    double avg=0, var=0;
    int nReps = 6;

    ameb.setEBparams(a1,b1,c1,a2,b2,c2);
    for( int i=0; i<nReps; ++i )
    {
      ++fevals;
      ameb.run();
      avg += ameb.discExpec;
      var += ameb.discExpec*ameb.discExpec;
    }
    avg /= nReps;
    var = var/(nReps-1) - avg*avg;

//err.format("expec= %.3f for params:\n",ameb.discExpec);
//err.format(" a1 %.4f b1 %.3f c1 %.3f, a2 %.4f b2 %.3f c2 %.3f\n",
//a1,b1,c1, a2,b2,c2);
    return new double[] {avg,var};
  }
//sss===
  double
  runGA(ArrayList runRecord)
  {
    int t = 1; //algorithm time for cooling/stopping purposes
    int z=0;  //current population size
    int noImprov=0; //no improvement count
    bestCrit.fitness = 0;

    double[] af = new double[1]; //for passback
    double temp; //temperature, using as RR living prob.

      //mark population slots free
    for( int i=0; i<POPSPACE; ++i )
    {
      c0[i] = new critterAP(this);
      c0[i].fitness = BUGINVALID;
    }

      //INITIALIZATION, create the population
    af[0]=0;
    z=0;
    for(int i=0; i<STDZ; ++i )
    {
        //random initialization
      c0[i].a1 = aBot + (aTop-aBot)*rng.nextDouble();
      c0[i].b1 = bBot + (bTop-bBot)*rng.nextDouble();
      c0[i].c1 = cBot + (cTop-cBot)*rng.nextDouble();
      c0[i].a2 = aBot + (aTop-aBot)*rng.nextDouble();
      c0[i].b2 = bBot + (bTop-bBot)*rng.nextDouble();
      c0[i].c2 = cBot + (cTop-cBot)*rng.nextDouble();
      c0[i].setFitness();
//err.format("init@ i %d, fitness= %.4f\n",i,c0[i].fitness);
      z = merge(-1,i,af,z);  /* sense = -1 for maximization */
    }

      //EVOLUTION LOOP
    int nmate = MAXMATE; //currently 6
    int nmute = MAXMUTE; //current 6
    int dltz=0;

    int ni, k; //int helpers
    int alIndex = 0;

    outerLoop: //block label
    for( t=1; t <= nIterations; ++t)
    {
      if( vfreq!=0 && t%vfreq == 0 )
      {
        report(1,0,t); //1=>only print value, 0=>only for best

        double min =
      ((Double)amerputGAweb.minPerformance.get(alIndex)).doubleValue();
        double max =
      ((Double)amerputGAweb.maxPerformance.get(alIndex)).doubleValue();

//err.format(
//"before: t= %2d, alIndex= %3d, min= %.3f, max= %.3f, best= %.3f\n",
//t,alIndex,min,max,bestCrit.fitness);

        if( bestCrit.fitness < min )
      amerputGAweb.minPerformance.set(alIndex,(Double)bestCrit.fitness);
        if( bestCrit.fitness > max )
      amerputGAweb.maxPerformance.set(alIndex,(Double)bestCrit.fitness);

        ++alIndex;
      }

//hhh
      for( int i=0; i<nmute; ++i)
      { 
        if( (ni = getslot()) < 0 ) break; //no free slots
          //pick random one to mutate
        k = rng.nextInt(z);
        z = mutate(c0[role[k]],ni,af,z);
      } 

        //mating (cross-over) step
        // mate top half, find mate at random (among current pop.)
      for( int m=0; m<nmate; ++m )
      {
        if((ni = getslot()) < 0 ) break; //no free slots
          //random population member
        k = rng.nextInt(z); 
        z = mate(new critterAP[] {c0[role[m]],c0[role[k]]},ni,af,z);
      }

//err.format("t= %d, bestfitness now %.3f, best of population %.3f\n",
//t,bestCrit.fitness,c0[role[0]].fitness);
        //check for improvement
      if( c0[role[0]].fitness > bestCrit.fitness )
      {
          //evaluate a second time, take the average of the two
        c0[role[0]].copy(rerunCrit);
        rerunCrit.setFitness();
        System.out.format("Update best: fitness_1= %.3f, fitness_2= %.3f\n",
         c0[role[0]].fitness,rerunCrit.fitness);

        rerunCrit.fitness = (c0[role[0]].fitness+rerunCrit.fitness)/2;

          //if space available, re-merge
        if( (ni = getslot()) >= 0 )
        {
          rerunCrit.copy(c0[ni]);
          z=removal(0,af,z);
          z=merge(-1,ni,af,z);
        }

        if( rerunCrit.fitness > bestCrit.fitness )
        { 
            //copy c0[role[0]] into best
          c0[role[0]].copy(bestCrit);
//err.format("new best: fitness= %.3f\n",bestCrit.fitness);
          noImprov=0;
            //check whether to print solution
          if( logProg && bestCrit.fitness > logBase )
          {
//report(1,0,t); //1=>only print value, 0=>only for best
            logBase = bestCrit.fitness;
          }
        }
      }
  
        //check for no change exit
      if( ++noImprov == NCEXIT ) break;

        //select down the population size
      for( dltz=STDZ-z; dltz<0; ++dltz)
        z = removal(gauntlet(tempSTART,z),af,z);

      runRecord.add((Double)bestCrit.fitness);
    }

//err.println("#GA generations done= "+t);
//report(2,0,t); //2=>print value & matrix, 0=>only for best
//System.out.println("\n#fevals "+fevals);

    return bestCrit.fitness;
  }
//===
  public int
  merge(int sense, int i, double[] af, int cs) 
  /*::::::::
    merge critter i into role according to rank via fitness
    sense= +1, small to big, i.e. minimization
    sense= -1, big to small, maximization
    i= index of this critter
    af= average fitness, cs= present population size
    return new population size
  :::::::*/
  {
    int m, k, j;
    int z = cs;

      //insert newbug according to its fitness
    for( m=0; ; ++m) {
      if( m == z ) break;
      if( sense*c0[role[m]].fitness > sense*c0[i].fitness ) break;
    }

      //goes into slot m, see if room and update avg fitness
    af[0] *= z;  //sum of fitnesses
    if( z < POPSPACE ) //have room, so add
    {
      ++z;
      af[0] = (af[0] + c0[i].fitness)/z;
    }
    else { //no room, knock off last
      af[0] += c0[i].fitness - c0[POPSPACE-1].fitness;
      af[0] = (af[0] + c0[i].fitness)/z;
    }

    for( j=z-2; j>=m; j--) role[j+1] = role[j];
    role[m] = i;

    return z;
  }
//===
  public int
  getslot()
  {
  int i;

    for(i=0; i<POPSPACE; ++i )
      if( c0[i].fitness == BUGINVALID ) return i;
    return -1;
  }
//rrr===
  void
  report( int rLevel, int todo, long its)
  {
    /*::::::
      rLevel= extent to report:
        0 minimal (just dots),
        1 just fitness values
        2 fitness values and entire matrix
      todo= # from top of gene pool to print
      its iteration count
    ::::::*/
    if( rLevel == 0 ) { System.out.print("."); return;}

    if( rLevel == 1 )
    {
      System.out.format(
        "best= %.3f, #its= %d\n",bestCrit.fitness,its);
      for( int i=0; i<todo; ++i)
        System.out.format("critter %4d fitness= %.3f\n",i,c0[i].fitness);
      return;
    }

    if( rLevel == 2 )
    {
      System.out.format("best= %.3f, #its= %d\n",bestCrit.fitness,its);
      for( int i=0; i<todo; ++i)
      {
        err.format(
      "critter %4d fitness= %.3f\n",i,c0[i].fitness);
      }
    }
  }
//===
    /**
    * crossover style mating, keep one offspring randomly
    */
  int //return updated value of z
  mate( critterAP[] cArray, int ni, double[] af, int z)
  {
     int k;
      //uniform cross-over among the 6 parameters
    for( int i=0; i<6; ++i )
    {
      if( rng.nextDouble() < 0.5 )
        k=0;//pick from critter 0
      else
        k=1;
      switch( i ) {
      case 0:
        c0[ni].a1 = cArray[k].a1;
        break;
      case 1:
        c0[ni].b1 = cArray[k].b1;
        break;
      case 2:
        c0[ni].a2 = cArray[k].a2;
        break;
      case 3:
        c0[ni].b2 = cArray[k].b2;
        break;
      case 4:
        c0[ni].c1 = cArray[k].c1;
        break;
      case 5:
        c0[ni].c2 = cArray[k].c2;
        break;
      }
    }

    c0[ni].setFitness();

//err.format("  off statex= %.4f, statey= %.4f fit: %.7f\n",
//c0[ni].statex,c0[ni].statey,c0[ni].fitness);

    return merge(-1,ni,af,z);
  }
//===
    /**
    * perturbate x and y up to +/- mutLimit
    */
  int
  mutate(critterAP crit, int ni, double[] af, int z)
  {
    double a, b;

      //copy crit into slot ni
    c0[ni].a1 = crit.a1;
    c0[ni].b1 = crit.b1;
    c0[ni].c1 = crit.c1;
    c0[ni].a2 = crit.a2;
    c0[ni].b2 = crit.b2;
    c0[ni].c2 = crit.c2;

      //mutate one of the 6
    int k = rng.nextInt(6);
    switch( k ) {
      case 0:
        a = Math.max(aBot,crit.a1-mutLimit);
        b = Math.min(aTop,crit.a1+mutLimit);
        c0[ni].a1 = a + (b-a)*rng.nextDouble();
        break;
      case 1:
        a = Math.max(bBot,crit.b1-mutLimit);
        b = Math.min(bTop,crit.b1+mutLimit);
        c0[ni].b1 = a + (b-a)*rng.nextDouble();
        break;
      case 2:
        a = Math.max(aBot,crit.a2-mutLimit);
        b = Math.min(aTop,crit.a2+mutLimit);
        c0[ni].a2 = a + (b-a)*rng.nextDouble();
        break;
      case 3:
        a = Math.max(bBot,crit.b2-mutLimit);
        b = Math.min(bTop,crit.b2+mutLimit);
        c0[ni].b2 = a + (b-a)*rng.nextDouble();
        break;
      case 4:
        a = Math.max(cBot,crit.c1-mutLimit);
        b = Math.min(cTop,crit.c1+mutLimit);
        c0[ni].c1 = a + (b-a)*rng.nextDouble();
        break;
      case 5:
        a = Math.max(cBot,crit.c2-mutLimit);
        b = Math.min(cTop,crit.c2+mutLimit);
        c0[ni].c2 = a + (b-a)*rng.nextDouble();
        break;
      default:
        err.println("mutate exception, parameter k= "+k+" non-existent");
        System.exit(1);
    }
    c0[ni].setFitness();

      //merge rank, return new pop. size
    return merge(-1,ni,af,z);
  }
//===
  int //return rank for reaping
  gauntlet(double livprob,int z)
/*::::::
  "Russian" roulette scheme, start at the bottom rank and work up.
  Independent of max/min sense
::::::*/
  {
  int kr;
 
      //Establish bug at risk, start at bottom of role
    for( kr = z-1; ; )
    {
      if( rng.nextDouble() > livprob) break; //livprob= prob. of surviving
      if( kr == 0) 
      {
          //cycled, tap random from bottom 2/3
        kr = (int)(0.66*z*rng.nextDouble()); //first top 2/3
        break;
      }
      kr--;
    }
   
    return kr;
  }
//===
  int  //return new population size
  removal(int kr, double[] af, int cs)
/*:::::
  sense independent, works for maximization or minimization, uses
  rank according to the role.
  kr;  rank for reaping
  cs;  colony size
:::::*/
  {
  int ki,i;
  int z=cs;

    ki = role[kr];  // killed bug index, fitness rank kr

      //remove bug ki = role[kr]
    af[0] *= z;  /* this is total fitness */
    --z;
      //following should work even if fitness is integer
    af[0] = (af[0] - c0[ki].fitness)/z;
      //must re-merge
    c0[ki].fitness = BUGINVALID; //invalid mark for index ki
    for( i=kr; i< z; ++i )  role[i] = role[i+1];

    return z;
  }
//mmm===
  public static void
  main(String[] args)
  {
 
    System.out.println("\nUsing:\n  "+
       "S0= "+dfS0+", X= "+dfX+", T(days)= "+dfnDays+
       ", rfr= "+dfrfr+", volatility= "+dfvol);

      //dtINdays is transmitted via static parameter
    amerputGAweb apga = new amerputGAweb(dfS0,dfrfr,dfvol,dfX,dfnDays);

  }
}
//kkk===
class critterAP //american put critter
{
  double a1,b1,c1;
  double a2,b2,c2;
  double fitness=0;
  double variance=0;
  amerputGAweb apg=null; //so this class can calculate the fitness

//ccc===
  public
  critterAP(amerputGAweb s)
  {
    apg = s;
  }
//===
  public
  critterAP(amerputGAweb s,
      double a1, double b1, double c1,
      double a2, double b2, double c2
    )
  {
    apg = s;
    this.a1 = a1; this.b1 = b1; this.c1 = c1;
    this.a2 = a2; this.b2 = b2; this.c2 = c2;
    setFitness();
  }
//===
  void
  setFitness()
  {
    double[] retvalue = apg.objective(a1,b1,c1,a2,b2,c2);
    fitness = retvalue[0]; variance = retvalue[1];
  }
//===
  void
  copy(critterAP clone)
  {
    clone.a1 = a1; clone.b1 = b1; clone.c1 = c1;
    clone.a2 = a2; clone.b2 = b2; clone.c2 = c2;
    clone.fitness = fitness; clone.apg = apg;  
    clone.variance = variance;
  }
}