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

static short
CompareChainID(CandStructPtr new_cand, ChainIDPtr chain_id, short lmax, int num_t);
static double
CalculateTPOQijQjiRatio(CandStructPtr new_cand, CandStructPtr current_cand, PoolStructPtr pool);
static double
CalculateConserveTPOQijQjiRatio(CandStructPtr new_cand, CandStructPtr current_cand, PoolStructPtr pool);

ChainIDPtr
ChainIDFree(ChainIDPtr chain_id)
{
	if(chain_id==NULL){return NULL;}

   if(chain_id->pos_vector){
   	chain_id->pos_vector = MemFree(chain_id->pos_vector);
   }
   chain_id = MemFree(chain_id);

	return chain_id;
}

ChainIDPtr
ChainIDNew(short nsum)
{
	ChainIDPtr chain_id;

   chain_id = (ChainIDPtr)calloc(1,sizeof(ChainID));
   if(chain_id==NULL){return NULL;}
   chain_id->pos_vector = (short *)calloc(nsum,sizeof(short));
   if(chain_id->pos_vector==NULL){
   	chain_id = ChainIDFree(chain_id); return NULL;
   }
   chain_id->num_visit=0;

   return chain_id;
}


StockIDPtr
StockIDFree(StockIDPtr stock_id)
{
   int  i;
   ChainIDPtr PNTR chain_id_array;

	if(stock_id==NULL)return NULL;

   chain_id_array = stock_id->chain_id_array;
   for(i=0;i<stock_id->max_idnum;i++){
   	chain_id_array[i] = ChainIDFree(chain_id_array[i]);
   }
   chain_id_array = MemFree(chain_id_array);

   stock_id = MemFree(stock_id);

   return stock_id;
}


StockIDPtr
StockIDAllocate(long int max_idnum, short nmax, short select_num)
{
   long int        i;
	StockIDPtr stock_id;
   ChainIDPtr PNTR chain_id_array;

   stock_id = (StockIDPtr)malloc(sizeof(StockID));
   if(stock_id==NULL){return NULL;}

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

   chain_id_array = (ChainIDPtr PNTR)malloc(max_idnum*sizeof(ChainIDPtr));
   if(chain_id_array==NULL){
   	stock_id = StockIDFree(stock_id);return NULL;
   }
   for(i=0;i<max_idnum;i++){
   	chain_id_array[i]= ChainIDNew((short)(2*nmax));
      if(chain_id_array[i]==NULL){stock_id = StockIDFree(stock_id);return NULL;}
   }
   stock_id->chain_id_array = chain_id_array;

   return stock_id;
}


static short
CompareChainID(CandStructPtr new_cand, ChainIDPtr chain_id, short lmax, int num_t)
{
	short   ni, nd;
   int     i, temp, product;
   short  *ins_array, *del_array, *order, *t_vector;

   if(chain_id==NULL){printf("chain_id=NULL\n");getch(); return -3;}

   ni = new_cand->ni;
   nd = new_cand->nd;
   ins_array = new_cand->ins_array;
   del_array = new_cand->del_array;
   /*i_order*/
   temp = 0;
   product = 1;
   for(i=0;i<ni;i++){
   	temp += ins_array[i]*product;
      product *= (lmax+1);
   }
   if(chain_id->i_order!=temp){return 0;}
   /*d_order*/
   temp = 0;
   product = 1;
   for(i=0;i<nd;i++){
   	temp += del_array[i]*product;
      product *= (lmax+1);
   }
   if(chain_id->d_order!=temp){return 0;}

   /*order_index*/
   order = new_cand->order;
   temp=0;
   product = 1;
   for(i=0;i<(ni+nd);i++){
   	temp += order[i]*product;
      product *= 2;
   }
   if(chain_id->order_index!=temp){return 0;}
   /*t_index*/
   t_vector = new_cand->t_vector;
   temp = 0;
   product = 1;
   for(i=0;i<num_t;i++){
   	temp += t_vector[i]*product;
      product *= 2;
   }
   if(chain_id->t_index!=temp){return 0;}
   /*pos_vector*/
   if(chain_id->pos_vector==NULL){printf("Warning!pos_vector=NULL(press any key)\n");getch();return -1;}
	for(i=0;i<(ni+nd);i++){
   	if(*(chain_id->pos_vector+i)!=*(new_cand->pos_vector+i)){return 0;}
   }
	return 1;

}

