#include "pseud.h"
#include <math.h>

static void
ExchangeAlignInf(AlignIDPtr temp_align, int i, int j);
static short
CompareAlign(AlignIDPtr	temp_align, AlignIDPtr align_id);
static void
CopyAlignID(AlignIDPtr dist_align, AlignIDPtr src_align);
static void
SumAlignID(AlignIDPtr align_id, AlignIDPtr temp_align);

AlignIDPtr
AlignIDFree(AlignIDPtr align_id)
{
	if(align_id==NULL){return NULL;}
	if(align_id->gap_length){
   	align_id->gap_length = MemFree(align_id->gap_length);
   }
	if(align_id->pos_vector){
   	align_id->pos_vector = MemFree(align_id->pos_vector);
   }
   align_id = MemFree(align_id);

   return align_id;
}

AlignIDPtr
AlignIDNew(short nmax)
{
	AlignIDPtr	align_id;

   align_id = (AlignIDPtr)calloc(1,sizeof(AlignID));
   if(align_id==NULL){return NULL;}

   align_id->gap_length = (short *)calloc(nmax, sizeof(short));
   if(align_id->gap_length==NULL){
   	align_id = MemFree(align_id);  return NULL;
   }
   align_id->pos_vector = (short *)calloc(nmax, sizeof(short));
   if(align_id->pos_vector==NULL){
   	align_id = AlignIDFree(align_id);  return NULL;
   }
   align_id->num_visit = 0;

	return align_id;
}

StockAlignPtr
StockAlignDestruct(StockAlignPtr stock_align)
{
	long int   i;
   AlignIDPtr PNTR align_id_array;

	if(stock_align==NULL){return NULL;}

   if(stock_align->align_id_array){
   	align_id_array = stock_align->align_id_array;
      for(i=0;i<(stock_align->max_idnum);i++){
      	align_id_array[i] = AlignIDFree(align_id_array[i]);
      }
      align_id_array = MemFree(align_id_array);
   }

   stock_align = MemFree(stock_align);

   return stock_align;
}


StockAlignPtr
StockAlignAllocate(long int max_idnum, short nmax, short select_num)
{
   long int      i;
	StockAlignPtr stock_align;
   AlignIDPtr PNTR align_id_array;

	stock_align = (StockAlignPtr)calloc(1,sizeof(StockAlign));
   if(stock_align==NULL){return NULL;}

   stock_align->max_idnum = max_idnum;
   stock_align->current_id = 0;
   stock_align->total_accepted = 0;
   stock_align->select_num = select_num;
   stock_align->select_id0 = select_num;

   align_id_array = (AlignIDPtr PNTR)calloc(max_idnum,sizeof(AlignIDPtr));
   if(align_id_array==NULL){stock_align = MemFree(stock_align);  return NULL;}
   for(i=0;i<max_idnum;i++){
   	align_id_array[i] = AlignIDNew(nmax);
   	if(align_id_array[i]==NULL){stock_align = MemFree(stock_align);  return NULL;}
   }
   stock_align->align_id_array = align_id_array;

   return stock_align;

}

/*group summarized id with same align*/
static void
ExchangeAlignInf(AlignIDPtr temp_align, int i, int j)
{
	short temp;
   short *gap_length, *pos_vector;

   gap_length = temp_align->gap_length;
   pos_vector = temp_align->pos_vector;
   /*exchange pos*/
   temp = pos_vector[j];
   pos_vector[j] = pos_vector[i];
   pos_vector[i] = temp;
   /*exchange gap_length*/
   temp = gap_length[j];
   gap_length[j] = gap_length[i];
   gap_length[i] = temp;

   return;
}

