/* This program performs proposing of new (c,t,p), thereby conducting the
MCMC-cycles.*/  

#include "pseud.h"

/*functions which are used in the loop of candidate ins_array, del_array*/
static void
CountNumCand(short ni, short nd, short lmax, short nmax, int target_d, PoolStructPtr pool);

/*select vector for each (ni,nd)*/
static short
SelectOrder(short *order, short ni, short nd, short flag);
static short
SelectTvector(short *t_vector, int num_t, short nsum, short flag);

static void
RandomSelectPosVector(CandStructPtr new_cand, int seqlen, int num_t);
static void
SelectNeighbourPos(CandStructPtr new_cand, int seqlen, int num_t, int no);
static void
SelectPosVector(CandStructPtr new_cand, int seqlen, int num_t, short flag);
static void
SelectIndelArray(PoolStructPtr pool, CandStructPtr new_cand, short flag);
static short
SearchCandidate(int select, short lmax, int target_d, short ni, short nd, short *temp_ins, short *temp_del);
static void
SearchIndelArray(PoolStructPtr pool, CandStructPtr new_cand, short flag);

static short
CheckNewCandidate(CandStructPtr new_cand, CandStructPtr current_cand, int num_t);
static short
CountNumPosCand(CandStructPtr new_cand, int seqlen, int num_t);


/* This function ensures that the given (c,t,p) generate a sequence with the proper length */
int
CheckLengthChange(short *ins_array, short ni, short *del_array, short nd)
{ 
   int  i, d;

   d=0;
	for(i=0;i<ni;i++){d += ins_array[i];}
	for(i=0;i<nd;i++){d -= del_array[i];}
   return d;

}

void
CheckCarryUp(short *array, int n, int lmax)
{
	int i;

   for(i=n-1;i>0;i--){
   	if(array[i]<=lmax){break;}
      else{array[i-1]++;  array[i]=1;}
   }
	return;
}


/* De-allocates the memory used for PoolStruct */ 
PoolStructPtr
PoolDestruct(PoolStructPtr pool)
{
	short i, max_setnum;
   double **rij;
  	int   **i_pool, **d_pool;

   if(pool==NULL) return NULL;
   max_setnum = pool->max_setnum;
   if(pool->rij){
      rij=pool->rij;
   	for(i=0;i<max_setnum;i++){rij[i]=MemFree(rij[i]);}
      rij = MemFree(rij);
   }
   if(pool->num_cand){pool->num_cand = MemFree(pool->num_cand);}
   if(pool->limit_nind){pool->limit_nind = MemFree(pool->limit_nind);}
   if(pool->i_pool){
   	i_pool = pool->i_pool;
      for(i=0;i<max_setnum;i++){
      	i_pool[i] = MemFree(i_pool[i]);
      }
      i_pool = MemFree(i_pool);
   }
   if(pool->d_pool){
   	d_pool = pool->d_pool;
      for(i=0;i<max_setnum;i++){
      	d_pool[i] = MemFree(d_pool[i]);
      }
      d_pool = MemFree(d_pool);
   }
   if(pool->num_order){pool->num_order = MemFree(pool->num_order);}
   if(pool->num_t_cand){pool->num_t_cand = MemFree(pool->num_t_cand);}
   pool = MemFree(pool);

   return pool;

}


/* The number of (c,t,p)s belonging to each (ni,nd) group is counted by this function */
static void
CountNumCand(short ni, short nd, short lmax, short nmax, int target_d, PoolStructPtr pool)
{
	int   i, d, nind_index,product;
   long int   num_passed;
   int         max_iter, iter;
	short *ins_array, *del_array;
   		/*ins_array is the array of the length of insertions
           del_array is the array of the length of deletions*/
   double  *num_cand;
	int   i_order, d_order;
  	int   **i_pool, **d_pool;

   if((ni==0)&&(nd==0)){
   	nind_index = (nmax+1)*ni + nd;
   	i_pool = pool->i_pool;
   	d_pool = pool->d_pool;
   	num_cand = pool->num_cand;
   	num_cand[nind_index]=0.0;
   	i_pool[nind_index]=(int *)calloc(1,sizeof(int));
   	if(i_pool==NULL){printf("out of memory:i_pool(press any key)\n");getch();}
   	d_pool[nind_index]=(int *)calloc(1,sizeof(int));
   	if(d_pool==NULL){printf("out of memory:d_pool(press any key)\n");getch();}
		i_pool[nind_index][0]=0;
		d_pool[nind_index][0]=0;
      return;
   }


   ins_array = (short *)calloc(nmax,sizeof(short));
   del_array = (short *)calloc(nmax,sizeof(short));
   if(del_array==NULL){printf("out of memory(press any key)\n");getch();}

   num_passed=0;
   max_iter = 1;
   for(i=0;i<ni;i++){max_iter *= lmax;}
   for(i=0;i<nd;i++){max_iter *= lmax;}
   /*printf("max_iter=%d\n",max_iter);  */

   nind_index = (nmax+1)*ni + nd;
   i_pool = pool->i_pool;
   d_pool = pool->d_pool;
   i_pool[nind_index]=(int *)calloc(pool->max_stored,sizeof(int));
   if(i_pool==NULL){printf("out of memory:i_pool(press any key)\n");getch();}
   d_pool[nind_index]=(int *)calloc(pool->max_stored,sizeof(int));
   if(d_pool==NULL){printf("out of memory:d_pool(press any key)\n");getch();}

   /*initialization*/
   for(i=0;i<ni;i++){ins_array[i]=1;}
   for(i=0;i<nd;i++){del_array[i]=1;}
   iter=0;
   while(iter<max_iter){
      d = CheckLengthChange(ins_array, ni, del_array, nd);
      if(d==target_d){
         if(num_passed<(pool->max_stored)){
         	i_order = 0;
            product = 1;
		   	for(i=0;i<ni;i++){
		      	i_order += ins_array[i]*product;
               product *= (lmax+1);
            }
         	i_pool[nind_index][num_passed]=i_order;
         	d_order = 0;
            product = 1;
		   	for(i=0;i<nd;i++){
		      	d_order += del_array[i]*product;
               product *= (lmax+1);
   			}
         	d_pool[nind_index][num_passed]=d_order;
         }
      	num_passed++;
         if(num_passed>=(pool->max_stored)){
				i_pool[nind_index]=MemFree(i_pool[nind_index]);
				d_pool[nind_index]=MemFree(d_pool[nind_index]);
         }
      }
      /*loop of del*/
      if(nd>0){
      	del_array[nd-1]++;
      	CheckCarryUp(del_array, nd, lmax);
      }
      if((del_array[0]>lmax)||(nd==0)){
      	if(ni>0){
      		ins_array[ni-1]++;
      		CheckCarryUp(ins_array, ni, lmax);
   			for(i=0;i<nd;i++){del_array[i]=1;}
         }
         if((ins_array[0]>lmax)||(ni==0)){break;}
      }
      iter++;
   }
   /*printf("total_checked=%d\n",iter);  */
   num_cand = pool->num_cand;
   if(num_passed==0){
      num_cand[nind_index]=-1.0;
   	realloc(i_pool[nind_index],1*sizeof(int)); i_pool[nind_index][0]=0;
   	realloc(d_pool[nind_index],1*sizeof(int)); d_pool[nind_index][0]=0;
   }
   else{
      num_cand[nind_index]=log10((double)num_passed);
      if(num_passed<(pool->max_stored)){
   		realloc(i_pool[nind_index],num_passed*sizeof(int));
   		realloc(d_pool[nind_index],num_passed*sizeof(int));
      }
   }
   free(ins_array); free(del_array);

   return;

}