long int
CheckRevisit(CandStructPtr new_cand, StockIDPtr stock_id, OptionsBlkPtr options)
{
	long int i;
   ChainIDPtr PNTR chain_id_array;
   short visited;

   chain_id_array = stock_id->chain_id_array;
	for(i=0;i<stock_id->current_id;i++){
		visited = CompareChainID(new_cand, chain_id_array[i], options->lmax, options->num_t);
      if(visited==-3){return visited;}
      if(visited){return i;}
   }
   return -1;
}

static double
CalculateTPOQijQjiRatio(CandStructPtr new_cand, CandStructPtr current_cand, PoolStructPtr pool)
{
	double tpo_ratio, t_ratio, p_ratio, o_ratio;
   int      i, nsum_i, nsum_j;


   nsum_i = current_cand->ni + current_cand->nd;
   nsum_j = new_cand->ni + new_cand->nd;

   t_ratio = 0.0;
   p_ratio = 0.0;
   if(nsum_j>nsum_i){
   	for(i=0;i<(nsum_j-nsum_i);i++){
      	t_ratio += log10((double)((pool->num_t)-(nsum_j-i)+1));
         t_ratio -= log10((double)(nsum_j-i));
      }
      for(i=0;i<(nsum_j-nsum_i);i++){
      	p_ratio += log10((double)(pool->orig_seqlen));
      }

   }
	if(nsum_j<nsum_i){
   	for(i=0;i<(nsum_i-nsum_j);i++){
      	t_ratio -= log10((double)((pool->num_t)-(nsum_i-i)+1));
         t_ratio += log10((double)(nsum_i-i));
      }
      for(i=0;i<(nsum_i-nsum_j);i++){
      	p_ratio -= log10((double)(pool->orig_seqlen));
      }
   }

   o_ratio=1.0;
   if((new_cand->ni+new_cand->nd)>0)for(i=(new_cand->ni+new_cand->nd);i>0;i--){o_ratio += log10((double)i);}
   if(new_cand->ni>0)for(i=new_cand->ni;i>0;i--){o_ratio -= log10((double)i);}
   if(new_cand->nd>0)for(i=new_cand->nd;i>0;i--){o_ratio -= log10((double)i);}
   if((current_cand->ni+new_cand->nd)>0)for(i=(current_cand->ni+current_cand->nd);i>0;i--){o_ratio -= log10((double)i);}
   if(current_cand->ni>0)for(i=current_cand->ni;i>0;i--){o_ratio += log10((double)i);}
   if(current_cand->nd>0)for(i=current_cand->nd;i>0;i--){o_ratio += log10((double)i);}

	tpo_ratio = 0.0;
   tpo_ratio += t_ratio;
   tpo_ratio += p_ratio;
   tpo_ratio += o_ratio;

   return tpo_ratio;
}


