// Un generateur d'echantillon de recombinaisons intercross...



//----------------- INCLUDE files standard ------------------------------------

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

  int tr[256][256];
//------------------- Les classes utilisees -----------------------------------
typedef char* gamete;

class marqueur
{
public:
    double position;
    double frecomb;
  
  marqueur()  { position = frecomb = 0.0;};
  ~marqueur() {};
};
// Un marqueur est defini par sa position sur le chromosome et la FR entre lui et son precedent 

class individu
{
public:
  gamete pere;
  gamete mere;
  
  individu() {pere=NULL; mere=NULL;};
  void init(int n) {pere=new char[n+1];memset(pere,'0',n);pere[n]=0;                               mere=new char[n+1];memset(mere,'0',n);mere[n]=0;};
  ~individu() {delete[] pere; delete[] mere;}; 
  void affiche (ofstream & out){out << pere <<"\t"<< mere <<"\t";}; 
};
// Chaque individu est cree a partir d'un gamete paternel et un maternel


//------------- Generateur aléatoire de Simon ---------------------------------

double randomax ()
{ 
    static int premiere_fois = 0;
    char *demo;
    long heure;

// Le code qui suit initialise le generateur de nombre aleatoire  a une           valeur precise pour obtenir la meme suite de nombres aleatoires si DEMO        est present dans l'environnement d'appel a l'execution. Mais cette liste       de nombre est identique pour chaque creation de processus, ce qui donne        un caractere tres repetitifs aux taches. Ainsi, en ne definissant pas DEMO,    on initialise le random avec une valeur pseudo-aleatoire qui est l'heure       calendaire multipliee par le pid du processus 

  if (!premiere_fois) 
      {
	  premiere_fois = 1;
	  if (!(demo=getenv("DEMO")))
	    {
	      heure = (long) time(NULL);
	      srand48(heure * getpid());
	    }
	  else
	    {
	      srand48(atol(demo));
	    }
      }
    return drand48();
} 


// ---------------- Comparer la position de deux marqueurs --------------------

int compmarq(const void* x, const void* y)
{
    marqueur* a = (marqueur*)x;
    marqueur* b = (marqueur*)y;

    
    if ((a->position) > (b->position))
	return 1;
    else if ((a->position) < (b->position))
	return -1;
    else 
        return 0;
}


// ---------------- Calcule la f. de recomb. a partir de la distance ----------

// Fonction de Haldane, Independance + Poisson 
double Haldane(double distance)
{
    return ((1-exp(-2*fabs(distance)))/2);
}

// Fonction de Kosambi. 
double Kosambi(double distance)
{
    return (0.5*tanh(2*fabs(distance)));
}


// ---------------- Générateur de marqueurs -----------------------------------

marqueur* genmarq(int nombre, double length)
{
    marqueur* Eux = new marqueur[nombre];

    int i;
    for (i=0; i<nombre; i++)  
       {
	Eux[i].position = randomax()*length;
       }
    
    qsort(Eux,nombre,sizeof(marqueur),compmarq);
    // Apres avoir cree les marqueurs en leur attribuant une position aleatoire       sur le chromosome, on les range selon leur ordre sur le chromosome

    // Puis on calcule la f. de recomb via Haldane
    double distance;
    for (i=1; i<nombre; i++)
      {
	distance = Eux[i-1].position-Eux[i].position;
	Eux[i].frecomb = Haldane(distance);
	// Eux[i].frecomb = Kosambi(distance);
      }
    
    return Eux;
    // On obtient en fait une carte genetique
}


// -------- Générateur de carte uniforme (d constante entre marq.succ.) -------

marqueur* cartunif(int nombre, double length)
{
    marqueur* Eux = new marqueur[nombre];

    int i;
    for (i=0; i<nombre; i++)  
       {
	Eux[i].position = i*length/(nombre-1);
       }
    
    // Les marqueurs sont espaces de maniere reguliere, dans l'ordre de numerotation usuel croissant le long du chromosome

    // Puis on calcule la f. de recomb via Haldane
    double distance;
    for (i=1; i<nombre; i++)
      {
	distance = Eux[i-1].position-Eux[i].position;
	Eux[i].frecomb = Haldane(distance);
	// Eux[i].frecomb = Kosambi(distance);
      }
    
    return Eux;
    // On obtient en fait une carte genetique
}


//----------- Generateur de gamète --------------------------------------------