/* PoolStruct  is a look-up table of (c,t,p) candidates */
PoolStructPtr
SetupPoolStruct(short lmax, short nmax, int target_d, int orig_seqlen, int num_t)
{
   PoolStructPtr pool;
	short i, j, ni, nd, max_setnum;
   double **rij;
   double *num_cand;
   double *limit_nind;
   int   **i_pool, **d_pool;
   double  *num_order, *num_t_cand;
   int   max_d, min_d, nind_index;
   short i_ni, j_ni, i_nd, j_nd, dist_ni, dist_nd, dist;

   pool = (PoolStructPtr)calloc(1,sizeof(PoolStruct));
   if(pool==NULL){return NULL;}

	pool->lmax = lmax;
   pool->nmax = nmax;
   pool->target_d = target_d;
   pool->num_t = num_t;
   pool->orig_seqlen = orig_seqlen;

   max_setnum = (short)((nmax+1)*(nmax+1));
   pool->max_setnum = max_setnum;
   rij = (double **)calloc(max_setnum,sizeof(double *));
   if(rij==NULL){pool = MemFree(pool);return NULL;}
   for(i=0;i<max_setnum;i++){
   	rij[i] = (double *)calloc(max_setnum,sizeof(double));
      if(rij[i]==NULL){pool = MemFree(pool);return NULL;}
   }
   pool->rij = rij;

   num_cand = (double *)calloc(max_setnum,sizeof(double));
   if(num_cand==NULL){pool = PoolDestruct(pool);return NULL;}

   pool->num_cand = num_cand;

   i_pool = (int **)calloc(max_setnum,sizeof(int *));
   if(i_pool==NULL){pool = PoolDestruct(pool); return NULL;}
   d_pool = (int **)calloc(max_setnum,sizeof(int *));
   if(d_pool==NULL){pool = PoolDestruct(pool); return NULL;}
   pool->i_pool = i_pool;
   pool->d_pool = d_pool;
	pool->max_stored = 500;

   /*count candidate*/
   for(ni=0;ni<=nmax;ni++){
      for(nd=0;nd<=nmax;nd++){
       if((ni+nd)>nmax){
            nind_index = (nmax+1)*ni+nd;
            num_cand[nind_index]=-1.0;
            /*printf("num_cand[%d]=%d\n",nind_index,num_cand[nind_index]); */
   			i_pool[nind_index]=(int *)calloc(1,sizeof(int));
   			i_pool[nind_index][0]=0;
   			d_pool[nind_index]=(int *)calloc(1,sizeof(int));
   			d_pool[nind_index][0]=0;
       }
       else{
         /*check minimum requirement*/
         min_d = ni*1-nd*lmax;
         max_d = ni*lmax-nd*1;
         if((min_d<=target_d)&&(max_d>=target_d)){
            CountNumCand(ni, nd, lmax, nmax, target_d, pool);
         }
         else{
            nind_index = (nmax+1)*ni+nd;
            num_cand[nind_index]=-1.0;
            /*printf("num_cand[%d]=%d\n",nind_index,num_cand[nind_index]); */
   			i_pool[nind_index]=(int *)calloc(1,sizeof(int));
   			i_pool[nind_index][0]=0;
   			d_pool[nind_index]=(int *)calloc(1,sizeof(int));
   			d_pool[nind_index][0]=0;
         }
       }
      }
   }
   /*calculate rij*/
   for(i=0;i<max_setnum;i++){
   	i_ni = (short)(i/(nmax+1));
   	i_nd = (short)(i%(nmax+1));
   	for(j=0;j<max_setnum;j++){
   		j_ni = (short)(j/(nmax+1));
   		j_nd = (short)(j%(nmax+1));
         dist_ni = (short)(j_ni-i_ni);
         dist_nd = (short)(j_nd-i_nd);
         dist = (short)(dist_ni*dist_ni + dist_nd*dist_nd);
         if((j_ni+j_nd)>nmax){rij[i][j]=0.0;}
         if(dist==0){rij[i][j]=0.0;}
         else{rij[i][j]=100.0/dist;}
      }
   }
   /*setup limit_nind*/
   limit_nind = (double *)calloc(max_setnum+1,sizeof(double));
   if(limit_nind==NULL){printf("no memory for limit_nind\n");getch();return pool;}
   pool->limit_nind = limit_nind;
   /*the case of even limit_nind*/
   limit_nind[0]=0.0;
   for(i=0;i<max_setnum;i++){
      if(num_cand[i]>=0){limit_nind[i+1]=limit_nind[i]+1.0; }
      else{limit_nind[i+1]=limit_nind[i];}
   }
   /*the case of uneven limit_nind*/
   /*limit_nind[0]=0.0;
   for(i=0;i<max_setnum;i++){
      limit_nind[i+1]=limit_nind[i]+num_cand[i];
   }
   */
   /*count num_order*/
   pool->num_order = (double *)calloc(max_setnum,sizeof(double));
   if(pool->num_order==NULL){
   	printf("No memory for num_order\n");getch();  return pool;
   }
   num_order = pool->num_order;
   for(i=0;i<max_setnum;i++){
   	ni = (short)(i / (nmax+1));
      nd = (short)(i % (nmax+1));
      /*(ni+nd)!/(ni!*nd!)*/
      num_order[i]=0.0;
      for(j=(ni+nd);j>0;j--){num_order[i] += log10((double)j);}
      for(j=ni;j>0;j--){num_order[i] -= log10((double)j);}
      for(j=nd;j>0;j--){num_order[i] -= log10((double)j);}
   }
   /*count num_t_cand---num_tCnsum*/
   pool->num_t_cand = (double *)calloc(nmax+1, sizeof(double));
   if(pool->num_t_cand==NULL){
   	printf("No memory for num_t_cand\n");getch();  return pool;
   }
   num_t_cand = pool->num_t_cand;
   num_t_cand[0]=1;
   for(i=1;i<=nmax;i++){
   	num_t_cand[i]=0.0;
      for(j=0;j<i;j++){num_t_cand[i]+=log10((double)(pool->num_t - j));}
      for(j=i;j>0;j--){num_t_cand[i]-=log10((double)j);}
   }

   return pool;

}