static double
CalculateConserveTPOQijQjiRatio(CandStructPtr new_cand, CandStructPtr current_cand, PoolStructPtr pool)
{
	double   tpo_ratio, t_ratio, p_ratio, o_ratio;
   int      i, nsum_i, nsum_j;
   short    *order;


   nsum_i = current_cand->ni + current_cand->nd -1;
   nsum_j = new_cand->ni + new_cand->nd -1;

   t_ratio = 0.0;
   p_ratio = 0.0;
   if(nsum_j>nsum_i){
   	for(i=0;i<(nsum_j-nsum_i);i++){
      	t_ratio += log10((double)((pool->num_t)-(nsum_j-i)+1));
         t_ratio -= log10((double)(nsum_j-i));
      }
      for(i=0;i<(nsum_j-nsum_i);i++){
      	p_ratio += log10((double)(pool->orig_seqlen));
      }
   }
	if(nsum_j<nsum_i){
   	for(i=0;i<(nsum_i-nsum_j);i++){
      	t_ratio -= log10((double)((pool->num_t)-(nsum_i-i)+1));
         t_ratio += log10((double)(nsum_i-i));
      }
      for(i=0;i<(nsum_i-nsum_j);i++){
      	p_ratio -= log10((double)(pool->orig_seqlen));
      }
   }

   o_ratio=0.0;
   order = new_cand->order;
   if(order[0]==1){
   	if((new_cand->ni+new_cand->nd-1)>0)for(i=(new_cand->ni+new_cand->nd-1);i>0;i--){o_ratio += log10((double)i);}
   	if((new_cand->ni-1)>0)for(i=(new_cand->ni-1);i>0;i--){o_ratio -= log10((double)i);}
   	if(new_cand->nd>0)for(i=new_cand->nd;i>0;i--){o_ratio -= log10((double)i);}
   	if((current_cand->ni+new_cand->nd-1)>0)for(i=(current_cand->ni+current_cand->nd-1);i>0;i--){o_ratio -= log10((double)i);}
   	if((current_cand->ni-1)>0)for(i=(current_cand->ni-1);i>0;i--){o_ratio += log10((double)i);}
   	if(current_cand->nd>0)for(i=current_cand->nd;i>0;i--){o_ratio += log10((double)i);}
   }
   else{
   	if((new_cand->ni+new_cand->nd-1)>0)for(i=(new_cand->ni+new_cand->nd-1);i>0;i--){o_ratio += log10((double)i);}
   	if(new_cand->ni>0)for(i=new_cand->ni;i>0;i--){o_ratio -= log10((double)i);}
   	if((new_cand->nd-1)>0)for(i=(new_cand->nd-1);i>0;i--){o_ratio -= log10((double)i);}
   	if((current_cand->ni+new_cand->nd-1)>0)for(i=(current_cand->ni+current_cand->nd-1);i>0;i--){o_ratio -= log10((double)i);}
   	if(current_cand->ni>0)for(i=current_cand->ni;i>0;i--){o_ratio += log10((double)i);}
   	if((current_cand->nd-1)>0)for(i=(current_cand->nd-1);i>0;i--){o_ratio += log10((double)i);}

   }
	tpo_ratio = 0.0;
   tpo_ratio += t_ratio;
   tpo_ratio += p_ratio;
   tpo_ratio += o_ratio;

   /*printf("tpo_ratio = %lf\n",tpo_ratio); */

   return tpo_ratio;
}


