/*
Primer classification genetic algorithm.  Using finite state machines 
to divide primers into two classes.

Modified to perform hybridization.
Modified to use two classifications.
*/

#include <iostream.h>
#include <fstream.h>
#include <string.h>
#include <ctype.h>
#include <stdlib.h>
#include <stdio.h>
#include <time.h>
#include <math.h>

#include "vaut.h"

//data file for training data
#define TrainingData "primer2000.dat"

//population size
#define Popsize 600

//number of generations
#define gen 1000

//number of runs
#define runs 100

//maximum number of training examples
#define texm 4200

//number of states per automaton
#define NS 64

//number of goodies
#define hybrids 100

//payoff[i][j] for automata of type i, state j 
int pay[2][3]={{ 1,-1, 0},  //pay[good,bad][good,bad,duh?]  
               {-1, 1, 0}
              };

int nt;             //number of training examples
aut pop[Popsize];   //population of variable automata
int fit[Popsize];   //fitness of population members
int dx[Popsize];    //sorting index
int tp[texm];       //primer type
int *pr[texm];      //primer 
int pl[texm];       //primer length
int bdex;           //index of best classifier

int CGAT(char c);           //numerate a character
void ReadData(char *fn);    //read in a data set with NS sequences of length SL
int fitness(aut &a);        //compute the quality of an arrangement
void initialize(int run);   //initialize population
void mutate(aut &a);        //glue together the vaut primitive mutations
void breed();               //breed configurations
void report(ostream &aus);  //report results
void classification(aut &a, //run the classifier on the inputs
                    ostream &aus
		    );

main(){

int g,r;
fstream best,aus;
char fn[60];

  srand48(121461);
  ReadData(TrainingData);
  best.open("classAut.dat",ios::out);
  for(r=0;r<runs;r++){
    initialize(r);
    sprintf(fn,"run%d.dat",r);
    aus.open(fn,ios::out);
    for(g=0;g<gen;g++){//generation loop
      cout << r << " " << g << " ";
      breed();
      report(aus);
      classification(pop[bdex],cout);
    }
    aus.close();
    best << pop[bdex];
    classification(pop[bdex],best);
    best << fit[bdex] << endl;
  }
  best.close();
}

void ReadData(char *fn){

fstream inp;
char buf[300];
int asmb[50];
int i,j,q;
char tr[5]="CGAT";

  inp.open(fn,ios::in);
  nt=0;
  while(!inp.eof()){
    inp.getline(buf,299);
    if(strlen(buf)>0){
      tp[nt]=atoi(buf);
      j=0;
      for(i=0;i<strlen(buf);i++){
	        q=CGAT(buf[i]);
        if(q<4)asmb[j++]=q;      
      }
      pl[nt]=j;
      pr[nt]=new int[j];
      for(i=0;i<j;i++)pr[nt][i]=asmb[i];
      nt++;
    }
  }
  inp.close();
  cout << "found " << nt << " training examples" << endl;
}

void initialize(int run){

int i,j;
fstream hyb;
char buf[256];

  if(run==0)for(i=0;i<Popsize;i++){
    pop[i].create(NS);
    dx[i]=i; 
  } else for(i=0;i<Popsize;i++)pop[i].recreate();

  if(hybrids>0){
    hyb.open("BestClA.dat",ios::in);
    for(i=0;i<hybrids;i++){
      hyb >> pop[i];
      //cout << pop[i];
      for(j=0;j<3;j++)hyb.getline(buf,255);
    }
    hyb.close();
  }

  for(i=0;i<Popsize;i++){
    fit[i]=fitness(pop[i]);
    //cerr << fit[i] << " ";
  }
  //cout << endl;
}

int fitness(aut &a){//Compute the fitness

int i,j,ttl,dum;

  //cerr << "Fit eval" << endl;
  ttl=0;

  for(i=0;i<nt;i++){
    a.reset();
    //for(j=0;j<pl[i];j++)dum=a.run(pr[i][j]); //burn in supressed
    for(j=0;j<pl[i];j++){
      ttl+=pay[tp[i]][a.run(pr[i][j])];
    }
  }

  return(ttl);
  
}

int CGAT(char c){//numerate a character

  switch(c){
  case 'C':
  case 'c': return(0);
            break;
  case 'G':
  case 'g': return(1);
            break;
  case 'A':
  case 'a': return(2);
            break;
  case 'T':
  case 't': return(3);
            break;
  }
  //cout << "Alert!!!!!" << endl;
  return(4);
}

void swap(int i,int j){

int sw;

   sw=dx[i];
   dx[i]=dx[j];
   dx[j]=sw;

}

int better(int i,int j){

   if(fit[i]>fit[j])return(1);
   return(0);
}

void shuffle4(int *fit,int * dx,int popz){

int i,rp;

   for(i=0;i<Popsize;i++){  //randomize 
      rp=lrand48()%Popsize;
      swap(i,rp);
   }
   for(i=0;i<Popsize;i+=4){  //sort into groups of four
      if (better(dx[i],dx[i+3])  ) swap(i,i+3);
      if (better(dx[i+1],dx[i+2])) swap(i+1,i+2);
      if (better(dx[i],dx[i+1])  ) swap(i,i+1);
      if (better(dx[i+2],dx[i+3])) swap(i+2,i+3);
      if (better(dx[i+1],dx[i+2])) swap(i+1,i+2);
   }
}

void mutate(aut &a){

  switch(lrand48()%10){
  case 0: a.initialmutation();
          break;
  case 1: 
  case 2:
  case 3: a.mutatetransition();
          break;
  case 4: 
  case 5:
  case 6: 
  case 7:
  case 8: 
  case 9: a.mutateaction();
          break;

  }

}

void breed(){

int i,j,cp;
float g[5],h[5],sw;

  shuffle4(fit,dx,Popsize);
  for(i=0;i<Popsize;i+=4){//for each tournament
    copy(pop[dx[i]],pop[dx[i+2]]);    
    copy(pop[dx[i+1]],pop[dx[i+3]]);
    tpc(pop[dx[i]],pop[dx[i+1]]);
    mutate(pop[dx[i]]);
    mutate(pop[dx[i+1]]);
    fit[dx[i]]=fitness(pop[dx[i]]);
    fit[dx[i+1]]=fitness(pop[dx[i+1]]);
  }
}

void report(ostream &aus){

int i,best;
double mu,sg,val;

  best=0;
  mu=0.0;
  sg=0.0;
  for(i=0;i<Popsize;i++){
    val=((double)fit[i]);
    mu+=val;
    sg+=val*val;
    if(best<fit[i]){
      best=fit[i];
      bdex=i;
    }
  }
  mu/=Popsize;
  sg=sg/Popsize-mu*mu;
  if(sg<0)sg=0; else sg=sqrt(sg);
  cout << mu << " " << sg << " " << best << endl;
  aus << mu << " " << sg << " " << best << endl;

}

void classification(aut &a,
                    ostream &aus
                   ){

int i,j,pt;
int wh[3];
int C[2][2];
char Q[50];

  for(i=0;i<2;i++)for(j=0;j<2;j++)C[i][j]=0;
  for(i=0;i<nt;i++){
    for(j=0;j<3;j++)wh[j]=0;
    a.reset();
    for(j=0;j<pl[i];j++)wh[a.run(pr[i][j])]++;
    pt=0;
    for(j=1;j<3;j++)if(wh[j]>wh[pt])pt=j;
    C[tp[i]][pt]++; 
  }
  for(i=0;i<2;i++){
    for(j=0;j<2;j++){
      sprintf(Q,"%5d ",C[i][j]);
      aus << Q;
    }
    aus << endl;
  }
}