void
DisplayPoolStruct(PoolStructPtr pool, FILE *fout)
{
   int i, ni, nd;
	double  *num_cand;  /*log10(number of candidate of the set(ni, nd))*/
   /*int    j, num_passed, **i_pool, **d_pool; */

   fprintf(fout,"pool\n");
   fprintf(fout,"lmax=%d, nmax=%d\n",pool->lmax, pool->nmax);
   fprintf(fout,"max_stored=%d\n",pool->max_stored);
   num_cand=pool->num_cand;
   /*i_pool = pool->i_pool;
   d_pool = pool->d_pool; */
   for(i=0;i<pool->max_setnum;i++){
      ni=i/(pool->nmax+1);
      nd=i%(pool->nmax+1);
   	fprintf(fout,"(ni,nd)=(%d,%d), ",ni,nd);
      fprintf(fout,"num_cand[]=%d\n",round(pow(10.0,num_cand[i])));
      /*num_passed = round(pow(10.0,num_cand[i]));
      if((num_passed>=0)&&(num_passed<(pool->max_stored))){
      	for(j=0;j<num_passed;j++){
         	fprintf(fout,"(%d,%d)\n",i_pool[i][j],d_pool[i][j]);
         }
      }
      */
   }

   return;
}


/* Once the number (and lengths) of insertions and deletions are decided, this function
arranges the order of these events */
static short
SelectOrder(short *order, short ni, short nd, short flag)
{
	int  i, nsum;
   short  num_i, num_d;
   short  ok_flag;
   int    r1, r2;
   int    order_index, new_order, product, max_order, num_order, temp;
   int   *order_cand;

  if(flag<=0){   /*random selection*/
  	if(ni*nd==0){
   	if(flag==-1){
      	if(ni==0){for(i=0;i<nd;i++)order[i]=0;}
      	if(nd==0){for(i=0;i<ni;i++)order[i]=1;}
      	return 0;
      }
      else{return 1;}
   }
   order_index=0;
   product = 1;
   for(i=0;i<(ni+nd);i++){
   	order_index += order[i]*product;
      product *= 2;
   }
   nsum = ni + nd;
   max_order = 1;
   for(i=0;i<nsum;i++){max_order *= 2;}
   order_cand = (int *)calloc(max_order,sizeof(int));
   num_order=0;
   for(new_order=0;new_order<max_order;new_order++){
      num_i=0;
      temp = new_order;
   	for(i=0;i<nsum;i++){
      	if(temp%2){num_i++;}
         temp /= 2;
      }
      if(num_i==ni){order_cand[num_order]=new_order;num_order++;}
   }
   /*printf("num_order=%d(press any key)\n",num_order); getch(); */

   if((num_order==1)&&(flag==0)){order_cand = MemFree(order_cand); return 1;}

   ok_flag = 1;
   while(ok_flag){
   	new_order = floor(rndu()*num_order);
      new_order = order_cand[new_order];
      if((order_index!=new_order)||(flag==-1)){
   		for(i=0;i<nsum;i++){
         	order[i]=(short)(new_order % 2);
            new_order /= 2;
         }
      	ok_flag = 0;
      }
   }
   order_cand = MemFree(order_cand);
   return ok_flag;
  }
  if(flag==1){/*select neighboring vector*/
  	if(ni*nd==0){/*printf("cannot change order\n");getch();*/return 1;}
   else{
   	r1 = ceil(rndu()*ni);
      r2 = ceil(rndu()*ni);
   	nsum = ni + nd;
      num_i=0; num_d=0;
      for(i=0;i<nsum;i++){
      	if(order[i]==1){num_i++;if(num_i==r1){r1=i;}}
      	if(order[i]==0){num_d++;if(num_d==r2){r2=i;}}
      }
      order[r1]=0;order[r2]=1;/*exchange*/
      return 0;
   }
  }
  return 1;

}
static short
SelectTvector(short *t_vector, int num_t, short nsum, short flag)
{
	int i,num, r1, r2, num_1, direction, count;
   short  ok_flag;
   int    t_index, new_t, product, temp;
   int    max_index, *index_cand, num_index;

  if(flag<=0){   /*random selection*/
	t_index = 0;
   product = 1;
   for(i=0;i<num_t;i++){
   	t_index += t_vector[i]*product;
      product *= 2;
   }

   max_index = 1;
   for(i=0;i<num_t;i++){max_index *= 2;}
   index_cand = (int *)calloc(max_index,sizeof(int));
   num_index=0;
   for(new_t=0;new_t<max_index;new_t++){
      num_1=0;
      temp = new_t;
   	for(i=0;i<num_t;i++){
      	if(temp%2){num_1++;}
         temp /= 2;
      }
      if(num_1==nsum){index_cand[num_index]=new_t;num_index++;}
   }
   /*printf("num_index=%d(press any key)\n",num_index); getch(); */

   if((num_index==1)&&(flag==0)){index_cand = MemFree(index_cand); return 1;}

   ok_flag = 1;
   while(ok_flag){
   	new_t = floor(rndu()*num_index);
      new_t = index_cand[new_t];
      if((t_index!=new_t)||(flag==-1)){
   		for(i=0;i<num_t;i++){
         	t_vector[i]=(short)(new_t % 2);
            new_t /= 2;
         }
      	ok_flag = 0;
      }
   }
   index_cand = MemFree(index_cand);
   return ok_flag;
  }
  /*select neighbour t-unit*/
  if(flag==1){
  	r1 = ceil(rndu()*nsum);
   num=0;
   for(i=0;i<num_t;i++){
   	if(t_vector[i]==1){num++;if(num==r1){r1=i;}}
   }
   if(r1==0){
   	r2=r1;
   	while(t_vector[r2]==1){
      	r2++;
      	if(r2==(num_t-1)){return 1;}
      }
      t_vector[r1]=0; t_vector[r2]=1;  return 0;
   }
   if(r1==(num_t-1)){
   	r2=r1;
      while(t_vector[r2]==1){
      	r2--;
      	if(r2==0){return 1;}
      }
      t_vector[r1]=0; t_vector[r2]=1;  return 0;
   }
   direction = ((rndu()-0.5)>=0? 1 : -1);
   r2=r1;
   count=0;
   while(t_vector[r2]==1){
   	r2+=direction;  count++;
   	if((r2==0)||(r2==(num_t-1))){r2=r1;direction*=-1;}
      if(count>=num_t){return 1;}
   }
   t_vector[r1]=0; t_vector[r2]=1;  return 0;
  }
  return 1;
}


/*  id_array is a array holding information of an insertion or a deletion occurring
for a given time period.  SetupIDarray summarize such information and store into id_array */
void
SetupIDarray(CandStructPtr cand, int num_t)
{
	int i,j,k,p;
   short  *ins_array, *del_array;
   short  *order;
   short  *id_array;
   short  *t_vector;

   ins_array = cand->ins_array;
   del_array = cand->del_array;
   order = cand->order;
   id_array = cand->id_array;
   t_vector = cand->t_vector;

   /*id_array*/
   j=0; k=0; p=0;
   for(i=0;i<num_t;i++){
   	if(t_vector[i]==1){
      	if(order[p]==1) {id_array[i]=ins_array[j]; j++;}
         else {id_array[i]=(short)(-del_array[k]); k++;}
         p++;
      }
      else{id_array[i]=0;}
   }

   return;

}