void
ArrayToAlign(CandStructPtr temp_cand, AlignIDPtr temp_align, int num_t)
{
   int    t, i, j, p, nsum;
   short  *id_array, *gap_length;
   short  *pos_vector, *pos_align;

   id_array = temp_cand->id_array;
   pos_vector = temp_cand->pos_vector;
   gap_length = temp_align->gap_length;
   pos_align = temp_align->pos_vector;

   nsum = temp_cand->ni + temp_cand->nd;
   temp_align->num_gap = nsum;
   /*printf("id_array:");
   for(t=0;t<num_t;t++){printf("%d,",id_array[t]);}
   printf("\n");
   printf("pos_vector:");
   for(i=0;i<nsum;i++){printf("%d,",pos_vector[i]);}
   printf("\n"); */

   p=0;
   for(t=0;t<num_t;t++){
   	if(id_array[t]>0){
      	for(i=0;i<p;i++){
         	if(pos_vector[i]>pos_vector[p]){pos_vector[i]+=id_array[t];}
         }
         p++;
		}
   	if(id_array[t]<0){
      	for(i=(p+1);i<nsum;i++){
         	if(pos_vector[i]>pos_vector[p]){pos_vector[i]-=id_array[t];}
         }
         p++;
      }
   }
   p=0;
   for(t=0;t<num_t;t++){
   	if(id_array[t]!=0){
         gap_length[p]=id_array[t];
         pos_align[p]=pos_vector[p];
      	p++;
      }
   }
   if(p!=nsum){printf("Warning!p(%d)!=nsum(%d)\n",p,nsum);getch();}
   /*display temp_align*/
   /*printf("gap_length:");
   for(i=0;i<nsum;i++){printf("%d,",gap_length[i]);}
   printf("\n");
   printf("pos_align:");
   for(i=0;i<nsum;i++){printf("%d,",pos_align[i]);}
   printf("\n");
   getch(); */

   /*sort by position*/
   if((temp_align->num_gap)<=1){return;}
   for(i=1;i<(temp_align->num_gap);i++){
   	for(j=0;j<i;j++){
      	if(pos_align[i]<pos_align[j]){
         	ExchangeAlignInf(temp_align, i, j);
         }
      }
   }
   /*printf("sort by position\n");
   printf("gap_length:");
   for(i=0;i<nsum;i++){printf("%d,",gap_length[i]);}
   printf("\n");
   printf("pos_vector:");
   for(i=0;i<nsum;i++){printf("%d,",pos_align[i]);}
   printf("\n");
   getch();*/

   return;
}

void
IDdataToAlign(ChainIDPtr chain_id, AlignIDPtr temp_align)
{
   temp_align->lnL = chain_id->lnL;
   temp_align->num_visit = chain_id->num_visit;
   temp_align->lnProbAS = chain_id->lnProbAS;

	return;
}

static short
CompareAlign(AlignIDPtr	temp_align, AlignIDPtr align_id)
{
   int   i, num_gap;
	short  *gap_length1, *gap_length2;
   short  *pos_vector1, *pos_vector2;

   if((temp_align->num_gap)!=(align_id->num_gap)){return -1;}
   num_gap = temp_align->num_gap;
   gap_length1 = temp_align->gap_length;
   gap_length2 = align_id->gap_length;
   for(i=0;i<num_gap;i++){
   	if(gap_length1[i]!=gap_length2[i]){return -1;}
   }
   pos_vector1 = temp_align->pos_vector;
   pos_vector2 = align_id->pos_vector;
   for(i=0;i<num_gap;i++){
   	if(pos_vector1[i]!=pos_vector2[i]){return -1;}
   }

   return 1;
}

static void
CopyAlignID(AlignIDPtr dist_align, AlignIDPtr src_align)
{
   int   i;
   short *dist_gaplen, *src_gaplen;
   short *dist_pos, *src_pos;

   dist_align->num_gap = src_align->num_gap;
   dist_gaplen = dist_align->gap_length;
   src_gaplen = src_align->gap_length;
   dist_pos = dist_align->pos_vector;
   src_pos = src_align->pos_vector;
   for(i=0;i<(src_align->num_gap);i++){dist_gaplen[i]=src_gaplen[i];}
   for(i=0;i<(src_align->num_gap);i++){dist_pos[i]=src_pos[i];}
   dist_align->lnL = src_align->lnL;
   dist_align->num_visit = src_align->num_visit;

   return;
}

static void
SumAlignID(AlignIDPtr align_id, AlignIDPtr temp_align)
{
	double temp;

   temp = 0;
   /*temp += pow(10.0, align_id->lnL)*(align_id->num_visit);
   temp += pow(10.0, temp_align->lnL)*(temp_align->num_visit);
   temp /= (align_id->num_visit + temp_align->num_visit); */
   /*Modified on Feb.16,2001*/
   temp += pow(10.0, align_id->lnL);
   temp += pow(10.0, temp_align->lnL);
   align_id->lnL = log10(temp);

   align_id->num_visit += temp_align->num_visit;

   return;
}

