/**************************************************************************
 * Project: phylogenic tree model with insetions and/or deletions
 *
 * FILE:pseud_main.c
 *
 * File description: main program for the project
 *
 *
 *
 * Programmer:Manami and Kazuhisa Nishizawa
 *
 * Date: March 22, 2001
 *         modified on Feb.20,2003
 *
 *************************************************************************
*/

#include <stdio.h>
#include <stdlib.h>
#include <ctype.h>
#include <conio.h>
#include <time.h>
#include <alloc.h>
#include <string.h>
#include <math.h>
#include "pseud.h" 

/*----------------------------------------------*/
int main()
{
   char   *outfile, *infile, *subsfile, *seqfile;
   FILE   *fseq, *fout;
   char   title[15];
   int    num, num_pairs, count;
   OptionsBlkPtr    options;
	StockIDPtr       stock_id;
	StockIDPtr       summarized_id;
	SubsMatrixPtr    subs_matrix;
   short            failed;
	StockAlignPtr    stock_align;
   double           sum_AS;

   outfile = (char *)calloc(NameSize,sizeof(char));
   infile = (char *)calloc(NameSize,sizeof(char));
   subsfile = (char *)calloc(NameSize,sizeof(char));
   seqfile = (char *)calloc(NameSize,sizeof(char));


   strcpy(subsfile,"subs.txt");
   printf("Enter sequence file name(sample1.txt):");
   scanf("%s",seqfile);
   printf("Enter output file name(out.txt):");
   scanf("%s",outfile);
   strcpy(infile,"common.txt"); 

   if((fseq = fopen(seqfile,"rt"))==NULL){
   	printf("Cannot open sequence file(press any key)\n");getch(); return 1;
   }
   free(seqfile);
   if((fout = fopen(outfile,"wt"))==NULL){
   	printf("Cannot open output file(press any key)\n");getch(); return 1;
   }
   free(outfile);

   /*Read sequence and parameters*/
   fprintf(fout,"Common optional parameters\n");
   options = GetOptions(infile, fout);
   free(infile);
   if(options==NULL){
   	printf("Get Options failed(press any key)\n");getch();return 1;
   }


   SetSeed ((unsigned int)time(NULL));

	subs_matrix = CalculateSubsMatrix(options, subsfile);
   if(subs_matrix==NULL){
   	printf("CalculateSubsMatrix failed(press any key)\n");getch();return 1;
   }
   printf("CalculateSubsMatrix done\n");

   stock_id = StockIDAllocate(options->max_idnum, options->nmax, options->select_num);
   if(stock_id==NULL){
   	printf("Stock ID Allocate failed(press any key)\n");getch();
      return 1;
   }
   printf("Stock ID successfully allocated\n");
	summarized_id = StockIDAllocate((options->max_idnum)/2, options->nmax, options->select_num);
   if(summarized_id==NULL){
   	printf("No memory for summarized_id\n");getch();return 1;
   }
   printf("Summarized ID successfully allocated\n");

   stock_align = StockAlignAllocate(summarized_id->max_idnum, options->nmax, options->select_num);
   if(stock_align==NULL){
   	printf("No memory for stock_align\n");getch();return 1;
   }

   fscanf(fseq, "%[^=]=%d, ",title, &num_pairs);
   count = 0;
   sum_AS = 0.0;
   for(num=0;num<num_pairs;num++){
   	fprintf(fout,"###seq_no.=%d\n",(num+1));
		ReadOneSequencePair(fseq, options, fout);
		ReCalculateOptionsForNewPair(options, fout);

		RecalculateSubsMatrix(subs_matrix, options/*, subsfile*/);
      printf("re-calculate subs matrix done\n");
   	InitializeStockID(stock_id, (short)(2*options->nmax));
   	InitializeStockID(summarized_id, (short)(2*options->nmax));
		InitializeStockAlign(stock_align, (short)(2*options->nmax));
   	/*failed=ProcessMCMC(stock_id, options, subs_matrix, fout);*/
		/*failed = ModProcessMCMC(stock_id, options, subs_matrix, options->target_flag, fout);*/

      failed=ExModProcessMCMC(stock_id, options, subs_matrix, fout);
   	if(failed){printf("MCMC failed(press any key)\n");getch();}
      else{
   		SortStockID(stock_id, options->nmax);
			/*DisplayStockID(stock_id, options, fout); */

   		SummarizedStockID(summarized_id, stock_id, options->nmax, options->num_t);
   		SortStockID(summarized_id, options->nmax);
         fprintf(fout,"current_id=%ld, \n",stock_id->current_id);
         fprintf(fout,"current_id=%ld, \n",summarized_id->current_id);
			DisplayStockID(summarized_id, options, fout);

			failed = AlignSummarizedID(stock_align, summarized_id, options->lmax, options->nmax, options->num_t);
   		if(failed){printf("stock_align failed\n");getch();return 1;}
			SortStockAlign(stock_align, options->nmax);
			DisplayStockAlign(stock_align, fout); 
      	count++;
      	sum_AS += stock_align->Ave_lnAS;
   		fprintf(fout,"No.%d:Ave_lnAS = %10.4f\n",count,stock_align->Ave_lnAS);
      }
   }
   printf("loop done\n");
   printf("The number of successfully calculated pairs is %d\n",count);
   printf("sum_AS = %10.4f\n",sum_AS);
   fprintf(fout,"The number of successfully calculated pairs is %d\n",count);
   fprintf(fout,"sum_AS = %10.4f\n",sum_AS);

   fclose(fout);fclose(fseq);

	options = OptionsBlkFree(options);
   /*stock_id = StockIDFree(stock_id);*/
   /*summarized_id = StockIDFree(summarized_id);*/
   /*stock_align = StockAlignDestruct(stock_align); */
   subs_matrix = SubsMatrixDestruct(subs_matrix); 

   /*free(subsfile); */

   return 0;

}