/* Specifies the positions of insertions and deletions by random number */
static void
RandomSelectPosVector(CandStructPtr new_cand, int seqlen, int num_t)
{
	int    p, t, current_seqlen;
   short  *id_array, *pos_vector;

   id_array = new_cand->id_array;
   pos_vector = new_cand->pos_vector;

   current_seqlen = seqlen;
   p=0;
   for(t=0;t<num_t;t++){
   	if(id_array[t]>0){
      	pos_vector[p]=(short)floor(rndu()*(current_seqlen+1)); p++;
         current_seqlen += id_array[t];
      }
      if(id_array[t]<0){
      	pos_vector[p]=(short)floor(rndu()*(current_seqlen+id_array[t]+1)); p++;
         current_seqlen += id_array[t];
      }
   }

	return;

}



/*  Upon proposal of a new (c,t,p), a (c,t,p) specifying a 'pos'b different but close to the 
old one is occasionally chosen using the following function */
static void
SelectNeighbourPos(CandStructPtr new_cand, int seqlen, int num_t, int no)
{
	int    p, t, current_seqlen, direction;
   short  *id_array, *pos_vector;

   id_array = new_cand->id_array;
   pos_vector = new_cand->pos_vector;

   current_seqlen = seqlen;
   p=0;
   for(t=0;t<num_t;t++){
   	if(id_array[t]>0){
         if(p==no){
         	if(pos_vector[p]==0){pos_vector[p]=1; return;}
         	if(pos_vector[p]==current_seqlen){pos_vector[p]=(short)(current_seqlen-1); return;}
   			direction = ((rndu()-0.5)>=0? 1 : -1);
            pos_vector[p] = (short)(pos_vector[p]+direction); return;
         }
         p++;
         current_seqlen += id_array[t];
      }
      if(id_array[t]<0){
         if(p==no){
         	if(pos_vector[p]==0){pos_vector[p]=1; return;}
         	if(pos_vector[p]==(current_seqlen+id_array[t])){pos_vector[p]=(short)(current_seqlen+id_array[t]-1); return;}
   			direction = ((rndu()-0.5)>=0? 1 : -1);
            pos_vector[p] = (short)(pos_vector[p]+direction); return;
         }
      	p++;
         current_seqlen += id_array[t];
      }
   }

	return;

}

static void
SelectPosVector(CandStructPtr new_cand, int seqlen, int num_t, short flag)
{
   int    i, p, nsum, no, current_seqlen;
   short  *pos_vector, *old_pos;
   short  ok_flag;
   short  *id_array;

   if(flag==2){/*check pos is in the limit*/
		pos_vector = new_cand->pos_vector;
      id_array = new_cand->id_array;
      current_seqlen = seqlen;
      p=0;
      for(i=0;i<num_t;i++){
      	if(id_array[i]>0){
         	if(pos_vector[p]<0){pos_vector[p]=0;}
            if(pos_vector[p]>current_seqlen){pos_vector[p]=(short)current_seqlen;}
            current_seqlen += id_array[i];
            p++;
         }
      	if(id_array[i]<0){
         	if(pos_vector[p]<0){pos_vector[p]=0;}
            if(pos_vector[p]>(current_seqlen+id_array[i])){pos_vector[p]=(short)(current_seqlen+id_array[i]);}
            current_seqlen += id_array[i];
            p++;
         }
      }
      return;
   }

   if(flag==1){
      nsum=new_cand->ni + new_cand->nd;
      no = (int)floor(rndu()*nsum);
		SelectNeighbourPos(new_cand, seqlen, num_t, no);
      return;
   }

   if(flag==0){
      nsum=new_cand->ni + new_cand->nd;
      old_pos = (short *)calloc(nsum,sizeof(short));
   	memcpy(old_pos, new_cand->pos_vector, nsum*sizeof(short));
   }


   ok_flag=1;
   while(ok_flag){
      RandomSelectPosVector(new_cand, seqlen, num_t);
   	if(flag==0){
         pos_vector = new_cand->pos_vector;
         for(i=0;i<nsum;i++){
            if(old_pos[i]!=pos_vector[i]){ok_flag=0;break;}
         }
      }
      else{ok_flag=0;}
   }
   if(flag==0){old_pos = MemFree(old_pos);}

   return;
}

static void
SelectIndelArray(PoolStructPtr pool, CandStructPtr new_cand, short flag)
{
	int   i, max_num;
   short ni,nd, nind_index, nmax, lmax;
   double  *num_cand;
   int   select, distance, dist_limit;
   int   **i_pool, **d_pool;
   int   i_order, d_order;
   short  *ins_array, *del_array;
   short  *temp_ins, *temp_del;

   ni = new_cand->ni;
   nd = new_cand->nd;
   nmax = pool->nmax;
   lmax = pool->lmax;
   nind_index = (short)((nmax+1)*ni+nd);
   num_cand = pool->num_cand;
   i_pool = pool->i_pool;
   d_pool = pool->d_pool;
   ins_array = new_cand->ins_array;
   del_array = new_cand->del_array;

   if(flag<=0){dist_limit=30000;}/*large enough*/
   if(flag==1){dist_limit=10;}
   temp_ins = (short *)calloc(ni,sizeof(short));
   temp_del = (short *)calloc(nd,sizeof(short));

   distance=dist_limit;
   if(num_cand[nind_index]<0){printf("Warning!num_cand[%d] = 0\n",nind_index);getch();return;}
   else{max_num = round(pow(10.0,num_cand[nind_index]));}
   while((distance>=dist_limit)||((flag>=0)&&(distance==0))){
	   if(max_num==1){select=0;}
      else{select = floor(rndu()*max_num);}
   	i_order = i_pool[nind_index][select];
      d_order = d_pool[nind_index][select];
      /*printf("max_num=%d, select=%d, i_order=%d, d_order=%d\n",max_num,select,i_order,d_order);*/
      distance=0;
      for(i=0;i<ni;i++){
      	temp_ins[i]=(short)(i_order%(lmax+1));
         i_order /= (lmax+1);
         distance += (ins_array[i]-temp_ins[i])*(ins_array[i]-temp_ins[i]);
      }
      for(i=0;i<nd;i++){
      	temp_del[i]=(short)(d_order%(lmax+1));
         d_order /= (lmax+1);
         distance += (del_array[i]-temp_del[i])*(del_array[i]-temp_del[i]);
      }
   }
   for(i=0;i<ni;i++){ins_array[i]=temp_ins[i];}
   for(i=0;i<nd;i++){del_array[i]=temp_del[i];}

   temp_ins = MemFree(temp_ins);
   temp_del = MemFree(temp_del);


   return;
}