short
AlignSummarizedID(StockAlignPtr stock_align, StockIDPtr summarized_id, short lmax, short nmax, int num_t)
{
   long int i,j,p,max_idnum2,found;
   AlignIDPtr PNTR align_id_array;
   AlignIDPtr      temp_align;
   ChainIDPtr PNTR chain_id_array;
   ChainIDPtr      chain_id;
	CandStructPtr temp_cand;

   max_idnum2 = stock_align->max_idnum;

	temp_cand = CandStructAllocate(nmax, num_t);
   if(temp_cand==NULL){
   	printf("no memory for temp_cand\n");getch();return 1;
   }
   temp_align = AlignIDNew(nmax);
   if(temp_align==NULL){
   	printf("no memory for temp_align\n");getch();
   	return 1;
   }

   align_id_array = stock_align->align_id_array;
   chain_id_array = summarized_id->chain_id_array;
   p=0;
   for(i=0;i<(summarized_id->current_id);i++){
      /*align with chain_id*/
   	chain_id = chain_id_array[i];
		IDToArray(chain_id, temp_cand, lmax, nmax, num_t);
      ArrayToAlign(temp_cand, temp_align, num_t);
      IDdataToAlign(chain_id, temp_align);
      /*Compare Align*/
      if(p==0){found=-1;}
      else{
         found = -1;
      	for(j=0;j<p;j++){
         	found = CompareAlign(temp_align, align_id_array[j]);
         	if(found>=0){found = j; break;}
         }
      }
      if(found<0){
      	CopyAlignID(align_id_array[p], temp_align);
         p++;
         if(p>=max_idnum2){
         	printf("Warning!currrent_id(%d)>max_idnum(%d)\n",p,max_idnum2);getch();
   			temp_cand = CandDestruct(temp_cand);
   			temp_align = AlignIDFree(temp_align);
            stock_align->current_id = p; return 0;
         }
      }
      else{
      	SumAlignID(align_id_array[found], temp_align);
      }
   }
   stock_align->current_id = p;

   temp_cand = CandDestruct(temp_cand);
   temp_align = AlignIDFree(temp_align);

   return 0;
}

void
DiplayAlignID(AlignIDPtr align_id, FILE *fout)
{
   int    i;
	short  *gap_length;
   short  *pos_vector;

   fprintf(fout,"num_gap=%d\n",align_id->num_gap);
   gap_length = align_id->gap_length;
   pos_vector = align_id->pos_vector;
   fprintf(fout,"gap_length:");
   for(i=0;i<(align_id->num_gap);i++){
   	fprintf(fout,"%d,",gap_length[i]);
   }
   fprintf(fout,"\t");
   fprintf(fout,"pos_vector:");
   for(i=0;i<(align_id->num_gap);i++){
   	fprintf(fout,"%d,",pos_vector[i]);
   }
   fprintf(fout,"\n");
   fprintf(fout,"lnL=%10.4f, num_visit=%d\n",align_id->lnL, align_id->num_visit);

   return;
}
void
DisplayStockAlign(StockAlignPtr stock_align, FILE *fout)
{
	int   i;
   AlignIDPtr PNTR align_id_array;
   AlignIDPtr      align_id;
   short         select_num;

   align_id_array = stock_align->align_id_array;
   if((stock_align->select_num)>(stock_align->current_id)){select_num=(short)(stock_align->current_id);}
   else{select_num = stock_align->select_num;}
   for(i=0;i<select_num;i++){
   	fprintf(fout,"No.%ld:",(i+1));
      DiplayAlignID(align_id_array[i], fout);
   }

   /*probability of A->S1*/
   for(i=0;i<select_num;i++){
      align_id = align_id_array[i];
		fprintf(fout,"No. %d: ln(prob) = %10.4f\n",(i+1), align_id->lnProbAS);
   }
   fprintf(fout,"total_accepted = %ld\n",stock_align->total_accepted);
   fprintf(fout,"current_id = %ld\n",stock_align->current_id);
   fprintf(fout,"select_num = %d\n",stock_align->select_num);
   fprintf(fout,"ave_ln(prob) = %10.4f\n",stock_align->Ave_lnAS);

  	return;
}