gamete gengamete(int nombre, marqueur* modele)
{
    gamete parent;
    int parite;
    int i; 
	
          // On cree un gamete aleatoire d'apres le modele de carte genere 

          parent=new char[nombre+1];
          memset(parent,'0',nombre);parent[nombre]=0;
          if (randomax() < 0.5)
              {parent[0] ='1';}
          parite = parent[0]-'0';  
          for(i=1; i<nombre; i++)
	     {
	         if (randomax() < modele[i].frecomb)
		   parite = (1-parite);
		 parent[i] = '0'+parite;
                 // Si randomax < FR il y a eu recombinaison
	     }
          return parent;
}
 

//---------------- Generateur modele ------------------------------------------

individu* gensample(int taille, int nombre, marqueur* modele)
{
    individu* sample =new individu[taille];
    int i;

    for (i=0; i< taille; i++)
      {
       sample[i].pere=gengamete(nombre, modele);
       sample[i].mere=gengamete(nombre, modele);
       // Pour chaque individu de l'echantillon, on cree un gamete paternel et           un maternel 
     } 
    return sample;
}


//---------------- Generateur d'échantillon -----------------------------------

individu* incomplet(int taille, int nombre,individu *echantillon, double proba,                    double probae)
{
    individu* sample =new individu[taille];
    int i,j;

    for (i=0; i< taille; i++)
      {
       sample[i].init(nombre); 
       for (j=0; j<nombre; j++)
         {
          if (randomax()<proba)
	        {
	         if (randomax() < probae)
		   sample[i].pere[j]='1'- echantillon[i].pere[j];
                 else 
                   sample[i].pere[j]=echantillon[i].pere[j];
		}
          else  sample[i].pere[j]='-';
          if (randomax()<proba)
	        {
	         if (randomax() < probae)
		   sample[i].mere[j]='1'-echantillon[i].mere[j];
                 else 
                   sample[i].mere[j]=echantillon[i].mere[j];
		}
          else  sample[i].mere[j]='-';
         // Pour chaque individu, on introduit erreurs et manques
	  }
     } 
    return sample;
}


//--------------------- Dump du modele ----------------------------------------

void pontemodele(int nombre,marqueur* marqs,ofstream & out)
{
  int i;
  
  out << "Carte:\n";
  for (i=0; i<nombre; i++)
    {
      out << "M" << i+1  << " à " << 100*marqs[i].position << " cM";
      if (i) 
	{
	 out << "\tdistance : " << 100*(marqs[i].position-marqs[i-1].position)              <<  " cM";
	 out << "\tF.Rec.: " << marqs[i].frecomb;
	}
      out << endl;
    }
}


//-------------------- Codage intercross --------------------------------------

int code(individu* echant, int numero, int marq, int* col)
{
  int codage[6][3]={{5,1,1},
		    {5,4,3},
                    {2,4,2},
                    {5,5,5},
                    {4,4,4},
                    {0,0,0}};
  int colonne;
  int ligne;

  

  colonne=col[marq];
  ligne=tr[echant[numero].pere[marq]][echant[numero].mere[marq]];
  return codage[ligne][colonne];  
}


//------------------- Dump de l'echantillon -----------------------------------

void pontesample(int taille,int nombre,individu* echantillon,int* colonne,                      ofstream & out)
{
  int i,j;

  out << "data type f2 intercross\n";
  out << taille << " " << nombre << " 0 symbols 1=A 2=B 3=H 4=C 5=D 0=-\n";
    
   for (i=0; i<nombre; i++)
     {
       out << "*M" << i+1 << "\t";
       
       
       for (j=0;j<taille; j++)
	 out << code(echantillon,j,i,colonne);
       out << endl;
     }
}


//------------- Vecteur dominance ---------------------------------------------

int* dominance(int nombre,double codomin)
{
 int i;
 int* colonne=new int[nombre];
 
 for (i=0; i<nombre; i++)      // La dominance est propre a chaque marqueur 
   {
    if (randomax()<codomin)  // Les deux alleles sont codominants
          colonne[i]=2;
    else if (randomax()<0.5) 
	   colonne[i]=1;        // L'allele note '1' est dominant
    else colonne[i]=0;          // L'allele note '0' est dominant
   }
 return colonne;
}


//------------- Le generateur de chromo lui-meme-------------------------------