short
CompareCandidate(CandStructPtr new_cand, CandStructPtr current_cand, PoolStructPtr pool, PoolConsPtr pool_cons, short flag, FILE *fout)
{
   double   alpha_ij, q_ratio, chance;
   double   **rij;
	double    *num_cand;
   int      index_i, index_j;
   int      nsum_i, nsum_j;
   double   t_ratio, p_ratio, o_ratio;
   double  *limit_nind; /*probability of each (ni,nd) set selected*/
   double  *num_order;  /*number of candidates of order*/
   double  *num_t_cand; /*number of candidates of t_vector*/
	int     *num_cons;

   /*(0,0) redoing*/
   if(((new_cand->ni==0)&&(new_cand->nd==0))&&((current_cand->ni==0)&&(current_cand->nd==0))){
   	return 1;
   }

   alpha_ij = (new_cand->lnL)-(current_cand->lnL);
   /*printf("lnL(j)-lnL(i)=%5.1f\n",alpha_ij); */
   /*fprintf(fout,"lnL(j)-lnL(i)=%5.1f\n",alpha_ij);  */
   index_i = (current_cand->ni)*(pool->nmax+1)+(current_cand->nd);
   index_j = (new_cand->ni)*(pool->nmax+1)+(new_cand->nd);
   num_cand = pool->num_cand;

if(flag==0){
   if(new_cand->chain_no==4){
   	q_ratio = 0.0;/*the ratio of (ni,nd)i<->(ni,nd)j is 1.0*/
      num_cons = pool_cons->num_cons;
      /*printf("num_cons[%d]=%d, num_cons[%d]=%d\n",index_j, num_cons[index_j], index_i, num_cons[index_i]);*/
      if(num_cons[index_j]<=0){printf("num_cons[%d]<=0(%d)\n",index_j,num_cons[index_j]);getch();num_cons[index_j]=1;}
      if(num_cons[index_i]<=0){printf("num_cons[%d]<=0(%d)\n",index_i,num_cons[index_i]);getch();num_cons[index_i]=1;}
      q_ratio += log10((double)num_cons[index_j]);
      q_ratio -= log10((double)num_cons[index_i]);
      /*printf("q_ratio = %lf\n",q_ratio);getch(); */
      q_ratio += CalculateConserveTPOQijQjiRatio(new_cand, current_cand, pool);
      /*the end of modified part*/
      /*printf("ln(qji-qij)=%5.1f\n",q_ratio); */
      /*fprintf(fout,"ln(qji-qij)=%5.1f\n",q_ratio);*/
      alpha_ij += q_ratio;
   	/*printf("ln(alpha_ij)=%5.1f\n",alpha_ij);*/
   	/*fprintf(fout,"ln(alpha_ij)=%5.1f\n",alpha_ij); */
   }

   if(new_cand->chain_no==5){
   	rij = pool->rij;
      q_ratio = rij[index_j][index_i]/rij[index_i][index_j];
      if(num_cand[index_j]<0){printf("num_cand[%d]<0(%d)\n",index_j, num_cand[index_j]);getch();num_cand[index_j]=0.0;}
      if(num_cand[index_i]<0){printf("num_cand[%d]<0(%d)\n",index_i, num_cand[index_i]);getch();num_cand[index_i]=0.0;}
      q_ratio += num_cand[index_j];
      q_ratio -= num_cand[index_i];
      /*printf("q_ratio=%lf\n",q_ratio); */
      /*modified on Dec. 13, 2000*/
      /*t_ratio is the ratio of the number of t_vector candidates*/
      /*p_ratio is the ratio of the number of pos_vector candidates*/
      nsum_i = current_cand->ni + current_cand->nd;
      nsum_j = new_cand->ni + new_cand->nd;

  		if((pool->num_t_cand==NULL)||(new_cand->num_pos_cand<0)||(pool->num_order==NULL)){
         q_ratio += CalculateTPOQijQjiRatio(new_cand, current_cand, pool);
		}
      else{
      	num_t_cand = pool->num_t_cand;
      	num_order = pool->num_order;
  			t_ratio = num_t_cand[nsum_j]-num_t_cand[nsum_i];
         /*printf("t_ratio=%lf\n",t_ratio);
         printf("num_pos_cand=%lf:%lf\n",new_cand->num_pos_cand,current_cand->num_pos_cand);
         */
      	p_ratio = (new_cand->num_pos_cand)-(current_cand->num_pos_cand);
         /*printf("p_ratio=%lf\n",p_ratio); */
  			o_ratio = num_order[index_j]-num_order[index_i];
         /*printf("o_ratio=%lf\n",o_ratio); */
      	q_ratio += t_ratio;
      	q_ratio += p_ratio;
      	q_ratio += o_ratio;
  		}

      /*the end of modified part*/
      /*printf("ln(qji-qij)=%5.1f\n",q_ratio); */
      /*fprintf(fout,"ln(qji-qij)=%5.1f\n",q_ratio);*/
      alpha_ij += q_ratio;
   	/*printf("ln(alpha_ij)=%5.1f\n",alpha_ij); */
   	/*fprintf(fout,"ln(alpha_ij)=%5.1f\n",alpha_ij); */
   }
}
else{
   if(index_i!=index_j){
   	limit_nind = pool->limit_nind;
      q_ratio = (limit_nind[index_i+1]-limit_nind[index_i])/(limit_nind[index_j+1]-limit_nind[index_j]);
      if(q_ratio<=0){printf("q_ratio=%1.3e(press any key)\n");getch();q_ratio = MinlnL;}
      else{q_ratio = log10(q_ratio); }
      q_ratio += num_cand[index_j];
      q_ratio -= num_cand[index_i];
   }
   else{q_ratio=0.0;}
   /*modified on Dec. 13, 2000*/
   /*t_ratio is the ratio of the number of t_vector candidates*/
   /*p_ratio is the ratio of the number of pos_vector candidates*/
   nsum_i = current_cand->ni + current_cand->nd;
   nsum_j = new_cand->ni + new_cand->nd;

   if((pool->num_t_cand==NULL)||(new_cand->num_pos_cand==0)||(pool->num_order==NULL)){
   	q_ratio += CalculateTPOQijQjiRatio(new_cand, current_cand, pool);
   }
   else{
   	num_t_cand = pool->num_t_cand;
      num_order = pool->num_order;
      t_ratio = num_t_cand[nsum_j]-num_t_cand[nsum_i];
      p_ratio = (new_cand->num_pos_cand)-(current_cand->num_pos_cand);
      o_ratio = num_order[index_j]-num_order[index_i];
      q_ratio += t_ratio;
      q_ratio += p_ratio;
      q_ratio += o_ratio;
   }
   /*the end of modified part*/
   /*printf("ln(qji-qij)=%5.1f\n",q_ratio); */
   /*fprintf(fout,"ln(qji-qij)=%5.1f\n",q_ratio);*/
   alpha_ij += q_ratio;
   /*printf("ln(alpha_ij)=%5.1f\n",alpha_ij); */
   /*fprintf(fout,"ln(alpha_ij)=%5.1f\n",alpha_ij); */

}

   if(alpha_ij>=0){return 1;}
   else{
   	chance = rndu();
      if(chance < pow(10.0, alpha_ij)){return 1;}
      else{return 0;} 
   }

}

