#include"reconstruction.h"
#include<unistd.h>
#include<string.h>
     
/*
rabid:/scratch/xiang/ima_new/magicpool/bcarab_1st/fqs/*.tmp2
rabid:/scratch/xiang/ima_new/magicpool/bcarab_2nd/fqs/*.tmp2
rabid:/scratch/xiang/ima_new/magicpool/bcarab_3rd/recon/*.tmp2

reconstruction -a /data/www/19genomes/variants.tables -m ./transposons.gff -c ./MAGIC/

*/


int main( int argc, char **argv ) {

  ALLELE_DATA **add;
  int index;
  int c;
  int diploid=0;
  int forback=0;
  int penalty=10;
  int backcross=-1;
  int format=1, nchr=0;
  char *alleledir = NULL;
  char *outdir = NULL;
  char *writedir = NULL;
  char *calldir = NULL;
  char *balleledir = NULL;
  char *maskfile = NULL;
  char *zfile = NULL;
  char *Penalty = NULL;
  char *Species = strdup("mouse");
  char *opterr = 0;
  char *postdir = NULL;
     
  while ((c = getopt (argc, argv, "da:b:m:c:k:o:p:s:uw:z:fg:")) != -1)
    switch (c) {	
    case 'd':
      diploid=1;
      break;
    case 'b':
      sscanf( optarg, "%d", &backcross); // this shoudl be the index in the alleles data of the strain to use for backcrossing
      break;
    case 'a':
      alleledir = optarg;
      break;
    case 'c':
      calldir = optarg;
      break;
    case 'k':
      balleledir = optarg;
      break;
    case 'm':
      maskfile = optarg;
      break;
    case 'o':
      outdir = optarg;
      break;
    case 'p':
      Penalty = optarg;
      break;
    case 's':
      Species = optarg;
      break;
    case 'u':
      format = 0;
      break;
    case 'w':
      writedir = optarg;
      break;
    case 'z':
      zfile = optarg;
      break;
    case 'f':
      forback=1;
      break;
    case 'g':
      postdir = optarg;
      break;
    }
  if ( Penalty )
    sscanf( Penalty, "%d" , &penalty);
  if ( outdir == NULL )
    outdir = strdup("./");
  if (postdir == NULL)
    postdir = strdup("./");
  if ( zfile == NULL ) 
    zfile = strdup("mosaic.txt"); // outputfile of mosaic

  printf("calldir %s\n", calldir );
  printf("balleledir %s\n", balleledir );
  printf("penalty = %d\n", penalty); 	
  printf("species = %s\n", Species );
  char **chroms = TheChromosomes( Species, &nchr );
  
  if ( ( alleledir || balleledir ) && calldir && nchr > 0  ) {
    int s, chr, nseq, ii, jj=0;
    ALLELE_DATA ** add;
    if ( alleledir ) 
      add = ReadGenomeAlleleData( alleledir, maskfile, format, nchr, chroms ) ;
    else if ( balleledir )
      add = ReadAllBinaryAlleleData(balleledir, 0, nchr, chroms );


    printf("ReadSeqCalls");

    SEQ_CALL **call = ReadSeqCalls( calldir , "", add, nchr, chroms, &nseq );
    FILE *fp = fopen(zfile,"w");

    if ( writedir != NULL ) 
      WriteSeqCalls( call, nseq, writedir );

    if (!forback)
      fprintf(fp, "magic	chr	acc	from.bp	to.bp	from.site	to.site	len.bp	sites	errors	error.site	error.bp\n");

    printf("Main function call \n");
    
    if ( forback )
      for( s=0;s<nseq;s++) {
        printf("Genome %s out of %d \n", call[s]->name, nseq);
        char filename[500]; 
        char *ending = ".posteriors";
        char version[10];
        int vv = sprintf(version, "%d", diploid);
        int ret = sprintf(filename, "%s%s%s%s", postdir, call[s]->name, ending, version);
        FILE *fp = fopen(filename, "w");
        for(chr=0;chr<nchr;chr++) {
	  if ( diploid ) 
	    ForwardBackwardDiploid( call[s], penalty, chr, chroms[chr], fp );
	  else
	    ForwardBackward( call[s], penalty, chr, chroms[chr], fp);
        }
      }
    else
      for( s=0;s<nseq;s++) {
        for(chr=0;chr<nchr;chr++) {
	  if ( diploid ) 
	    ReconstructMosaicDiploid( call[s], penalty, chr, chroms[chr], fp );
	  else 
	    ReconstructMosaic( call[s], penalty, chr, chroms[chr], fp , backcross);
        } 
      }
    fclose(fp);
    ImputeGenomes( add, nchr, chroms, call, nseq, outdir );
  }
}