/*
 * gen.c
 * Exemplo simples de algoritmo genetico com gene real
 * Autor: Sadao Massago
 *   http://www.dm.ufscar.br/~sadao
 * Data: 2008-10-20
 * Licenca: GNU GPL v2 ou posterior
 * http://www.gnu.org/licenses/gpl-2.0.html
 * http://www.gnu.org/copyleft/gpl.html
 * */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <assert.h>
#include <time.h>

/* tamanho da populacao */
#define POP_SIZE 4

/* numero de genes no cromossomo */
#define N 2

/* intervalo para individuo inicial eh [A0, B0] */
#define A0 -3.0
#define B0 3.0

/* largura da faixa para busca (perturbacao) na mutacao eh 2*L
   (L para cima e para baixo) */
#define L 0.5

/* numero de mutacao */
#define NUM_MUTACAO 1
/* numero de cruzamentos */
#define NUM_CRUZAMENTO 1

/* somente para steady-state */
/* numero de novos */
#define NUM_NOVOS_STEADY_STATE 1

/* somente para generacional */
/* numero de elite (sobrevivente) >= 1 */
#define NUM_ELITE 1

/* Necessario que
   No generacional: NUM_MUT + NUM_CR + 1 <= POP_SIZE
   No steady-state: NUM_MUT + NUM_CR  <= POP_SIZE  */

/* numero de iteracao para teste (no programa principal) */
#define NUM_ITER 10

/* De quanto em quanto (iteracao0 efetua a limpeza */
#define NUM_FREQ_LIMPAR 5

/* funcao objetiva que pretende minimizar */
/* por enquanto, eh funcao quadratica (com minimo na origem) */
double f(double x[])
{
    double soma = 0.0;
    int i;
    for(i=0; i<N; i++)
      soma += x[i]*x[i];
    return soma;
}


/* numero aleatorio */
/* retorna em [0, 1] */
double drand()
{
    return ( (double)rand() ) / ( (double) RAND_MAX  );
}

/* individuo */
/* cromossomo eh apenas um numero real. */
typedef struct {
  double x[N];
  double custo;
} ind;

/* Comparacao de custo
  Retorna valor inteiro:
  custo de a < custo de b => negativo
  custo de a == custo de b => 0
  custo de a > custo de b => positivo

  Para usar o sort() da biblioteca, recebe void *.
  Tera que converter para ind * antes de acesar elementos
*/

/*  a eh melhor que b deve retornar -1 (a vem antes de b) */
int compara(const void* a, const void * b)
{
  if( ((ind *)a)->custo < ((ind *)b)->custo) return -1;
  else if (((ind *)a)->custo ==  ((ind *)b)->custo) return 0;
  else return 1;
}
/* comparacao acima pode ser usado na ordenacao ou busca da biblioteca padrao
void qsort (void *array, size_t count, size_t size, comparison_fn_t compare)
void * bsearch (const void *key, const void *array, size_t count, size_t size, comparison_fn_t compare)
*/

/* comparacao de elementos: elementos sao iguais? 
 * Tambem pode efetaur comparacao do tipo proximos? */

int igual (const ind *item1, const ind *item2)
{
  int i;
  
  for(i=0; i<N; i++)
    /* aproximadamente igual */
    if(fabs(item1->x[i] - item2->x[i]) > 1E-10) 
      return 0; /* false */
  return 1; /* true */
}

/* inicializa */
void inicia(ind *item)
{
    int i;
    for(i=0; i<N; i++)
      item->x[i] = A0 + (B0 - A0) * drand();
    item->custo = f(item->x);
}

/* mutacao */
/* perturbacao uniforme em um ponto */
/* largura da mutacao eh no maximo L para cima e L para baixo */
/* espaco de busca eh vizinhanca retangular com centro no
individuo e largura L */

void mutacao(ind *item)
{
    int i;
    /* nao eh recomendado usar % para este caso, mas ... */
    i = rand() % N;
    item->x[i] += 2 * L * (drand()-0.5);
    item->custo = f(item->x);
}


