// Un generateur d'echantillon d'hybrides irradies diploides...



//----------------- 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>


//------------------- Les classes utilisees -----------------------------------
typedef char* gamete;

class marqueur
{
public:
    double position;
    double cassure;
  
  marqueur()  { position = cassure = 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,'A',n);pere[n]=0;                               mere=new char[n+1];memset(mere,'A',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 proba de cassure a partir de la distance en ray8000 -------

double Ray2Theta(double ray)
{
  return  (1.0-exp(-ray));
}


// ---------------- 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 proba de cassure
    double ray;
    for (i=1; i<nombre; i++)
      {
	ray = Eux[i].position-Eux[i-1].position;
	Eux[i].cassure = Ray2Theta(ray);
      }
    
    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 proba de cassure
    double ray;
    for (i=1; i<nombre; i++)
      {
	ray = Eux[i].position-Eux[i-1].position;
	Eux[i].cassure = Ray2Theta(ray);
      }    
    
    return Eux;
    // On obtient en fait une carte genetique
}


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

gamete gengamete(int nombre, marqueur* modele, double retention)
{
    gamete fragment;
    char prec;
    int i; 
       
         fragment=new char[nombre+1];
         memset(fragment,'A',nombre);fragment[nombre]=0;         
         // On cree un chromosome aleatoire d'apres le modele de carte genere 
          if (randomax() < retention)
            fragment[0] ='H';
          prec = fragment[0];  
          
          for(i=1; i<nombre; i++)
  	     {
	      if (randomax() < modele[i].cassure)
	         {
                  if (randomax() < retention)
                      prec='H';
                  else
                      prec='A';  
	         }
              fragment[i]=prec;                                     
	     }
      return fragment;
}
 

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

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

    for (i=0; i< taille; i++)
      {
       sample[i].pere=gengamete(nombre, modele, retention );
       sample[i].mere=gengamete(nombre, modele, retention );
       // 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]=((echantillon[i].pere[j]=='A')? 'H':'A'); 
                 else 
                   sample[i].pere[j]=echantillon[i].pere[j];
		}
          else  sample[i].pere[j]='-';
          if (randomax()<proba)
	        {
	         if (randomax() < probae)
		   sample[i].mere[j]=((echantillon[i].mere[j]=='A')? 'H':'A');
                 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,double retention,marqueur* marqs,ofstream & out)
{
  int i;
  
  out << "CARTE:            proba de retention : "<<retention << "\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 << "\tcassure : " << marqs[i].cassure;
	}
      out << endl;
    }
}


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

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

  out << "data type radiated hybrid\n";
  out << taille << " " << nombre << " 0 0\n";
    
   for (i=0; i<nombre; i++)
     {
       out << "*M" << i+1 << "\t";
       
       
       for (j=0;j<taille; j++)
	 {
         if ((echantillon[j].pere[i]=='H')||(echantillon[j].mere[i]=='H'))
             out << "H";
         else if ((echantillon[j].pere[i]=='-')&&(echantillon[j].mere[i]=='-'))
             out << "-";
         else 
             out << "A";
         } 
       out << endl;
     }
}


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

int main(int argc, char *argv[])
{
  // Si erreur d'appel,on explique a l'utilisateur la liste des param. a donner
  if ((argc <  6) ||  (argc > 9)) 
    { 
      cout << "Usage: " << argv[0] << " choix n m r f [p] [p'] [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 << "       r est la probabilite de retention\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 << "       l la longueur totale du chromosome,  par defaut 20ray\n";
      return(1);
    }
  
  int NMarq = atoi(argv[2]);
  
  if (NMarq <= 1)
    { 
      cout << "Usage: " << argv[0] << " choix n m r f [p] [p'] [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 << "       r est la probabilite de retention\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 << "       l la longueur totale du chromosome,  par defaut 20ray\n";
      return(2);
    }
  
  int Taille = atoi(argv[3]);
  double Retention = strtod(argv[4],(char **)NULL);
  string OutName = argv[5];
  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 Length = 20.0;
  
  if (argc >= 7)  Prob = atof(argv[6]);
  if (argc >= 8) ProbE = 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 r f [p] [p'] [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 << "       r est la probabilite de retention\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 << "       l la longueur totale du chromosome,  par defaut 20ray\n";
      return(3);
    }

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

   Raw.open(RawName.c_str());
  Ram.open(RamName.c_str());
  Trace.open(TraceName.c_str());
  
  individu* Echantillon = gensample(Taille,NMarq,Marqueurs,Retention);
  individu* JDincomp = incomplet(Taille,NMarq,Echantillon,Prob,ProbE);
  //  On a cree la carte et l'echantillon complet
 
  pontemodele(NMarq,Retention,Marqueurs,Trace);
  Trace.close();

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

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