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

static void
CalculateTindex(short *t_clump, short num_size, int num_t);

static short
CompareChainPos(ChainIDPtr chain_id, ChainIDPtr sum_id, short nmax);
static void
CopyChainPos(ChainIDPtr dist_id, ChainIDPtr src_id, short nmax, short *t_clump, short num_size);
static void
SumChainPos(ChainIDPtr sum_id, ChainIDPtr chain_id);

/*clump t_index*/
void
CalculateTindex(short *t_clump, short num_size, int num_t)
{
	int i, j, k, p, num, temp, product;
   int dev, remind;
   short *t_vector;

   t_clump[0]=0;
   t_vector = (short *)calloc(num_t,sizeof(short));
   for(num=1;num<=num_size;num++){
      dev = num_t/(num+1);
      remind = num_t%(num+1);
      for(i=0;i<num_t;i++){t_vector[i]=0;}
      for(i=0;i<num;i++){
      	j=i/2+1; k=i%2;
         if(!k){p=j*dev-1;if((remind)&&(j>1)){p++;remind--;}}/*even*/
         else{p=num_t-j*dev-1;if((remind)&&(j>1)){p--;remind--;}}/*odd*/
         t_vector[p]=1;
      }
   	temp = 0;
   	product = 1;
   	for(i=0;i<num_t;i++){
   		temp += t_vector[i]*product;
      	product *= 2;
   	}
      t_clump[num]=(short)temp;
   }
   free(t_vector);

   return;
}

/*compare chain and position(t_vector can be different)*/
static short
CompareChainPos(ChainIDPtr chain_id, ChainIDPtr sum_id, short nmax)
{
   int   i, nsum;
   short *pos_vector1, *pos_vector2;

   if(chain_id->nind_index!=sum_id->nind_index) return -1;
   if(chain_id->i_order!=sum_id->i_order) return -1;
   if(chain_id->d_order!=sum_id->d_order) return -1;
   if(chain_id->order_index!=sum_id->order_index) return -1;

   pos_vector1 = chain_id->pos_vector;
   pos_vector2 = sum_id->pos_vector;
   nsum = chain_id->nind_index/(nmax+1)+chain_id->nind_index%(nmax+1);
   for(i=0;i<nsum;i++){
   	if(pos_vector1[i]!=pos_vector2[i]){return -1;}
   }

   return 1;
}

static void
CopyChainPos(ChainIDPtr dist_id, ChainIDPtr src_id, short nmax, short *t_clump, short num_size)
{
   int   i, nsum;
   short *dist_pos, *src_pos;

   dist_id->nind_index = src_id->nind_index;
   dist_id->i_order = src_id->i_order;
   dist_id->d_order = src_id->d_order;
   dist_id->order_index = src_id->order_index;
   dist_pos = dist_id->pos_vector;
   src_pos = src_id->pos_vector;
   nsum = src_id->nind_index/(nmax+1)+src_id->nind_index%(nmax+1);
   for(i=0;i<nsum;i++){dist_pos[i]=src_pos[i];}
   dist_id->lnL = src_id->lnL;
   dist_id->num_visit = src_id->num_visit;

   /*clump t_vector*/
   if(nsum>num_size){printf("Warning!nsum(%d)>num_size(%d)\n",nsum,num_size);getch();dist_id->t_index=0;}
   else{dist_id->t_index = t_clump[nsum];}

   return;
}

static void
SumChainPos(ChainIDPtr sum_id, ChainIDPtr chain_id)
{
	double temp;

   temp = 0;
   /*temp += pow(10.0, sum_id->lnL)*(sum_id->num_visit);
   temp += pow(10.0, chain_id->lnL)*(chain_id->num_visit);
   temp /= (sum_id->num_visit + chain_id->num_visit); */
   /*Modified on Feb 16,2001*/
   temp += pow(10.0, sum_id->lnL);
   temp += pow(10.0, chain_id->lnL);
   if(temp<=0){printf("temp=%1.3e(press any key)\n");/*getch();*/sum_id->lnL = MinlnL;}
   else{sum_id->lnL = log10(temp); }

   sum_id->num_visit += chain_id->num_visit;

   return;

}



/*Summarize the "chain-t-pos"es with different t-vector and the same chain-pos*/

void
SummarizedStockID(StockIDPtr summarized_id, StockIDPtr stock_id, short nmax, int num_t)
{
   long int i,j,p,max_idnum2,found;
   ChainIDPtr PNTR chain_id_array;
   ChainIDPtr PNTR sum_id_array;
   ChainIDPtr      chain_id;
   short      *t_clump;

   max_idnum2 = summarized_id->max_idnum;
   /*calculate clump t_index*/
   t_clump = (short *)calloc(nmax+1,sizeof(short));
   if(t_clump==NULL){printf("No memory for t_clump\n");getch();return;}
   CalculateTindex(t_clump, nmax, num_t);

   chain_id_array = stock_id->chain_id_array;
   sum_id_array = summarized_id->chain_id_array;
   p=0;
   for(i=0;i<stock_id->current_id;i++){
   	chain_id = chain_id_array[i];
      if(p==0){found=-1;}
      else{
         found = -1;
      	for(j=0;j<p;j++){
         	found = CompareChainPos(chain_id, sum_id_array[j], nmax);
         	if(found>=0){found = j; break;}
         }
      }
      if(found<0){
      	CopyChainPos(sum_id_array[p], chain_id, nmax, t_clump, nmax);
         p++;
         if(p>=max_idnum2){
         	printf("Warning!currrent_id(%d)>max_idnum(%d)\n",p,max_idnum2);getch();
            summarized_id->current_id = p; return;
         }
      }
      else{
      	SumChainPos(sum_id_array[found], chain_id);
      }
   }
   summarized_id->current_id = p;

   return;
}