/********************************************************************************
 * Project: FUSS                                                                *
 * Module : FussDD.c                                                            *
 * Author : Marcus Hutter                                                       *
 * Content: Fittness uniform selection on D-dimensional functions               *
 * Source : http://www.hutter1.de/ai/fuss.htm                                   *
 * (c) 2000 by Marcus Hutter                                                    *
 ********************************************************************************/

#include "FussDD.h"

/********************************************************************************/
/*                        D e s c r i p t i o n                                 */
/********************************************************************************

This file (together with FussDD.h) contains a compact and easy to understand 
(but somewhat dirty) implementation of FUSS for discretized fitness values. 

In evolutionary algorithms, the fitness of a population increases with time by
mutating and recombining individuals and by a biased selection of more fit
individuals. The right selection pressure is critical in ensuring sufficient
optimization progress on the one hand and in preserving genetic diversity to be
able to escape from local optima on the other. The code contains a new
selection scheme, which is uniform in the fitness values. It generates
selection pressure towards sparsely populated fitness regions, not necessarily
towards higher fitness, as is the case for all other selection schemes. This
selection scheme can be much more effective than standard selection schemes on
hard optimization problems.

See http://www.hutter1.de/ai/fuss.htm for a more detailed description: 
"Fitness Uniform Selection to Preserve Genetic Diversity" 
Proceedings of the 2002 Congress on Evolutionary Computation (CEC-2002) 783-788

Other versions:
FussTSP.cpp/h contains FUSS for the Traveling Salesman Problem with
              random distance and random Euclidian examples and 2.5Opt mutation

 ********************************************************************************/
/*                              F u n c t i o n s                               */
/********************************************************************************

- CStat::CStat()        Initialize statistics
- CInv::CInv()          Create D-dimensional function/individual
- CInv::Fitness()       Compute fitness of individual = function value
- CInv::Mutate()        Mutate coordinates
- CrossOver()           Recombination of two invs (one-point cross-over)
- CPop::CPop()          Create empty population (zero individuals)
- CPop::AddInv()        Add individual with possible replacement if arrays are full
- CPop::FussSelect()    Fitness uniform selection
- CPop::TournSelect()   Tournament or linear ranking selection with ts individuals
- CPop::SelectInv()     Select individual
- GAInit()              Initialize evolutionary algorithm (create population)
- GAExit()              Exitialize evolutionary algorithm (delete population)
- GAMain()              Evolutionary algorithm (Init,Run,Exit)


/********************************************************************************/
/*                           S t a t i s t i c s                                */
/********************************************************************************/

/*------------------------------------------------------------------------------*
  Description:  Initialize statistics
 *-MH---------------------------------------------------------------------------*/
CStat::CStat()
/*------------------------------------------------------------------------------*/
{
  nselect = nmutate = nadd = 0.;
  nruns   = 0;
} // End of CInv::CInv()


/********************************************************************************/
/*                   I n d i v i d u a l   T S P   P a t h s                    */
/********************************************************************************/

/*------------------------------------------------------------------------------*
  Description:  Create random D-dimensional function/individual
 *-MH---------------------------------------------------------------------------*/
CInv::CInv(int mode)
/*------------------------------------------------------------------------------*/
{
  int   i;
  for(i=0;i<DINV;i++) v[i]=rrand();
} // End of CInv::CInv()

/*------------------------------------------------------------------------------*
  Description:  Create D-dimensional function from other individual (copy)
 *-MH---------------------------------------------------------------------------*/
CInv::CInv(const CInv &source)
/*------------------------------------------------------------------------------*/
{
  int   i;
  for(i=0;i<DINV;i++) v[i]=source.v[i];
  fitness=source.fitness;
} // End of CInv::CInv()

/*------------------------------------------------------------------------------*
  Description:  Compute fitness of individual = function value
  Remark:       set switch(#) to desired function or add ne ones
 *-MH---------------------------------------------------------------------------*/
