/********************************************************************************
 * Project: FUSS                                                                *
 * Module : FussTSP.c                                                           *
 * Author : Marcus Hutter                                                       *
 * Content: Fittness uniform selection for TSP                                  *
 * Source : http://www.hutter1.de/ai/fuss.htm                                   *
 * Copyright (c) 2000 by Marcus Hutter                                          *
 ********************************************************************************/

#include "fussTSP.h"

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

This file (together with FussTSP.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:
FussDD.cpp/h  contains FUSS for D-dimensional function optimization with
              the example as given in the paper.

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

- CCities::CCities()    Create cities and distances
- CInv::CInv()          Create TSP route from scratch or other individuals
- CInv::Fitness()       Compute fitness of TSP route from its length
- CInv::Flip2Links()    Insert city into (link1,link2) and remove it from (cprev,cnext),
- CInv::Mutate()        Mutate path (with 2.5Opt or 1Opt)
- 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
- InvTest()             Test consistency of incremental length update in when mutating individuals
- GAInit()              Initialize evolutionary algorithm (create population)
- GAExit()              Exitialize evolutionary algorithm (delete population)
- GAMain()              Evolutionary algorithm (Init,Run,Exit)

 ********************************************************************************/
/*         C i t i e s   w i t h   D i s t a n c i e s   f o r   T S P          */
/********************************************************************************/

/*------------------------------------------------------------------------------*
  Description:  Create cities and distances
  Parameters:   mode=1=random, mode=2=euklidian
 *-MH---------------------------------------------------------------------------*/
CCities::CCities(int mode)
/*------------------------------------------------------------------------------*/
{
  int   i,k;
  ncities = NCITIES;
  for(k=0;k<ncities;k++)
    for(i=k;i<ncities;i++)
      if (mode==1) dist[i][k]=dist[k][i]=rrand();
      else dist[i][k]=dist[k][i]=(abs((i%NCSQRT)-(k%NCSQRT))+abs((i/NCSQRT)-(k/NCSQRT)))/((float)NCITIES/2);
} // End of CCities::CCities()

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

/*------------------------------------------------------------------------------*
  Description:  Create TSP route from scratch
  Parameters:   mode=0=standard oder 0-1-2-... mode=1=random
  Remark:       To create random route, flip cities sufficiently often
 *-MH---------------------------------------------------------------------------*/
CInv::CInv(int mode)
/*------------------------------------------------------------------------------*/
{
  int   i;
  length=0.;
  for(i=1;i<cities->ncities;i++) { next[i-1]=i; length+=cities->dist[i-1][i]; }
  next[cities->ncities-1]=0; length+=cities->dist[cities->ncities-1][0];
  if (mode>0) for(i=0;i<40*cities->ncities;i++) Mutate();
} // End of CInv::CInv()

/*------------------------------------------------------------------------------*
  Description:  Create TSP route from other individual (copy)
 *-MH---------------------------------------------------------------------------*/
CInv::CInv(const CInv &source)
/*------------------------------------------------------------------------------*/
{
  int   i;
  for(i=0;i<cities->ncities;i++) next[i]=source.next[i];
  length=source.length;
} // End of CInv::CInv()

/*------------------------------------------------------------------------------*
  Description:  Compute fitness of TSP route from its length
  Return Value: 2*NFITLVL/length
 *-MH---------------------------------------------------------------------------*/
double CInv::Fitness()
/*------------------------------------------------------------------------------*/
{
  return 2*NFITLVL/length;
} // End of CInv::Fitness()

/*------------------------------------------------------------------------------*
  Description:  Insert city into (link1,link2) and remove it from (cprev,cnext),
                i.e. apply 2.5 Opt operator to TSP route
  Return Value: Increase in length of route 
 *-MH---------------------------------------------------------------------------*/
double CInv::Flip2Links(int cprev, int link1)
/*------------------------------------------------------------------------------*/
{
  int city =next[cprev];
  int cnext=next[city];
  int link2=next[link1];
  
  if (link1==city||link2==city) return 0.; // overlap
  next[cprev]=cnext;
  next[link1]=city;
  next[city] =link2;
  double delta = (cities->dist[cprev][cnext]+cities->dist[link1][city]+cities->dist[city][link2]) -
                 (cities->dist[cprev][city]+cities->dist[city][cnext]+cities->dist[link1][link2]);
  length += delta;
  return delta;
} // End of CInv::Flip2Links()

/*------------------------------------------------------------------------------*
  Description:  Mutate path (with 2.5Opt or 1Opt)
  Return Value: Increase in length of route 
  Remark:       2.5 Opt currently active, 1Opt may be activated (see code)
  Remark:       Maybe more intelligent multiple mutation useful
 *-MH---------------------------------------------------------------------------*/
double CInv::Mutate()
/*------------------------------------------------------------------------------*/
{
  // 2.5 Opt
  return Flip2Links(irand(NCITIES),irand(NCITIES));
  // 1 Opt
  int c=irand(NCITIES); // test of a more local mutation
  return Flip2Links(c,(c+2)%NCITIES);
} // End of CInv::Mutate()

/********************************************************************************/
/*              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)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;
  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 (good)
//  int f=irand(2*(fmax-fmin))+fmin; f=min(fmax,f); // 50%FUSS+50%HILL (bad)
  
  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)
/*------------------------------------------------------------------------------*/
{
  switch(mode)
  {
    case 1: // Fuss
      return FussSelect();
      break;
    case 2: // Tournamont
      return TournSelect(TOURNSIZE,T2SSTRENGTH);
      break;
  }
  return 0;
} // End of CPop::SelectInv()

/*------------------------------------------------------------------------------*
  Description:  Test consistency of incremental length update in when mutating individuals
  Return Value: 1=Ok, 0=error
 *-MH---------------------------------------------------------------------------*/
int InvTest(void)
/*------------------------------------------------------------------------------*/
{
  int i;
  CInv *inv = new CInv();

  for(i=0;i<100;i++) 
  {
    inv->Mutate();
//    printf("length=%lf\n",inv->length);
  }
  double l=0.;
  for(i=0;i<cities->ncities;i++) l+=cities->dist[i][inv->next[i]];

  if (fabs(l-inv->length)>1e-9) 
  { printf("Inconsistency in TSP length\n"); delete inv; return 0; }
  else 
  { printf("Incremental length update is perfect\n"); delete inv; return 1; }
} // End of InvTest()


/********************************************************************************/
/*                 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;
  cities = new CCities(CITYGEOMETY);
  InvTest();
  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 cities; cities=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;

  printf("Fitness uniform Selection (Fuss) in GA for TSP\n\n");
  for(i=0;/*i<NROUNDS*/;i++)
  {
    inv= *(pop->SelectInv(SELECTMODE));
    inv.Mutate();
    pop->AddInv(&inv);
//    double lmin2=1e99;
//    for(k=0;k<pop->nipf[pop->fmax];k++)
    f=inv.Fitness();
    if (f>flast) 
    { flast=f; 
      printf("R%-6d, BestLength=%lf, BestFit=%lf, ninv=%d\n",i,inv.length,f,pop->ninvs);
    }
  }

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

/*------------------------------------------------------------------------------*
  Description:  Start evolutionary algorithm
  Return Value: 1
 *-MH---------------------------------------------------------------------------*/
int main(void)
/*------------------------------------------------------------------------------*/
{
  GAMain();
  return 1;                
} // End of main()

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