/* cruzamento */
/* Metodo "Blending": 
 * a media ponderada de cada gene dos cromossomos
   peso muda de gene por gene e tem como espaco de busca no
   retangulo N dimensional com vertices nos pai1 e pai 2
   Nota: Se o peso for mesmo para todo gene, espaco de busca
   sera segmento de reta que liga pai1 e pai2 (bem reduzido)
 */
void cruzamento(ind *item, const ind *pai1, const ind *pai2)
{
    int i;
    double t;
    for(i=0; i<N; i++)
    {
      t = drand();
      item->x[i] = t * (pai1->x[i]) + (1.0 - t) * (pai2->x[i]);
    }
    item->custo = f(item->x);
}

/* genetico */

/* inicializa a populacao */
void inicia_pop(ind *pop)
{
  int i;
  for(i=0; i<POP_SIZE; i++)
    inicia(& (pop[i]) );
  qsort (pop, POP_SIZE, sizeof(ind), compara);
}

/*usado somente pelo steady-state */
/* insere, so se for melhor. */
void insere(const ind *item, ind *pop)
{
        /* se melhor, entra na populacao */
        /* nota: bsearch retorna null se nao encontrado
              como queremos inserir no caso em que nao existe, 
              nao podemos usar o bsearch
        pos = (ind *)bsearch (&item, pop, N, sizeof(ind), compara);
        */

 int i, j;

 /* busca linear */
 for(i=0; i<POP_SIZE; i++)
   /* se achou posicao (melhor que alguem), pare */
   if( compara(item, &(pop[i]) )  <= 0)
       break;
        
 /* Se achou posicao e nao existe ainda */
 if( (i < POP_SIZE) && !igual(item, &(pop[i]) ) )
 {
     /* desloca para frente, puxando os itens de tras para frente */
     for(j=POP_SIZE-1; j>i; j--)
       pop[j] = pop[j-1];
     /* insere */
     pop[i] = *item;
 }
}

/****************************************/
/* steady-state: so entra se for melhor */
/****************************************/
/* efetua uma iteracao do steady-state */
void steady_state(ind *pop)
{
    int i, j1, j2;
    ind item;

    /* muta NUM_MUTACAO elementos */
    for(i=0; i<NUM_MUTACAO; i++) {
        /* seleciona um individuo */
        /* nao eh recomendado usar % para este caso, mas ... */
        j1 = rand() % POP_SIZE; 
        /*printf("DEBUG: steady-sate:mutacao: j1 = %d\n", j1);*/
        item = pop[j1];
        /* muta */
        mutacao(&item);
        /* insere se eh melhor: teste eh dentro do insere */
        insere(&item, pop);
    }

    /* cruzamento */
    for(i=0; i<NUM_CRUZAMENTO; i++)
    {
        /* nao eh recomendado usar % para este caso, mas ... */
        j1 = rand() % POP_SIZE;
        /* j2 nao deve ser igual a j1 */
        /* eliminando j1, tera POP_SIZE-1 elementos */
        /* nao eh recomendado usar % para este caso, mas ... */
        j2 = rand() % (POP_SIZE-1);  
        /* se for depois de j1, posicao real de j2 eh um para frente */
        if(j1 <= j2) j2++;
        cruzamento( &item, &(pop[j1]), &(pop[j2]));
        /* insere se eh melhor: teste eh dentro do insere */
        insere(&item, pop); 
    }

    /* novos */
    for(i=0; i<NUM_NOVOS_STEADY_STATE; i++) {
        inicia(&item);
        /* insere se eh melhor: teste eh dentro do insere */
        insere(&item, pop);
    }
} /* steady_state() */

/****************************************/
/* Usado somente pelo generacional:
 * troca a populacao */
/****************************************/

void troca_pop(ind *pop1, ind *pop2)
{
    ind temp;
    int i;

    for(i=0; i<POP_SIZE; i++)
    {
      temp = pop1[i];
      pop1[i] = pop2[i];
      pop2[i] = temp;
    }
}

/***************************************
 * generacional: renova toda populacao
 * *************************************/