double CInv::Fitness()
/*------------------------------------------------------------------------------*/
{
  const double SE[] = { 0.3679, 0.3183, 0.5314, 0.1897}; // a,b,c, of simple example
  int i,e1,e2,e3;                    // feature flags
  const double SED=0.0001;                  // Delta of simple example
  fitness=0;                               // default fitness (should not occur)
  
  switch(4)
  {
    case 1: // Simple 3 valued example
    {
      e2=1; e3=1;
      for(i=0;i<DINV;i++)
      {
        e1=(v[i]>=SE[i] && v[i]<=SE[i]+SED);
        e3 &= e1; e2 &= !e1;
      }
      if      (e2) fitness=0.6;  // 2
      else if (e3) fitness=0.9;  // 3
      else         fitness=0.4;  // 1
      break;
    }
    case 2: // Simple continuous valued example
    { 
      double x2=(v[0]-SE[0])*(v[0]-SE[0]);
      double y2=(v[1]-SE[1])*(v[1]-SE[1]);
      double d2=SED;
      fitness=(x2/(x2+d2))*(y2/(y2+y2)/2 + exp(-(x2+y2)/d2));
      break;
    }
    case 3: // Simple 5 valued 3D example with 1->123, 2->4, 3->5
    {
      e1=(v[0]>=SE[0] && v[0]<=SE[0]+SED);
      e2=(v[1]>=SE[1] && v[1]<=SE[1]+SED);
      e3=(v[2]>=SE[2] && v[2]<=SE[2]+SED);
      if      (e1&&e2&&e3) fitness=0.9;  // 5
      else if (e1) fitness=0.6;  // 3
      else if (e2) fitness=0.4;  // 2
      else if (e3) fitness=0.2;  // 1
      else         fitness=0.75;  // 4
      break;
    }
    case 4: // Simple 4 valued 2D example with 1->12, 2->3, 3->4
    {
      e1=(v[0]>=SE[0] && v[0]<=SE[0]+SED);
      e2=(v[1]>=SE[1] && v[1]<=SE[1]+SED);
      if      (e1&&e2) fitness=0.9;  // 4
      else if (e1) fitness=0.5;  // 2
      else if (e2) fitness=0.3;  // 1
      else         fitness=0.7;  // 3
      break;
    }
    case 5: // Simple d+2 valued dD example (not yet functional)
    {
      int E[DINV];
      e1=1;
      for(i=0;i<DINV;i++)
      {
        E[i]=(v[i]>=SE[i] && v[i]<=SE[i]+SED);
        e1 &= E[i];
      }
      if      (e2) fitness=0.6;  // 2
      else if (e3) fitness=0.9;  // 3
      else         fitness=0.4;  // 1
      break;
    }
    default: printf("No valid Fitness() function");
  }
  return fitness;
} // End of CInv::Fitness()

/*------------------------------------------------------------------------------*
  Description:  Mutate coordinates
  Return Value: 0
  Remark:       Maybe more intelligent mutation useful
 *-MH---------------------------------------------------------------------------*/
double CInv::Mutate()
/*------------------------------------------------------------------------------*/
{
  v[irand(DINV)]=rrand();                  // one new random coordinate.
  /* ... */                             // bitflip
  /* ... */                             // local move
//  v[0]=rrand();  v[1]=rrand();          // two new random coordinates.
  stat.nmutate++;
//  fitness=Fitness();
  return 0.;
} // End of CInv::Mutate()

/*------------------------------------------------------------------------------*
  Description:  Recombination of two invs (one-point cross-over)
  Parameters:   2 invs
  Return Value: 0
 *-MH---------------------------------------------------------------------------*/
double CrossOver(CInv *inv1, CInv *inv2)
/*------------------------------------------------------------------------------*/
{
  int cp=irand(DINV-1);
  for(int i=0;i<=cp;i++)
  {
    double v=inv1->v[i]; inv1->v[i]=inv2->v[i]; inv2->v[i]=v;
  }
  return 0.;
} // End of CrossOver()

