#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 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;
     
  while ((c = getopt (argc, argv, "da:b:m:c:k:o:p:s:uw:z:")) != -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;
    }
  if ( Penalty )
    sscanf( Penalty, "%d" , &penalty);
  if ( outdir == NULL )
    outdir = 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 );

//    for (ii = 0;ii<20;ii++){
//        for(jj = 0; jj<24; jj++){
//            printf('%c', add[ii][jj]);
//        }
//    printf('\n');
//    }

    printf("ReadSeqCalls");

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

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

    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");
    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 );
  }
}