static short
SearchCandidate(int select, short lmax, int target_d, short ni, short nd, short *temp_ins, short *temp_del)
{
	int   i, d;
   int   max_iter, iter, num_passed;

   max_iter = 1;
   for(i=0;i<ni;i++){max_iter *= lmax;}
   for(i=0;i<nd;i++){max_iter *= lmax;}
   /*printf("max_iter = %d\n",max_iter); */

   /*initialization*/
   for(i=0;i<ni;i++){temp_ins[i]=1;}
   for(i=0;i<nd;i++){temp_del[i]=1;}
   iter=0;  num_passed=0;
   while(iter<max_iter){
      d = CheckLengthChange(temp_ins, ni, temp_del, nd);
      if(d==target_d){
      	num_passed++;
         if(num_passed==select){return 1;}
      }
      /*loop of del*/
      if(nd>0){
      	temp_del[nd-1]++;
      	CheckCarryUp(temp_del, nd, lmax);
      }
      if((temp_del[0]>lmax)||(nd==0)){
      	if(ni>0){
      		temp_ins[ni-1]++;
      		CheckCarryUp(temp_ins, ni, lmax);
   			for(i=0;i<nd;i++){temp_del[i]=1;}
         }
         if((temp_ins[0]>lmax)||(ni==0)){break;}
      }
      iter++;
   }

   if(num_passed==0){printf("No candidate found\n");/*getch();*/ }
   if(num_passed<select){return 0;}
   return 0;
}

static void
SearchIndelArray(PoolStructPtr pool, CandStructPtr new_cand, short flag)
{
	short i,ni,nd, nmax, lmax;
   int   nind_index, max_num1, max_num2;
   double *num_cand;
   int   select, select1, select2, distance, dist_limit;
   short *ins_array, *del_array;
   short *temp_ins, *temp_del;
   short success;

   ni = new_cand->ni;
   nd = new_cand->nd;
   nmax = pool->nmax;
   lmax = pool->lmax;
   nind_index = (nmax+1)*ni+nd;
   num_cand = pool->num_cand;
   ins_array = new_cand->ins_array;
   del_array = new_cand->del_array;
   temp_ins = (short *)calloc(nmax,sizeof(short));
   temp_del = (short *)calloc(nmax,sizeof(short));

   if(flag<=0){dist_limit=30000;}/*large enough*/
   if(flag==1){dist_limit=10;}
   distance=dist_limit;

   if(num_cand[nind_index]>4){
   	max_num1 = round(pow(10.0,num_cand[nind_index]-4.0)); max_num2 = round(pow(10.0,4.0));
   }
   else{
   	if(num_cand[nind_index]<0){printf("Warning!num_cand[%d]=0\n",nind_index);getch();}
      else{max_num1 = round(pow(10.0,num_cand[nind_index])); max_num2 = 1;}
   }
   while((distance>=dist_limit)||(distance==0)){
      if(max_num2==1){
      	select = floor(rndu()*max_num1);
      }
      else{
      	select1 = floor(rndu()*max_num1);
         select2 = floor(rndu()*max_num2);
         select = select1*10000+select2;
      }
      /*printf("select=%d\n",select); */
		success = SearchCandidate(select, lmax, pool->target_d, ni, nd, temp_ins, temp_del);
    if(success){
      distance=0;
      for(i=0;i<ni;i++){
         distance += (ins_array[i]-temp_ins[i])*(ins_array[i]-temp_ins[i]);
      }
      for(i=0;i<nd;i++){
         distance += (del_array[i]-temp_del[i])*(del_array[i]-temp_del[i]);
      }
    }
    else{distance=dist_limit;}
   }

   for(i=0;i<ni;i++){ins_array[i]=temp_ins[i];}
   for(i=0;i<nd;i++){del_array[i]=temp_del[i];}

   temp_ins = MemFree(temp_ins);
   temp_del = MemFree(temp_del);


   return;
}

short
ChangeIndelArray(PoolStructPtr pool, CandStructPtr new_cand, short flag)
{
	int nind_index;
   double  *num_cand;


   nind_index = (pool->nmax+1)*(new_cand->ni)+(new_cand->nd);
   num_cand = pool->num_cand;
   /*printf("num_cand[%d]=%d\n",nind_index, num_cand[nind_index]);  */
   if(num_cand[nind_index]<0){/*printf("cannot change indel array\n");*/return 1;}
   if(flag>=0){
   	if(num_cand[nind_index]==0.0){/*printf("cannot change indel array\n");*/return 1;}
   }

   if(num_cand[nind_index]<(log10((double)(pool->max_stored)))){
   	SelectIndelArray(pool, new_cand, flag);
   }
   else{
       /* if the number of candidates exceeds the 'max_stored', SearchIndelArray
is employed for ab initio choosing of cadidate. */
     
       SearchIndelArray(pool, new_cand, flag);
   }

   return 0;

}

void
MCSelectNIND(PoolStructPtr pool, int nind_index, CandStructPtr new_cand)
{
   short i, j, max_setnum;
   double **rij; /*relative possibility of set i->set j*/
   short ni, nd;
   double   *limit_rij, select;

   max_setnum = pool->max_setnum;
	rij = pool->rij;
   limit_rij = (double *)calloc(max_setnum+1,sizeof(double));
	limit_rij[0] = 0.0;
   for(i=0;i<max_setnum;i++){
   	limit_rij[i+1] = limit_rij[i]+rij[nind_index][i];
      /*printf("limit_rij[%d]=%d\n",(i+1), limit_rij[i+1]);getch(); */
   }

   select = rndu()*limit_rij[max_setnum];
   /*printf("select = %d\n",select); */
   j = -1;
   for(i=0;i<max_setnum;i++){
   	if((select>=limit_rij[i])&&(select<limit_rij[i+1])){j=i;break;}
   }
   if(j==-1){printf("Warning!select j is failed(press any key)\n");getch();}
   /*new nind_index=j*/
   ni = (short)(j / (pool->nmax+1));
   nd = (short)(j % (pool->nmax+1));
   /*printf("ni=%d, nd=%d\n",ni,nd); getch(); */
   new_cand->ni = ni;
   new_cand->nd = nd;

   return;
   
}