/********************************************************************************/
/*              P o p u l a t i o n   o f   I n d i v i d u a l s               */
/********************************************************************************/

/*------------------------------------------------------------------------------*
  Description:	Create empty population (zero individuals)
  Remark:       Use AddInv() to create non-empty populations
 *-MH---------------------------------------------------------------------------*/
CPop::CPop()
/*------------------------------------------------------------------------------*/
{
  int   i;
  ninvs = 0;
  fmax  = -1;
  fmin  = NFITLVL;
  maxnipf = NIPF;
  bestinv = NULL;
  for(i=0;i<NFITLVL;i++) nipf[i]=0;
} // End of CPop::CPop()


/*------------------------------------------------------------------------------*
  Description:  Add individual with possible replacement if arrays are full
  Parameters:   Individual inv
  Return Value: 1=Ok, 0=failure
  Remark:       Replacement currently only implemented for too many individual 
                per fitness but not if too many individuals overall
 *-MH---------------------------------------------------------------------------*/
int CPop::AddInv(CInv *inv)
/*------------------------------------------------------------------------------*/
{
  double fd = inv->Fitness();
  int f = (int)(NFITLVL*fd);

  if (f<0||f>=NFITLVL) { printf("fitness out of range in CPop::AddInv()"); return 0; }

  if (ninvs>=NINVS) 
  {
    printf("Max. number of individuals reached, deletion not yet implemented");
    return 0;
    //DecMaxNipf();
  }
  if (nipf[f]>=maxnipf)
  {
     int i=irand(nipf[f]);
//     if (fits[f][i]->Fitness() <= inv->Fitness()) // Test
     *fits[f][i] = *inv; // replace a random existing inv with same fitness
//     else printf("x");
  }else
  {
    fits[f][nipf[f]++] = &invs[ninvs];
    invs[ninvs++]= *inv;
    if (f>fmax) fmax=f;
    if (f<fmin) fmin=f;
  }
  if (bestinv && fd>bestinv->Fitness()) bestinv=inv;
  stat.nadd++;
  return 1;
} // End of CPop::AddInv()

/*------------------------------------------------------------------------------*
  Description:  Fitness uniform selection
  Return Value: pointer to individual
  Remark:    	Scale independent selection and FUSS+Hillclimbing may be activated
 *-MH---------------------------------------------------------------------------*/
CInv* CPop::FussSelect()
/*------------------------------------------------------------------------------*/
{
  int f=irand(fmax-fmin+1)+fmin; // FUSS-Selection
//  int f=fmax-(int)pow(fmax-fmin+1,rrand()); // SIF-Selection
//  int f=irand(2*(fmax-fmin))+fmin; f=min(fmax,f); // 50%FUSS+50%HILL
  
  while(nipf[f]==0) f++;
  //do { f=irand(fmax-fmin+1)+fmin; } while(nipf[f]==0);

  return fits[f][irand(nipf[f])];
} // End of CPop::FussSelect()

/*------------------------------------------------------------------------------*
  Description:  Tournament selection with ts individuals (=linear ranking for ts=2)
  Parameters:   p=0.5...1=probability of selecting better individual for ts=2
  Return Value: pointer to individual
 *-MH---------------------------------------------------------------------------*/
CInv* CPop::TournSelect(int ts, double p)
/*------------------------------------------------------------------------------*/
{
  int i1=irand(ninvs);

  if (ts==2)
  {
    int i2=irand(ninvs);
    int b1=invs[i1].Fitness()>invs[i2].Fitness();
    int b2=rrand()>=p;
    //if ( (invs[i1].Fitness()>invs[i2].Fitness()) ^ (rrand()>=p) )
    if (b1^b2) return &invs[i1]; else return &invs[i2];
  }else
  {
    for(int i=1;i<ts;i++)
    {
      int i2=irand(ninvs);
      if (invs[i1].Fitness()<invs[i2].Fitness()) i1=i2;
    }
    return &invs[i1];
  }
} // End of CPop::FussSelect()