void
StockCandidate(StockIDPtr stock_id, CandStructPtr new_cand, CandStructPtr current_cand, PoolStructPtr pool, long int visited)
{
	long int current_id;
   ChainIDPtr PNTR chain_id_array;
   ChainIDPtr      chain_id;
	double temp;

   chain_id_array = stock_id->chain_id_array;
   if(visited>=0){
   	chain_id = chain_id_array[visited];
	   temp = 0;
   	temp += pow(10.0, chain_id->lnL)*(chain_id->num_visit);
   	temp += pow(10.0, new_cand->lnL);
   	temp /= (chain_id->num_visit+1);

      chain_id->num_visit += 1;
	   if(temp<=0){printf("temp(%lf)<=0",temp);chain_id->lnL=MinlnL;}
	   else{
   		chain_id->lnL= log10(temp);
   	}
   }
   else{
   	current_id = stock_id->current_id;
      /*chain_id_array[current_id] = ChainIDNew((short)(new_cand->ni+new_cand->nd));
      if(chain_id_array[current_id]==NULL){
      	printf("Cannot allocate chain_id_array[%ld](press any key)\n",current_id); getch();
         return;
      } */
      ArrayToID(new_cand, chain_id_array[current_id], pool->lmax, pool->nmax, pool->num_t);
      stock_id->current_id++;
   }
   CopyCandidate(current_cand, new_cand, pool->num_t);

   return;

}

void
DisplayStockID(StockIDPtr stock_id, OptionsBlkPtr options, FILE *fout)
{
	long int   i/*, current_id*/;
   ChainIDPtr PNTR chain_id_array;
   ChainIDPtr      chain_id;
   /*short      ni, nd, j;*/
	CandStructPtr cand;
   short         select_num;

   /*current_id = stock_id->current_id; */
   chain_id_array = stock_id->chain_id_array;
   cand = CandStructAllocate(options->nmax, options->num_t);
   if(cand==NULL){printf("ID to cand failed(press any key)\n"); getch(); return;}
   if((stock_id->select_num)>(stock_id->current_id)){select_num=(short)(stock_id->current_id);}
   else{select_num = (stock_id->select_num);}
   for(i=0;i<select_num;i++){
   	fprintf(fout,"No.%ld:",i);
      chain_id = chain_id_array[i];
		IDToArray(chain_id, cand, options->lmax, options->nmax, options->num_t);
      DisplayCandidate(cand, options->num_t, fout, 1);
      /*fprintf(fout,"nind_index=%d, ",chain_id->nind_index);
      fprintf(fout,"i_order=%d, d_order=%d, ",chain_id->i_order, chain_id->d_order);
      fprintf(fout,"order_index=%d, t_index=%d, \n",chain_id->order_index, chain_id->t_index);
      ni = (short)(chain_id->nind_index / (options->nmax+1));
      nd = (short)(chain_id->nind_index % (options->nmax+1));
   	for(j=0;j<(ni+nd);j++){
   		fprintf(fout,"%d,",*(chain_id->pos_vector+j));
   	}
      fprintf(fout,"\n");
      */
      fprintf(fout,"lnL=%10.4f, num_visit=%d\n",chain_id->lnL, chain_id->num_visit);

   }
   cand = CandDestruct(cand);

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


  	return;
}