void
SortStockAlign(StockAlignPtr stock_align, short nmax)
{
   AlignIDPtr PNTR align_id_array;
   AlignIDPtr      align_id1;
   AlignIDPtr      align_id2;
   AlignIDPtr      align_temp;
   long int   i, j, k, current_id;
   short      select_num;
   double     prob, ave_prob, account, max_account;

   align_id_array = stock_align->align_id_array;
   current_id = stock_align->current_id;

   if(current_id==1){/*only one candidate, no need of sorting*/
   	align_id1 = align_id_array[0];
   	stock_align->total_accepted = align_id1->num_visit;
		align_id1->lnProbAS = align_id1->lnL;
      stock_align->Ave_lnAS = align_id1->lnL;
		stock_align->select_num = 1;
      return;
   }

   align_temp = AlignIDNew(nmax);

   for(i=1;i<current_id;i++){
   	align_id1 = align_id_array[i];
   	for(j=0;j<i;j++){
      	align_id2 = align_id_array[j];
         if((align_id1->num_visit)>(align_id2->num_visit)){
   			/*Copy align_id1 to align_temp*/
            CopyAlignID(align_temp, align_id1);
            for(k=i-1;k>=j;k--){
					CopyAlignID(align_id_array[k+1], align_id_array[k]);
            }
   			/*copy chain_temp to align_id2*/
            CopyAlignID(align_id2, align_temp);
         }
      }
   }

   /*count total number of accepted chain-t-pos*/
   stock_align->total_accepted = 0;
   for(i=0;i<current_id;i++){
   	align_id1 = align_id_array[i];
   	stock_align->total_accepted += align_id1->num_visit;
   }

   /*Calculate probability of A->S1*/
   if((stock_align->select_num)>current_id){
   	printf("select_num(%d)>current_id(%d)\n",(stock_align->select_num),current_id);
      select_num = (short)current_id;
   }
   else{select_num=stock_align->select_num;}

	/*modified on Apr.12, 2001*/
   /*The minimum number of the top ranked candidates which collectively
   account for 90 percent or more used to estimate A to S probability*/
   account = 0.0;
   max_account = (double)(stock_align->total_accepted)*0.90;
   ave_prob = 0;
   for(i=0;i<select_num;i++){
      align_id1 = align_id_array[i];
      prob = 1;
      prob /= align_id1->num_visit;
      prob *= stock_align->total_accepted;
      prob *= pow(10.0, align_id1->lnL);
      ave_prob += prob;
      if(prob<=0){printf("prob(%lf)<=0\n",prob);prob = MinlnL;}
      else{prob = log10(prob);}
      align_id1->lnProbAS = prob;
      account += (align_id1->num_visit);
      if(account>max_account){printf("the minimum number is %d\n",i); /*getch();*/select_num=(short)(i+1);break;}
   }
   ave_prob /= select_num;
   stock_align->select_num = select_num;
   if(ave_prob<=0){printf("ave_prob(%lf)<=0\n",ave_prob);stock_align->Ave_lnAS = MinlnL;}
   else{stock_align->Ave_lnAS = log10(ave_prob);}

   align_temp = AlignIDFree(align_temp);

	return;
}

void
InitializeStockAlign(StockAlignPtr stock_align, short nsum)
{
   AlignIDPtr PNTR align_id_array;
   AlignIDPtr      align_id;
	short  *gap_length;
   short  *pos_vector;
   int        i, j;

   align_id_array = stock_align->align_id_array;
   for(i=0;i<stock_align->current_id;i++){
      align_id = align_id_array[i];
   	align_id->num_gap=0;
		align_id->lnL=0;
   	align_id->num_visit=0;
   	align_id->lnProbAS=0;
      gap_length = align_id->gap_length;
      for(j=0;j<nsum;j++){gap_length[j]=0;}
      pos_vector = align_id->pos_vector;
      for(j=0;j<nsum;j++){pos_vector[j]=-1;}
   }
   stock_align->current_id=0;
   stock_align->select_num = stock_align->select_id0;

   return;
}