/*------------------------------------------------------------------------------*
  Description:  Select individual
  Parameters:   mode=1=Fuss,  mode=2=Tournament
  Return Value: pointer to individual
 *-MH---------------------------------------------------------------------------*/
CInv* CPop::SelectInv(int mode)
/*------------------------------------------------------------------------------*/
{
  stat.nselect++;
  switch(mode)
  {
    case 1: // Fuss
      return FussSelect();
      break;
    case 2: // Tournamont
      return TournSelect(TOURNSIZE,T2SSTRENGTH);
      break;
  }
  return 0;
} // End of CPop::SelectInv()


/********************************************************************************/
/*                 E v o l u t i o n a r y   A l g o r i t h m                  */
/********************************************************************************/

/*------------------------------------------------------------------------------*
  Description:  Initialize evolutionary algorithm (create population)
  Return Value: 1
 *-MH---------------------------------------------------------------------------*/
int GAInit(void)
/*------------------------------------------------------------------------------*/
{
  int nsinv;
  pop = new CPop();
  if (SELECTMODE==1) nsinv=1; else nsinv=NIPF*NFITLVL;    // # of initial individuals
  for(int i=0;i<nsinv;i++) pop->AddInv(new CInv(1));
  return 1;                
} // End of GAInit()

/*------------------------------------------------------------------------------*
  Description:  Exitialize evolutionary algorithm (delete population)
  Return Value: 1
 *-MH---------------------------------------------------------------------------*/
int GAExit(void)
/*------------------------------------------------------------------------------*/
{
  delete pop; pop=NULL;
  return 1;                
} // End of GAExit()

/*------------------------------------------------------------------------------*
  Description:  Evolutionary algorithm (Init,Run,Exit)
  Return Value: 1
 *-MH---------------------------------------------------------------------------*/
int GAMain(void)
/*------------------------------------------------------------------------------*/
{
  GAInit();

  double f,flast=-1e99;
  int i;
  CInv inv,inv1,inv2;

//  printf("Fitness uniform Selection (Fuss) in GA for DD\n\n");
  for(i=0;/*i<NROUNDS*/;i++)
  {
    // copy and mutate
    inv= *(pop->SelectInv(SELECTMODE));
    inv.Mutate();
    pop->AddInv(&inv);
    // copy and crossover
    if (rrand()<CROSSPROB)
    {
      inv1= *(pop->SelectInv(SELECTMODE));
      inv2= *(pop->SelectInv(SELECTMODE));
      CrossOver(&inv1,&inv2);
      pop->AddInv(&inv1);
      pop->AddInv(&inv2);
    }
//    double lmin2=1e99;
//    for(k=0;k<pop->nipf[pop->fmax];k++)
    f=inv.Fitness();
    if (f>flast) 
    { flast=f; 
//      printf("R%-6d, BestFit=%lf, ninv=%d\n",i,f,pop->ninvs);
    }
    if (f>0.89) break;
  }
  stat.nruns++;

  GAExit();
  return 1;                
} // End of GAMain()

/*------------------------------------------------------------------------------*
  Description:  Start evolutionary algorithm
  Return Value: 1
 *-MH---------------------------------------------------------------------------*/
int main(void)
/*------------------------------------------------------------------------------*/
{
  for(;;) 
  {
    GAMain();
    printf("\nRun%d, AvSelect=%lf, AvMutate=%lf, AvAdd=%lf\n",
    stat.nruns,stat.nselect/stat.nruns,stat.nmutate/stat.nruns,stat.nadd/stat.nruns);
  }
  return 1;                
} // End of main()

/*--------------------------End-of-File-FussDD.cpp------------------------------*/