void
CopyStockID(StockIDPtr stock_id, short nmax, long int id1, long int id2)
{
   ChainIDPtr PNTR chain_id_array;
   ChainIDPtr      chain_id1;
   ChainIDPtr      chain_id2;
   short      ni, nd;
   int        i;

   chain_id_array = stock_id->chain_id_array;
   chain_id1 = chain_id_array[id1];
   chain_id2 = chain_id_array[id2];

   /*copy chain_id2 to chain_id1*/
   chain_id1->nind_index=chain_id2->nind_index;
   chain_id1->i_order=chain_id2->i_order;
   chain_id1->d_order=chain_id2->d_order;
   chain_id1->order_index=chain_id2->order_index;
   chain_id1->t_index=chain_id2->t_index;
   ni = (short)(chain_id2->nind_index / (nmax+1));
   nd = (short)(chain_id2->nind_index % (nmax+1));
	for(i=0;i<(ni+nd);i++){
   	*(chain_id1->pos_vector+i)=*(chain_id2->pos_vector+i);
   }
	chain_id1->lnL=chain_id2->lnL;
   chain_id1->num_visit=chain_id2->num_visit;

   return;
}


void
SortStockID(StockIDPtr stock_id, short nmax)
{
   ChainIDPtr PNTR chain_id_array;
   ChainIDPtr      chain_id1;
   ChainIDPtr      chain_id2;
   ChainIDPtr      chain_temp;
   long int   i, j, k, current_id;
   short      ni, nd, p, select_num;
   double     prob, ave_prob, account, max_account;

   chain_id_array = stock_id->chain_id_array;
   current_id = stock_id->current_id;

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

   chain_temp = ChainIDNew(nmax);

   for(i=1;i<current_id;i++){
   	chain_id1 = chain_id_array[i];
   	for(j=0;j<i;j++){
      	chain_id2 = chain_id_array[j];
         if(chain_id1->num_visit>chain_id2->num_visit){
   			/*copy chain_id1 to chain_temp*/
   			chain_temp->nind_index=chain_id1->nind_index;
   			chain_temp->i_order=chain_id1->i_order;
   			chain_temp->d_order=chain_id1->d_order;
   			chain_temp->order_index=chain_id1->order_index;
   			chain_temp->t_index=chain_id1->t_index;
   			ni = (short)(chain_id1->nind_index / (nmax+1));
   			nd = (short)(chain_id1->nind_index % (nmax+1));
				for(p=0;p<(ni+nd);p++){
   				*(chain_temp->pos_vector+p)=*(chain_id1->pos_vector+p);
   			}
				chain_temp->lnL=chain_id1->lnL;
   			chain_temp->num_visit=chain_id1->num_visit;
            for(k=i-1;k>=j;k--){
					CopyStockID(stock_id, nmax, k+1, k);
            }
   			/*copy chain_temp to chain_id2*/
   			chain_id2->nind_index=chain_temp->nind_index;
   			chain_id2->i_order=chain_temp->i_order;
   			chain_id2->d_order=chain_temp->d_order;
   			chain_id2->order_index=chain_temp->order_index;
   			chain_id2->t_index=chain_temp->t_index;
   			ni = (short)(chain_temp->nind_index / (nmax+1));
   			nd = (short)(chain_temp->nind_index % (nmax+1));
				for(p=0;p<(ni+nd);p++){
   				*(chain_id2->pos_vector+p)=*(chain_temp->pos_vector+p);
   			}
				chain_id2->lnL=chain_temp->lnL;
   			chain_id2->num_visit=chain_temp->num_visit;

         }
      }
   }

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

   /*Calculate probability of A->S1*/
   if((stock_id->select_num)>current_id){
   	printf("select_num(%d)>current_id(%d)\n",(stock_id->select_num),current_id);
      select_num = (short)current_id;
   }
   else{select_num=stock_id->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_id->total_accepted)*0.90;
   ave_prob = 0;
   for(i=0;i<select_num;i++){
      chain_id1 = chain_id_array[i];
      prob = 1;
      prob /= chain_id1->num_visit;
      prob *= stock_id->total_accepted;
      prob *= pow(10.0, chain_id1->lnL);
      ave_prob += prob;
      if(prob<=0){printf("prob(%lf)<=0\n",prob);prob = MinlnL;}
      else{prob = log10(prob);}
      chain_id1->lnProbAS = prob;
      account += (chain_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_id->select_num = select_num;
   if(ave_prob<=0){printf("ave_prob(%lf)<=0\n",ave_prob);stock_id->Ave_lnAS = MinlnL;}
   else{stock_id->Ave_lnAS = log10(ave_prob);}

   chain_temp = ChainIDFree(chain_temp);

	return;
}

void
InitializeStockID(StockIDPtr stock_id, short nsum)
{
   ChainIDPtr PNTR chain_id_array;
   ChainIDPtr      chain_id;
   short           *pos_vector;
   int        i, j;

   chain_id_array = stock_id->chain_id_array;
   for(i=0;i<stock_id->current_id;i++){
   	chain_id = chain_id_array[i];
      chain_id->i_order=0;
      chain_id->d_order=0;
      chain_id->order_index=0;
      chain_id->t_index=0;
      chain_id->num_visit=0;
      chain_id->lnL=0;
      chain_id->lnProbAS=0;
      pos_vector = chain_id->pos_vector;
      for(j=0;j<nsum;j++){pos_vector[j]=-1;}
   }
   stock_id->current_id=0;
   stock_id->select_num = stock_id->select_id0;

   return;
}

void
StatisticalResult(StockIDPtr stock_id, int max_setnum, short nmax, FILE *fout)
{
	long int    i;
   int         nind_index, ni, nd;
   ChainIDPtr PNTR chain_id_array;
   ChainIDPtr      chain_id;
   double    *num_major, *max_lnL, *min_lnL;

   num_major = (double *)calloc(max_setnum,sizeof(double));
   if(num_major==NULL){printf("No memory for num_major\n");getch();return;}
   max_lnL = (double *)calloc(max_setnum,sizeof(double));
   if(max_lnL==NULL){printf("No memory for max_lnL\n");getch();return;}
   for(i=0;i<max_setnum;i++){
   	max_lnL[i]= MinlnL;
   }
   min_lnL = (double *)calloc(max_setnum,sizeof(double));
   if(min_lnL==NULL){printf("No memory for min_lnL\n");getch();return;}

   stock_id->total_accepted=0;
   chain_id_array = stock_id->chain_id_array;
   for(i=0;i<stock_id->current_id;i++){
   	chain_id = chain_id_array[i];
      nind_index = chain_id->nind_index;
      num_major[nind_index]+=chain_id->num_visit;
      stock_id->total_accepted += chain_id->num_visit;
      if(max_lnL[nind_index]<chain_id->lnL){max_lnL[nind_index]=chain_id->lnL;}
      if(min_lnL[nind_index]>chain_id->lnL){min_lnL[nind_index]=chain_id->lnL;}
   }
   fprintf(fout,"total_accepted=%ld\n",stock_id->total_accepted);
   for(i=0;i<max_setnum;i++){
   	ni = i/(nmax+1);
      nd = i%(nmax+1);
      if(num_major[i]>0){
      	num_major[i] /= stock_id->total_accepted;
      	fprintf(fout,"(%d,%d):num = %1.3e, max_lnL = %lf, min_lnL = %lf\n",ni, nd, num_major[i], max_lnL[i], min_lnL[i]);
   	}
   }
   return;
}

short
CheckMonopolyOfStockID(StockIDPtr stock_id, FILE *fout)
{
   ChainIDPtr PNTR chain_id_array;
   ChainIDPtr      chain_id;
   int             total_accepted, i;
   double          ratio;

   if(stock_id->current_id<=1){printf("only one chain-t-pos is accepted\n");return 1;}

   total_accepted=0;
   chain_id_array = stock_id->chain_id_array;
   for(i=0;i<stock_id->current_id;i++){
   	chain_id = chain_id_array[i];
      total_accepted += chain_id->num_visit;
   }
   /*check monopoly*/
   for(i=0;i<stock_id->current_id;i++){
   	chain_id = chain_id_array[i];
      ratio = (double)(chain_id->num_visit)/total_accepted;
      if(ratio>=0.99){printf("the destribution of that chain-t-pos is more than 0.99\n");return 1;}
   }
   return 0;

}

/*Add on March 22,2002*/
short
CompareID(ChainIDPtr chain_id1, ChainIDPtr chain_id2, short nmax)
{
	int i, ni, nd;
   short *pos1, *pos2;

   if(chain_id1->nind_index!=chain_id2->nind_index){return 0;}
   if(chain_id1->i_order!=chain_id2->i_order){return 0;}
   if(chain_id1->d_order!=chain_id2->d_order){return 0;}
   if(chain_id1->order_index!=chain_id2->order_index){return 0;}
   if(chain_id1->t_index!=chain_id2->t_index){return 0;}
   pos1 = chain_id1->pos_vector;
   pos2 = chain_id2->pos_vector;
   ni = (short)(chain_id1->nind_index / (nmax+1));
   nd = (short)(chain_id1->nind_index % (nmax+1));
	for(i=0;i<(ni+nd);i++){
   	if(pos1[i]!=pos2[i]){return 0;}
   }

   return 1;
}
void
AddID(ChainIDPtr chain_id2, ChainIDPtr chain_id1)
{
	double temp;

   temp = 0;
   temp += pow(10.0, chain_id1->lnL)*(chain_id1->num_visit);
   temp += pow(10.0, chain_id2->lnL)*(chain_id2->num_visit);
   chain_id2->num_visit += chain_id1->num_visit;
   temp /= (chain_id2->num_visit);
   if(temp<=0){printf("temp(%lf)<=0",temp);chain_id2->lnL=MinlnL;}
   else{
   	chain_id2->lnL= log10(temp);
   }

	return;
}
void
NewID(ChainIDPtr chain_id2, ChainIDPtr chain_id1, short nmax)
{
   int   i, ni, nd;

   /*copy chain_id1 to chain_id2*/
   chain_id2->nind_index=chain_id1->nind_index;
   chain_id2->i_order=chain_id1->i_order;
   chain_id2->d_order=chain_id1->d_order;
   chain_id2->order_index=chain_id1->order_index;
   chain_id2->t_index=chain_id1->t_index;
   ni = (short)(chain_id1->nind_index / (nmax+1));
   nd = (short)(chain_id1->nind_index % (nmax+1));
	for(i=0;i<(ni+nd);i++){
   	*(chain_id2->pos_vector+i)=*(chain_id1->pos_vector+i);
   }
	chain_id2->lnL=chain_id1->lnL;
   chain_id2->num_visit=chain_id1->num_visit;

   return;
}
void
AddEachIDtoStockID(StockIDPtr stock_id, StockIDPtr each_id, short nmax)
{
   long int   i, j, p, current_id1, current_id2;
   ChainIDPtr PNTR chain_id_array1;
   ChainIDPtr PNTR chain_id_array2;
   ChainIDPtr      chain_id1;
   ChainIDPtr      chain_id2;

   current_id1 = each_id->current_id;
   chain_id_array1 = each_id->chain_id_array;
   current_id2 = stock_id->current_id;
   chain_id_array2 = stock_id->chain_id_array;
   for(i=0;i<current_id1;i++){
   	chain_id1 = chain_id_array1[i];
      p=-1;
      for(j=0;j<current_id2;j++){
      	chain_id2 = chain_id_array2[j];
         if(CompareID(chain_id1,chain_id2,nmax)==1){
            p=j;  
         	AddID(chain_id2, chain_id1);break;
         }
      }
      if(p<0){
      	chain_id2 = chain_id_array2[stock_id->current_id];
         NewID(chain_id2, chain_id1, nmax);
         stock_id->current_id++;
      }
   }

   return;
}