short
SelectNewCandidate(PoolStructPtr pool, PoolConsPtr pool_cons, CandStructPtr new_cand, CandStructPtr current_cand, int current_id)
{
   int   i, nind_index, num_chain;
   short ni, nd, chain_no;
   double  *num_cand;
   double *prob_chain, *limit_chain;
   double chance;
	short failed;
   int   loop;

   if((new_cand==NULL)||(current_cand==NULL)){
   	printf("candidate is NULL\n");return 1;
   }
   if(current_id==0){
      /*printf("first candidate\n"); */
      /*fprintf(fout,"first candidate\n");*/
      new_cand->chain_no = 5;
      nind_index = 0;
   	failed = 1;
   	while(failed){
			MCSelectNIND(pool, nind_index, new_cand);
      	failed = ChangeIndelArray(pool, new_cand, -1);
         /*printf("change indel array done\n");*/
			if(!failed){
         	failed=SelectOrder(new_cand->order, new_cand->ni, new_cand->nd, -1);
            /*printf("select order done\n"); */
         }
			if(!failed){
         	failed=SelectTvector(new_cand->t_vector, pool->num_t, (short)(new_cand->ni+new_cand->nd), -1);
         	/*printf("select t vector done\n"); */
         }
			if(!failed){
         	SetupIDarray(new_cand, pool->num_t);
            /*printf("setup id array done\n"); */
         }
			if(!failed){
				SelectPosVector(new_cand, pool->orig_seqlen, pool->num_t, -1);
      		/*printf("select pos done\n");*/
         }
   		if(!failed){
         	failed=CountNumPosCand(new_cand, pool->orig_seqlen, pool->num_t);
         	if(failed)printf("seqlen<=0\n");
      	}
   	}
      return 0;
   }
   else{
      /*printf("new candidate is proposed\n"); */
      /*fprintf(fout,"new candidate is proposed\n");*/
      ni = current_cand->ni;
      nd = current_cand->nd;
      nind_index = (pool->nmax+1)*ni+nd;
      if(nind_index==0){
         chance = rndu();
         if(chance<0.1){/*re-propose (0,0) */
         	CopyCandidate(new_cand, current_cand, pool->num_t);
         	return 0;
			}
         else{chain_no=5;}
      }
      else{
      	num_chain = 8;
      	prob_chain = (double *)calloc(num_chain,sizeof(double));
      	/*
      	chain_no 0:change t vector
      	chain_no 1:change pos vector
      	chain_no 2:change order vector
      	chain_no 3:change ins,del array
         chain_no 4:conserve first component, change (ni,nd)
      	chain_no 5:change ni,nd
      	chain_no 6:select neighbour t vector
      	chain_no 7:select neughbour pos vector
     In the following, relative frequencies for different ways to propose
('jump' to) a new candidate are specified

*/
      	prob_chain[0]=0.5;
      	prob_chain[1]=4.0;
      	if(ni*nd==0) {prob_chain[2]=0.0;}
      	else{prob_chain[2]=0.1;}
      	/*check num_cand of indel_array*/
   		num_cand = pool->num_cand;
      	if(num_cand[nind_index]<=0.0){prob_chain[3]=0.0;}
      	else {prob_chain[3]=1.0;}
			prob_chain[4]=1.0;
      	prob_chain[5]=10.0;
      	prob_chain[6]=1.0;
      	prob_chain[7]=4.0;
      	limit_chain = (double *)calloc(num_chain+1, sizeof(double));
      	limit_chain[0]=0;
      	for(i=0;i<num_chain;i++){
      		limit_chain[i+1] = limit_chain[i]+prob_chain[i];
      	}
       	/*select chain*/
	   	chance = rndu()*limit_chain[num_chain];
      	for(i=0;i<num_chain;i++){
      		if((chance>=limit_chain[i])&&(chance<limit_chain[i+1])){chain_no=(short)i;break;}
      	}
         free(prob_chain); free(limit_chain);
      }

      CopyCandidate(new_cand, current_cand, pool->num_t);

      new_cand->chain_no=chain_no;
      switch(chain_no){
      	case 0: /*printf("case 0\n"); */
                 new_cand->chain_no=0;
                 failed = SelectTvector(new_cand->t_vector, pool->num_t, (short)(ni+nd), 0);
                 if(failed){printf("Warning!cannot change t_vector\n");/*getch();*/}
                 else{
                  SetupIDarray(new_cand, pool->num_t);
   					failed=CountNumPosCand(new_cand, pool->orig_seqlen, pool->num_t);
                  if(failed){printf("Warning!seqlen<=0\n");getch();return 1;}
                	break;
                 }
      	case 1: /*printf("case 1\n"); */
                 new_cand->chain_no=1;
                 SelectPosVector(new_cand, pool->orig_seqlen, pool->num_t, 0);
                 break;
      	case 2: /*printf("case 2\n");*/
                 new_cand->chain_no=2;
                 failed = SelectOrder(new_cand->order, ni, nd, 0);
                 if(!failed){
                 	SetupIDarray(new_cand, pool->num_t);
                  SelectPosVector(new_cand, pool->orig_seqlen, pool->num_t, 2);
                  failed=CountNumPosCand(new_cand, pool->orig_seqlen, pool->num_t);
         			if(failed)printf("seqlen<=0\n");
                 	if(!failed){break;}
                 }
      	case 3: /*printf("case 3\n"); */
                 new_cand->chain_no=3;
                 failed = ChangeIndelArray(pool, new_cand, 0);
                 if(failed){printf("Warning!no candidate\n");/*getch();*/}
                 else{
                 	SetupIDarray(new_cand, pool->num_t);
                  SelectPosVector(new_cand, pool->orig_seqlen, pool->num_t, 2);
   					failed=CountNumPosCand(new_cand, pool->orig_seqlen, pool->num_t);
         			if(failed)printf("seqlen<=0\n");
                 	if(!failed)break;
                 }
         case 4: /*printf("case 4\n"); */
                 new_cand->chain_no=4;
                 failed=SelectConserveChainTPos(pool_cons, pool, new_cand, current_cand);
                 if(!failed) {
                 	failed=CountNumPosCand(new_cand, pool->orig_seqlen, pool->num_t);
         			if(failed)printf("seqlen<=0\n");
						if(!failed)break;
                 }
                 /*if(failed) proceed to case 5*/
                 CopyCandidate(new_cand, current_cand, pool->num_t);
                 new_cand->chain_no=5;
         case 5: /*printf("case 5\n"); */
                 failed = 1;
                 loop=0;
         		  while(failed){
                  loop++;
                  if(loop>10){printf("loop=%d\n",loop); /*getch();*/}
                 	MCSelectNIND(pool, nind_index, new_cand);
                  /*printf("(ni,nd)=(%d,%d)\n",new_cand->ni,new_cand->nd); */
                  if((new_cand->ni==0)&&(new_cand->nd==0)){
                     /*printf("nind_index=0\n"); */
                     if(pool->target_d==0){SetCandToZero(new_cand, pool->nmax, pool->num_t); failed=0;return 0;}
                     else{failed=1;}
                  }
                  else{
      					failed = ChangeIndelArray(pool, new_cand, -1);
							if(!failed)failed=SelectOrder(new_cand->order, new_cand->ni, new_cand->nd, -1);
							if(!failed)failed=SelectTvector(new_cand->t_vector, pool->num_t, (short)(new_cand->ni+new_cand->nd), 0);
							if(!failed)SetupIDarray(new_cand, pool->num_t);
							if(!failed)SelectPosVector(new_cand, pool->orig_seqlen, pool->num_t, -1);
                 		if(!failed){
                     	failed=CountNumPosCand(new_cand, pool->orig_seqlen, pool->num_t);
         					if(failed)printf("seqlen<=0\n");
                     }
                  }
                 }
                 break;
      	case 6: /*printf("case 6\n");*/
                 new_cand->chain_no=6;
                 failed = SelectTvector(new_cand->t_vector, pool->num_t, (short)(ni+nd), 1);
                 if(failed){printf("Warning!cannot change t_vector\n");/*getch();*/}
                 else{
                  SetupIDarray(new_cand, pool->num_t);
   					failed=CountNumPosCand(new_cand, pool->orig_seqlen, pool->num_t);
         			if(failed){printf("Warning!seqlen<=0\n");return 1;}
                	break;
                 }
      	case 7: /*printf("case 7\n"); */
                 new_cand->chain_no=7;
                 SelectPosVector(new_cand, pool->orig_seqlen, pool->num_t, 1);
                 break;
         default:return 1;
      }
   }

   return 0;

}

