#include"reconstruction.h"
#include<unistd.h>
#include<string.h>
#include<stdlib.h>

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

  int c, max_annotations=10000;
  char phenotype_file[256];
  char phenotype_name[2256];
  char *phenotype_names[256], *pn;
  char dir[256];
  char alleledir[256];
  char outfile[256];
  char qtlfile[256];
  char larsfile[256];
  char annotationfile[256];
  char outdir[256];
  char estfile[256];
  char permutationfile[256];
  char species[256];
  BINARY_DATA **bd;
  double threshold=0, thresh=4;
  SCAN_DATA *qtl=NULL;
  PHENOTYPE_TABLE **ptt;
  ALLELE_DATA **ad=NULL;
  int uncompressed = 1;
  int haplotype = 0;
  int nphen=0, n, nperm=0, nchr=0;

  threshold=thresh;
  outfile[0] = 0;
  qtlfile[0] = 0;
  annotationfile[0] = 0;
  larsfile[0] = 0;
  permutationfile[0] = 0;
  strcpy(alleledir, "./");
  strcpy(dir, "./");
  strcpy(outdir, "./");
  strcpy(species, "arabidopsis");

  while ((c = getopt (argc, argv, "d:f:p:t:a:w:n:s:chH")) != -1)
    switch (c) {
    case 'a':
      strcpy(alleledir, optarg);
      break;
    case 'd':
      strcpy(dir, optarg);
      break;
    case 'f':
      strcpy(phenotype_file, optarg);
      break;
    case 'p':
      strcpy(phenotype_name, optarg);
      break;
    case 'w':
      strcpy( outdir, optarg );
      break;
    case 's':
      strcpy(species, optarg);
      break;
    case 't':
      sscanf( optarg, "%lf", &threshold);
      break;
    case 'h':
      haplotype = 1;
      break;
    case 'c':
      uncompressed = 0;
      break;
    case 'n':
      sscanf( optarg, "%d", &nperm );
      break;
    case 'H':
      printf( "genome_scan \n"
"-f  <filename>           tab-delimited file of phenotypes\n"
"-p  <phenotype>          phenotype column(s) in phenotype file to analyse (comma separated list)\n"
"-d  <./>                 directory with binary genotype files\n"
"-a  <./>                 directory with binary alleles file\n"
"-w  <./>                 output directory\n"
"-t  <5>                  logP threshold used in forward selection and annotation\n"
	      "-h                       use haplotype tests and files\n"
"-c                       use compressed sdp files\n"
     );
      break;
    }
  
  printf("uncompressed %d\n", uncompressed);
  printf("haplotype %d\n", haplotype );
  printf("species = %s\n", species );
  char **chroms = TheChromosomes( species, &nchr );
  
  pn = strtok(phenotype_name, ",");
  if ( pn == NULL) {
    nphen = 1;
    phenotype_names[0] = phenotype_name;
  }
  else {
    nphen = 1;
    phenotype_names[0] = pn;
    while( pn = strtok(NULL, ",") ) {
      printf( "nphen %d phen %s\n", nphen, pn );
      phenotype_names[nphen++] = pn;
    }
  }
  
  ptt = ReadPhenotypes(phenotype_file, phenotype_names, nphen);
  bd = ReadBinaryGenome( nchr, dir, uncompressed, haplotype );
  if ( (!haplotype) && alleledir ) 
    ad = ReadAllBinaryAlleleData( alleledir, 1, nchr, chroms);
  
  for(n=0;n<nphen;n++) {
    PERMUTATIONS *perm=NULL;
    SCAN_STATISTICS **stats;

    if ( nperm > 0 ) 
      perm = MakePermutations( ptt[n], bd[0], nperm );

    stats = GenomeScan( ptt[n], bd, nchr, perm );
    outfile[0] = qtlfile[0] = estfile[0] = annotationfile[0] = 0;
    if ( ! haplotype ) {
      sprintf(outfile, "%s/%s.logP.txt", outdir, phenotype_names[n]);
      sprintf(estfile, "%s/%s.estimate.txt", outdir, phenotype_names[n]);
      sprintf(qtlfile, "%s/%s.mult.txt", outdir, phenotype_names[n]);
      sprintf(annotationfile, "%s/%s.annotated.txt", outdir, phenotype_names[n]);
      sprintf(larsfile, "%s/%s.lars.txt", outdir, phenotype_names[n]);
    }
    else {
      threshold = 0.0;
      sprintf(outfile, "%s/%s.haplotype.logP.txt", outdir, phenotype_names[n]);
      sprintf(estfile, "%s/%s.haplotype.estimate.txt", outdir, phenotype_names[n]);
      sprintf(qtlfile, "%s/%s.haplotype.mult.txt", outdir, phenotype_names[n]);
    }
    
    WriteGenomeScan( stats, nchr, bd, outfile, threshold);

    WriteLarsData( ptt[n], stats, nchr, bd, larsfile, threshold );
    if ( (!haplotype) && ad) {
      WriteGenomeScanEstimate( stats, nchr, bd, ad, estfile, threshold);
      WriteAnnotatedScan( ptt[n], nchr, stats, bd, ad, thresh, max_annotations, annotationfile );
      qtl = ForwardSelection( ptt[n], nchr, bd, ad, stats, thresh, qtlfile );
    }
    else if ( haplotype ) {
      WriteHaplotypeGenomeScanEstimate( stats, nchr, bd, estfile, threshold);
    }
    if ( perm != NULL ) {
      int i;
      FILE *pfp = NULL;
      //      p = permutation_pvalue( logP, ptt[n], bd, nchr, nperm, &q50, &q90, &q95, &mx );
      printf( "%s mx %e pval %e q50 %e q90 %e q95 %e\n", ptt[n]->name, perm->mx, perm->pval, perm->q50, perm->q90, perm->q95 ); 
      if ( haplotype )
	sprintf(permutationfile, "%s/%s.perms_hap.txt", outdir, phenotype_names[n]);
      else 
	sprintf(permutationfile, "%s/%s.perms.txt", outdir, phenotype_names[n]);
      pfp = fopen(permutationfile,"w");
      fprintf(pfp, "%s mx %e pval %e q50 %e q90 %e q95 %e\n", ptt[n]->name, perm->mx, perm->pval, perm->q50, perm->q90, perm->q95 ); 
      for( i=0;i<perm->nperm;i++) {
	fprintf(pfp, "%d %e\n", i, perm->maxlogP[i]);
      }
      fclose(pfp);
      FreePermutations(perm);
    }
  }
}