int main(int argc, char *argv[])
{
  int i,j;
  // Si erreur d'appel,on explique a l'utilisateur la liste des param. a donner
  if ((argc <  5) ||  (argc > 9)) 
    { 
      cout << "Usage: " << argv[0] << " choix n m f [p] [p'] [d] [l], où\n\n";
      cout << "       choix est le choix de repartition des marqueurs sur le chromosome: uniforme(u)ou aleatoire(a) \n";
      cout << "       n est le nombre de marqueurs\n";
      cout << "       m est le nombre d'individus dans l'échantillon\n";
      cout << "       f est le nom de base des fichiers de sortie \n";
      cout << "       p est la proba d'etre informatif, par defaut 1.0\n";
      cout << "       p' est la probabilité d'erreur, par defaut 0.0\n";
      cout << "       d est la probabilité de codominance, par defaut 0.5\n";
      cout << "       l la longueur totale du chromosome,  par defaut 1.0\n";
      return(1);
    }
  
  int NMarq = atoi(argv[2]);
  
  if (NMarq <= 1)
    { 
      cout << "Usage: " << argv[0] << " choix n m f [p] [p'] [d] [l], où\n\n";
      cout << "       choix est le choix de repartition des marqueurs sur le chromosome: uniforme(u)ou aleatoire(a) \n";
      cout << "       n est le nombre de marqueurs\n";
      cout << "       m est le nombre d'individus dans l'échantillon\n";
      cout << "       f est le nom de base des fichiers de sortie \n";
      cout << "       p est la proba d'etre informatif, par defaut 1.0\n";
      cout << "       p' est la probabilité d'erreur, par defaut 0.0\n";
      cout << "       d est la probabilité de codominance, par defaut 0.5\n";
      cout << "       l la longueur totale du chromosome,  par defaut 1.0\n";
      return(2);
    }
  
  int Taille = atoi(argv[3]);
  string OutName = argv[4];
  string RawName = new char [OutName.length()+5];
  string RamName = new char [OutName.length()+5];
  string TraceName = new char [OutName.length()+5];
  
  RawName = OutName+".raw";
  TraceName = OutName+".map";
  RamName = OutName+".ram";
  
  double Prob = 1.0;
  double ProbE = 0.0;
  double Codomin = 0.5;
  double Length = 1.0;
  
  if (argc >= 6)  Prob = atof(argv[5]);
  if (argc >= 7) ProbE = atof(argv[6]);
  if (argc >= 8)  Codomin = atof(argv[7]);
  if (argc >= 9)  Length = atof(argv[8]);

  marqueur* Marqueurs;
  if (strcmp(argv[1],"a")==0)
    {Marqueurs = genmarq(NMarq,Length);}
  else if (strcmp(argv[1],"u")==0)
    {Marqueurs = cartunif(NMarq,Length);}
  else 
    { 
      cout << "Usage: " << argv[0] << " choix n m f [p] [p'] [d] [l], où\n\n";
      cout << "       choix est le choix de repartition des marqueurs sur le chromosome:                 uniforme(u)ou aleatoire(a) \n";

      cout << "       n est le nombre de marqueurs\n";
      cout << "       m est le nombre d'individus dans l'échantillon\n";
      cout << "       f est le nom de base des fichiers de sortie \n";
      cout << "       p est la proba d'etre informatif, par defaut 1.0\n";
      cout << "       p' est la probabilité d'erreur, par defaut 0.0\n";
      cout << "       d est la probabilité de codominance, par defaut 0.5\n";
      cout << "       l la longueur totale du chromosome,  par defaut 1.0\n";
      return(3);
    }

  ofstream raw,trace,ram;
  ofstream & Raw = raw;
  ofstream & Ram = ram;
  ofstream & Trace = trace;

  for (i=0; i<256;i++)
       {
        for (j=0;j<256;j++) 
	  {tr[i][j] = 0;}
       }
  tr['0']['0'] = 0;
  tr['0']['1'] = 1;
  tr['1']['0'] = 1;
  tr['1']['1'] = 2;
  tr['0']['-'] = 3;
  tr['1']['-'] = 4;
  tr['-']['0'] = 3;
  tr['-']['1'] = 4;
  tr['-']['-'] = 5;
 
  Raw.open(RawName.c_str());
  Ram.open(RamName.c_str());
  Trace.open(TraceName.c_str());
  
  individu* Echantillon = gensample(Taille,NMarq,Marqueurs);
  individu* JDincomp = incomplet(Taille,NMarq,Echantillon,Prob,ProbE);
  int* col = dominance(NMarq,Codomin);
  //  On a cree la carte et l'echantillon complet
 
  pontemodele(NMarq,Marqueurs,Trace);
  Trace.close();

  // On trace la carte modele
  pontesample(Taille,NMarq,Echantillon,col,Ram);
  Ram.close();
  // On renvoie le jeu de donnees cree sans erreurs ni manques
  
  pontesample(Taille,NMarq,JDincomp,col,Raw);
  Raw.close();
  // On renvoie le jeu de donnees cree avec erreurs et manques

  delete[] Echantillon;
  delete[] JDincomp;
  delete[] Marqueurs;
  delete[] col;
  
  return(0);
}