static short
CountNumPosCand(CandStructPtr new_cand, int seqlen, int num_t)
{
   int    t, current_seqlen;
   short  *id_array;

   id_array = new_cand->id_array;
   new_cand->num_pos_cand = 0.0;
   current_seqlen = seqlen;
   if(current_seqlen<=0){new_cand->num_pos_cand=-1;return 1;}
   for(t=0;t<num_t;t++){
   	if(id_array[t]>0){
      	new_cand->num_pos_cand += log10(current_seqlen+1);
         current_seqlen += id_array[t];
      }
      if(id_array[t]<0){
         if((current_seqlen+id_array[t])<=0){new_cand->num_pos_cand=-1;return 1;}
      	new_cand->num_pos_cand += log10(current_seqlen+id_array[t]+1);
         current_seqlen += id_array[t];
      }
   }
   if(new_cand->num_pos_cand<0){printf("Warning!num_pos_cand<0\n");return 1;}
   else{return 0;}

}


/* This ensures that new candidate is different from the old one. */
static short
CheckNewCandidate(CandStructPtr new_cand, CandStructPtr current_cand, int num_t)
{
   int    i, nsum;
   short  *id_array1, *id_array2;
   short  *pos_vector1, *pos_vector2;

   if(new_cand->ni!=current_cand->ni){return 0;}
   if(new_cand->nd!=current_cand->nd){return 0;}

   id_array1 = new_cand->id_array;
   id_array2 = current_cand->id_array;
   pos_vector1 = new_cand->pos_vector;
   pos_vector2 = new_cand->pos_vector;

   for(i=0;i<num_t;i++){
   	if(id_array1[i]!=id_array2[i]){return 0;}
   }
   nsum = new_cand->ni + new_cand->nd;
   for(i=0;i<nsum;i++){
   	if(pos_vector1[i]!=pos_vector2[i]){return 0;}
   }
   return 1;

}


/* At present, this function is not used .  */
short
ModSelectNewCandidate(PoolStructPtr pool, CandStructPtr new_cand, CandStructPtr current_cand, int current_id)
{
	int i, max_setnum, nind_index;
   double *limit_nind, r;
   short  failed;


   if((new_cand==NULL)||(current_cand==NULL)){
   	printf("candidate is NULL\n");return 1;
   }

   if(current_id>0){CopyCandidate(new_cand, current_cand, pool->num_t);}
   /*select (ni,nd)*/
   max_setnum = pool->max_setnum;
   limit_nind = pool->limit_nind;
   if(limit_nind==NULL){printf("ModSelectNewCandidate failed\n");getch();return 1;}

   nind_index=-1;
   while(nind_index<0){
   	r = rndu()*limit_nind[max_setnum];
   	for(i=0;i<max_setnum;i++){
   		if((r>=limit_nind[i])&&(r<limit_nind[i+1])){nind_index=i;break;}
   	}
   }
   new_cand->ni = (short)(nind_index / (pool->nmax+1));
   new_cand->nd = (short)(nind_index % (pool->nmax+1));
   /*printf("(ni,nd)=(%d,%d)\n",new_cand->ni,new_cand->nd);*/

   /*select indel*/
   failed = ChangeIndelArray(pool, new_cand, -1);
   if(failed){return 1;}
   failed=SelectOrder(new_cand->order, new_cand->ni, new_cand->nd, -1);
   if(failed){return 1;}
   failed=SelectTvector(new_cand->t_vector, pool->num_t, (short)(new_cand->ni+new_cand->nd), -1);
   if(failed){return 1;}
   SetupIDarray(new_cand, pool->num_t);
   SelectPosVector(new_cand, pool->orig_seqlen, pool->num_t, -1);
   CountNumPosCand(new_cand, pool->orig_seqlen, pool->num_t);
   if(current_id>0){
   	failed = CheckNewCandidate(new_cand, current_cand, pool->num_t);
      if(failed){return 1;}
   }
	return 0;

}

int
CheckMinLength(PoolStructPtr pool, CandStructPtr new_cand)
{
   int    min_len, t;
   short  *id_array;

   min_len = pool->orig_seqlen;
   id_array = new_cand->id_array;
   for(t=0;t<(pool->num_t);t++){
   	if(id_array[t]!=0){min_len += id_array[t];}
      if(min_len<=0){return min_len;}
   }
   return min_len;

}