/* efetua uma iteracao do generacional */
void generacional(ind *pop)
{
    int count, i, j1, j2;
    ind pop_ant[POP_SIZE];

    troca_pop(pop, pop_ant);

    /* elite sempre sobrevive */
    for(count=0; count<NUM_ELITE; count++)
        pop[count] = pop_ant[count];

    /* muta NUM_MUT elementos seguintes */
    for(i=0; i<NUM_MUTACAO; count++,i++) {
        /* seleciona um individuo */
        /* nao eh recomendado usar % para este caso, mas ... */
        j1 = rand() % POP_SIZE;  
        pop[count] = pop_ant[j1];
        /* muta */
        mutacao( &(pop[count])  );
    }

    /* cruzamento */
    for(i=0; i<NUM_CRUZAMENTO; count++,i++)
    {
        /* nao eh recomendado usar % para este caso, mas ... */
        j1 = rand() % POP_SIZE;  
        /* j2 nao deve ser igual a j1 */
        /* eliminando j1, tera POP_SIZE-1 elementos */
        /* nao eh recomendado usar % para este caso, mas ... */
        j2 = rand() % (POP_SIZE-1);
        /* se for depois de j1, posicao real de j2 eh um para frente */
        if(j1 <= j2) j2++;
        cruzamento( &(pop[count]), &(pop_ant[j1]), &(pop_ant[j2]));
    }

    /* restante eh preenchido pelo individuo novo */
    for(;count<POP_SIZE; count++)
       inicia(&(pop[count]) );
    /* reordena */
    qsort (pop, POP_SIZE, sizeof(ind), compara);
} /* generacional() */


/* programa de teste */

void imp_pop(const ind *pop)
{
    int i, j;

    for(j=0; j<N; j++)
      printf("x[%d]\t", j);
    printf(" f(x)\n");

    for(i=0; i<POP_SIZE; i++)
    {
      for(j=0; j<N; j++)
        printf("%f\t", pop[i].x[j]);
      printf(" %f\n", pop[i].custo);
    }
}

/* Como pode aparecer elementos repetidos, precisa de limpeza
 * Generacional costuma requerer maior frequencia */
void limpar(ind *pop)
{
    int i, j;
    for(i=0; i<POP_SIZE; i++)
    {
      /* tem alguem na frente que eh igual ? */
      for(j=i+1; j<POP_SIZE; j++)

      if( igual(& (pop[i]), & (pop[j]) ) ) /* se tem, gera um novo */
        inicia(&pop[j]);
      else /* se nao tem mais elemento igual a i-esimo */
        break; /* sai do laco de j (retorna ao laco de i) */
    }
    /* reordena */
    qsort (pop, POP_SIZE, sizeof(ind), compara);
}

int main(int argc, char *argv[])
{
  ind pop[POP_SIZE];

  int i;

  /*  inicializa com valor aleatorio do relogio */
  /*  srand(time(NULL)); */

  /* steady state */
  
  printf("========> Steady state <======\n");
  inicia_pop(pop);
  printf("Populacao inicial\n");
  imp_pop(pop);

  for(i=0; i<NUM_ITER; i++)
  {
     steady_state(pop);
     printf("i = %d\n", i);
     imp_pop(pop);
     if(((i+1) % NUM_FREQ_LIMPAR) == 0)
     {
        /* eliminar elementos repetidos */
        limpar(pop);
        printf("limpza efetuada\n");
        imp_pop(pop);
   }
  }


  /* generational */
  /* deve ter (NUM_ELITE+NUM_MUTACAO+NUM_CRUZAMENTO) <= POP_SIZE
     eh recomendavel que  (NUM_ELITE+NUM_MUTACAO+NUM_CRUZAMENTO) < POP_SIZE
       para ter um novo entrando, se tiver muito minimo local */
  assert( (NUM_ELITE + NUM_MUTACAO + NUM_CRUZAMENTO) <=  POP_SIZE );

  printf("\n");
  printf("==========> Generacional <==========\n");
  inicia_pop(pop);
  printf("Populacao inicial\n");
  imp_pop(pop);

  for(i=0; i<NUM_ITER; i++)
  {
     generacional(pop);
     printf("i = %d\n", i);
     imp_pop(pop);
     if(((i+1) % NUM_FREQ_LIMPAR) == 0)
     {
        /* eliminar elementos repetidos */
        limpar(pop);
        printf("limpza efetuada\n");
        imp_pop(pop);
   }
  }

  return 0; /* sucesso */
} /* main() */