short
ProcessMCMC(StockIDPtr stock_id, OptionsBlkPtr options, SubsMatrixPtr subs_matrix, FILE *fout)
{
   long int      iter, visited, current_state;
	CandStructPtr  current_cand;
   CandStructPtr  new_cand;
   PoolStructPtr   pool;
	PoolConsPtr	    pool_cons;
   int             seqlen, target_d, min_len;
   short           failed, accepted, end_flag;
   ChainIDPtr PNTR chain_id_array;
   ChainIDPtr      chain_id;

   current_cand = CandStructAllocate(options->nmax, options->num_t);
   if(current_cand==NULL){printf("current candidate allocation failure\n");return 1;}
   new_cand = CandStructAllocate(options->nmax, options->num_t);
   if(new_cand==NULL){printf("new candidate allocation failure\n");return 1;}

	chain_id_array = stock_id->chain_id_array;
   seqlen = options->orig_len;
   target_d = (options->target_len) - (options->orig_len);

   pool = SetupPoolStruct(options->lmax, options->nmax, target_d, seqlen, options->num_t);
   if(pool==NULL){printf("pool allocation failure\n");return 1;}
   /*DisplayPoolStruct(pool, fout);*/

   /*the struct for conserved (ni,nd) changing*/
   pool_cons = PoolConsAllocate(pool->max_setnum, 10000);
   if(pool_cons==NULL){return 1;}

   /*stock_id->current_id = 0; */
   current_state = -1;
   for(iter=0;iter<(options->max_iter);iter++){
      if(!(iter%500)){printf("mcmc %d\n",iter);}
   	failed = SelectNewCandidate(pool, pool_cons, new_cand, current_cand, stock_id->current_id);
      if(failed){
   		new_cand = CandDestruct(new_cand);
   		current_cand = CandDestruct(current_cand);
         return 1;
      }
      min_len = CheckMinLength(pool, new_cand);
      while(min_len<=0){
         printf("min_len=%d\n",min_len);
   		failed = SelectNewCandidate(pool, pool_cons, new_cand, current_cand, stock_id->current_id);
      	if(failed){
   			new_cand = CandDestruct(new_cand);
   			current_cand = CandDestruct(current_cand);
         	return 1;
      	}
      	min_len = CheckMinLength(pool, new_cand);
      }

      visited = -2;
      if(stock_id->current_id>0){
      	visited = CheckRevisit(new_cand, stock_id, options);
      }
      if(visited==-3){printf("CheckRevisit failed\n");getch(); return 1;}
      /*calculate likelihood in either visited case or not*/
      TrialOfNewCandidate(new_cand, options, subs_matrix, fout);

      if(visited==-2){accepted=1;}/*first candidate*/
      else{
         /*flag=0:SelectNewCandidate, flag=1:ModSelectNewCandidate*/
			accepted = CompareCandidate(new_cand, current_cand, pool, pool_cons, 0, fout);
      }
		if(accepted) {
      	StockCandidate(stock_id, new_cand, current_cand, pool, visited);
         if(visited<0){current_state = stock_id->current_id - 1;}
         else{current_state = visited;}
      }
   	else {
         if(current_state>=0){
         	chain_id = chain_id_array[current_state];
            chain_id->num_visit+=1;
         }
      }
      /*throw away first "num_throw" stock_id */
      if(iter==(options->num_throw)){
			InitializeStockID(stock_id, (short)(2*(options->nmax)));
      	ArrayToID(current_cand, chain_id_array[0], pool->lmax, pool->nmax, pool->num_t);
      	stock_id->current_id=1;
         current_state = 0;
      }
   }

   StatisticalResult(stock_id, pool->max_setnum, pool->nmax, fout);



   new_cand = CandDestruct(new_cand);
   current_cand = CandDestruct(current_cand);
   pool = PoolDestruct(pool);
   pool_cons = PoolConsDestruct(pool_cons);

   return 0;
}



/* This has been newly added to introduce the thermalization step, which 
occasionally accepts a candidate regardless of its likelihood*prior-probability */
short
ExModProcessMCMC(StockIDPtr stock_id, OptionsBlkPtr options, SubsMatrixPtr subs_matrix, FILE *fout)
{
   long int      iter, visited, current_state, max_iter;
	CandStructPtr  current_cand;
   CandStructPtr  new_cand;
   PoolStructPtr   pool;
	PoolConsPtr	    pool_cons;
   int             seqlen, target_d, min_len, num_throw;
   short           failed, accepted/*, end_flag*/;
   ChainIDPtr PNTR chain_id_array;
   ChainIDPtr      chain_id;
	int          num;
	StockIDPtr   each_id;

   current_cand = CandStructAllocate(options->nmax, options->num_t);
   if(current_cand==NULL){printf("current candidate allocation failure\n");return 1;}
   new_cand = CandStructAllocate(options->nmax, options->num_t);
   if(new_cand==NULL){printf("new candidate allocation failure\n");return 1;}

   seqlen = options->orig_len;
   target_d = (options->target_len) - (options->orig_len);

   pool = SetupPoolStruct(options->lmax, options->nmax, target_d, seqlen, options->num_t);
   if(pool==NULL){printf("pool allocation failure\n");return 1;}
   /*DisplayPoolStruct(pool, fout);*/

   /*the struct for conserved (ni,nd) changing*/
   pool_cons = PoolConsAllocate(pool->max_setnum, 10000);
   if(pool_cons==NULL){return 1;}

   each_id = StockIDAllocate(options->max_idnum, options->nmax, options->select_num);
   if(each_id==NULL){
   	printf("Each ID Allocate failed(press any key)\n");getch();
      return 1;
   }
	chain_id_array = each_id->chain_id_array;

for(num=0;num<10;num++){
	InitializeStockID(each_id, (short)(2*(options->nmax)));
   max_iter = options->max_iter;
   num_throw = options->num_throw;
   current_state = -1;
   for(iter=0;iter<max_iter;iter++){
      if(!(iter%500)){printf("mcmc %d\n",iter);}
   	failed = SelectNewCandidate(pool, pool_cons, new_cand, current_cand, each_id->current_id);
      if(failed){
   		new_cand = CandDestruct(new_cand);
   		current_cand = CandDestruct(current_cand);
         return 1;
      }
      min_len = CheckMinLength(pool, new_cand);
      while(min_len<=0){
         printf("min_len=%d\n",min_len);
   		failed = SelectNewCandidate(pool, pool_cons, new_cand, current_cand, each_id->current_id);
      	if(failed){
   			new_cand = CandDestruct(new_cand);
   			current_cand = CandDestruct(current_cand);
         	return 1;
      	}
      	min_len = CheckMinLength(pool, new_cand);
      }

      visited = -2;
      if(each_id->current_id>0){
      	visited = CheckRevisit(new_cand, each_id, options);
      }
      if(visited==-3){printf("CheckRevisit failed\n");getch(); return 1;}

      /*calculate likelihood in either visited case or not*/
		TrialOfNewCandidate(new_cand, options, subs_matrix, fout);

      if(visited==-2){accepted=1;}/*first candidate*/
      else{
         /*flag=0:SelectNewCandidate, flag=1:ModSelectNewCandidate*/
			accepted = CompareCandidate(new_cand, current_cand, pool, pool_cons, 0, fout);
      }
		if(accepted) {
      	StockCandidate(each_id, new_cand, current_cand, pool, visited);
         if(visited<0){current_state = each_id->current_id - 1;}
         else{current_state = visited;}
      }
   	else {
         if(current_state>=0){
         	chain_id = chain_id_array[current_state];
            chain_id->num_visit+=1;
         }
      }
      /*check and throw away the each_id of the burn-in period*/
      if(iter==num_throw){
         InitializeStockID(each_id, (short)(2*(options->nmax)));
         ArrayToID(current_cand, chain_id_array[0], pool->lmax, pool->nmax, pool->num_t);
         each_id->current_id=1;
         current_state = 0;
         if(current_cand->lnL<=Penalty){
         	num_throw += 200; max_iter += 200;
         }
      }
      /*if(iter==(num_throw+99)){
      	end_flag=CheckMonopolyOfStockID(each_id, fout);
         if(end_flag){printf("Monopoly\n");break;}
      } */
   }
   SortStockID(each_id, options->nmax);

   AddEachIDtoStockID(stock_id, each_id, options->nmax);
}
StatisticalResult(stock_id, pool->max_setnum, pool->nmax, fout);



   new_cand = CandDestruct(new_cand);
   current_cand = CandDestruct(current_cand);
   pool = PoolDestruct(pool);
   pool_cons = PoolConsDestruct(pool_cons);
   each_id = StockIDFree(each_id);

   return 0;
}