#include<stdlib.h>
#include<stdio.h>
#include<math.h>

#define PI 3.14159265

//do product recursively
//fixed x, n, g, p, but different mu, sigma for different p; need an index vector id[] (i_j in equation 53) for this
//remember both k and j start from 1 to n and p, respectively; need to subtract one for index consistency
double prod_of(long *n, long *g, long *p, long *k, long *j, long *id, double *x, double *mu, double *sigma){
    int i;
	double prod;
	long jnext;

	if ((*j)==1){
         prod=exp(0-0.5*(x[((*k)-1)*(*p)+(*j)-1]-mu[((*j)-1)*(*g)+id[(*j)-1]])*(x[((*k)-1)*(*p)+(*j)-1]-mu[((*j)-1)*(*g)+id[(*j)-1]])/sigma[((*j)-1)*(*g)+id[(*j)-1]])/sqrt(2*PI*sigma[((*j)-1)*(*g)+id[(*j)-1]]);      
	}
	else {
         jnext = (*j)-1;
         prod=prod_of(n,g,p,k,&jnext,id,x,mu,sigma)*exp(0-0.5*(x[((*k)-1)*(*p)+(*j)-1]-mu[((*j)-1)*(*g)+id[(*j)-1]])*(x[((*k)-1)*(*p)+(*j)-1]-mu[((*j)-1)*(*g)+id[(*j)-1]])/sigma[((*j)-1)*(*g)+id[(*j)-1]])/sqrt(2*PI*sigma[((*j)-1)*(*g)+id[(*j)-1]]);       
	}

	return(prod);
	/*printf("%lf \n", prod);*/
}

//make an index vector for each of combination of dimensional components
//remember j start from 1 to p; need to subtract one for index consistency
long index_of(long *g, long *j, long *id){
	long idx;
	long jnext;
	
	if ((*j)==1){
		idx=id[(*j)-1];
	}
	else {
         jnext=(*j)-1;
		 idx=index_of(g,&jnext,id)*(*g)+id[(*j)-1];
	}

	return(idx);
}

//calculate partial density for diagonal elements
double fpcdpart1_k(long *n, long *g, long *p, long *k, double *pi, double *alpha, double *x, double *mu, double *sigma){
   double total;
   int l,jj;
   long id[(*p)];
   total=0;

   for (l=0; l<(*g); l++){
        for (jj=0; jj<(*p); jj++){
   	         id[jj]=l;   	    
        }
   
   	    total=total+(pi[l]-(*alpha))*prod_of(n, g, p, k, p, id, x, mu, sigma);  	    
   }
   return(total);   	
}

//calculate partial density for non-diagonal elements
//do multiple sums using for and also recursion
//remember j start from 1 to p; need to subtract one for index consistency
double nest_of_fpcdpart2_k(long *n, long *g, long *p, long *k, long *j, long *id, double *x, double *mu, double *sigma){
       double total;
       int l;
       long jnext;

       l=0-1;       
       total=0;
      
       if ((*j)<=(*p)){
                    for(l=0;l<(*g);l++){
                                     id[(*j)-1]=l;
                                     jnext=(*j)+1;
                                     total=total+nest_of_fpcdpart2_k(n, g, p, k, &jnext, id, x, mu, sigma);                                     
                    }
       }
       else {
            total=prod_of(n, g, p, k, p, id, x, mu, sigma);
       }
       
 /*printf("%d %d %d %lf \n", *k, *j, l, total);*/

       return(total);
}



//calculate uhat which includes alpha
void nest_ez_k(long *n, long *g, long *p, long *k, long *j, long *id, double *uhat, double *pi, double *alpha, double *x, double *mu, double *sigma){
	double temp1, temp2, temp4, sumpi;
	int l;
    long jnext, idx, ii,jj,temp, temp3, id1[(*p)];
  	fflush(stdin);
    l=0-1;  
    jj=1;
    temp=pow((*g),(*p));
	if ((*j)<=(*p)) {
		for (l=0; l<(*g); l++){
		     id[(*j)-1]=l; 
		     jnext=(*j)+1;
			 nest_ez_k(n,g,p,k,&jnext, id, uhat, pi, alpha, x, mu, sigma);
		}
 	    /*for(ii=0; ii<(*p); ii++){
            printf("%d ", id[ii]);
    	}*/		
	}
	
 
	else { 
        temp1=fpcdpart1_k(n, g, p, k, pi, alpha, x, mu, sigma);
        temp2=nest_of_fpcdpart2_k(n, g, p, k, &jj, id1, x, mu, sigma);
        temp4=prod_of(n, g, p, k, p, id, x, mu, sigma);
        /*printf("%lf \n",temp1);*/
        temp3=index_of(g,p,id);
        /*printf("%d", temp3);
        printf("\n");*/
        uhat[((*k)-1)*temp+temp3]=((*alpha)*temp4)/(temp1+(*alpha)*temp2);
    	/*for(ii=0; ii<(*p); ii++){
            printf("%d ", id[ii]);
    	}
    	printf("\n");*/
	}
	
}


double nest_of_mu_jj_num_k(long *n, long *g, long *p, long *k, long *j, long *jj, long *ijj, long *id ,double *uhat, double *pi, double *x, double *mu, double *sigma){
       double total;
       int l;
       long jnext, temp1, temp;
       fflush(stdin);
       temp=pow((*g),(*p));

       l=0-1;       
       total=0;
       
       if ((*j)<=(*p)){
       	        if ((*j)==(*jj)){
                    for(l=0;l<(*g);l++){
                    	if (l==(*ijj)){
               	            id[(*j)-1]=l;
                            jnext=(*j)+1;
                            total=total+nest_of_mu_jj_num_k(n, g, p, k, &jnext, jj, ijj, id, uhat, pi, x, mu, sigma);   	
                    	/*printf("%d %d %d %d %d %lf \n", *k, *j, l,  *jj, *ijj, total);*/
						}  
                    	else{
                            total=total+0;   
                        /*printf("%d %d %d %d %d %lf \n", *k, *j, l, *jj, *ijj, total);*/
	                    	
	                    }
                        
                    }
                }
                else{
                	for(l=0;l<(*g);l++){
         	            id[(*j)-1]=l;
                        jnext=(*j)+1;
                        total=total+nest_of_mu_jj_num_k(n, g, p, k, &jnext, jj, ijj, id, uhat, pi, x, mu, sigma);   	 
                	    /*printf("%d %d %d %d %d %lf \n", *k, *j, l, *jj, *ijj, total);*/
					}
				}
       }
       else {
 	        if ((*j)!=(*jj)){   
		        temp1=index_of(g,p,id);
                total=uhat[((*k)-1)*temp+temp1]*x[((*k)-1)*(*p)+(*jj)-1];
 	        }
            /*else {
                  if (id[(*p)-1]==(*ijj)) {
                      temp1=index_of(g,p,id);
                      total=uhat[((*k)-1)*temp+temp1]*x[((*k)-1)*(*p)+(*jj)-1];    	              	
                  }      
            }*/
       }      
    /*printf("%d %d %d %lf \n", *k, *j, l, total);*/
       return(total);
}




double nest_of_mu_jj_den_k(long *n, long *g, long *p, long *k, long *j, long *jj, long *ijj, long *id, double *uhat, double *pi, double *x, double *mu, double *sigma){
       double total;
       int l;
       long jnext, temp1, temp;
       fflush(stdin);
       temp=pow((*g),(*p));

       l=0-1;       
       total=0;
       
       if ((*j)<=(*p)){
       	        if ((*j)==(*jj)){
                    for(l=0;l<(*g);l++){
                    	if (l==(*ijj)){
               	            id[(*j)-1]=l;
                            jnext=(*j)+1;
                            total=total+nest_of_mu_jj_den_k(n, g, p, k, &jnext, jj, ijj, id, uhat, pi, x, mu, sigma);   	
                    	/*printf("%d %d %d %d %d %lf \n", *k, *j, l,  *jj, *ijj, total);*/
						}  
                    	else{
                            total=total+0;   
                        /*printf("%d %d %d %d %d %lf \n", *k, *j, l, *jj, *ijj, total);*/
	                    	
	                    }
                        
                    }
                }
                else{
                	for(l=0;l<(*g);l++){
         	            id[(*j)-1]=l;
                        jnext=(*j)+1;
                        total=total+nest_of_mu_jj_den_k(n, g, p, k, &jnext, jj, ijj, id, uhat, pi, x, mu, sigma);   	 
                	    /*printf("%d %d %d %d %d %lf \n", *k, *j, l, *jj, *ijj, total);*/
					}
				}
       }
       else {
 	         if ((*j)!=(*jj)){   
		        temp1=index_of(g,p,id);
                total=uhat[((*k)-1)*temp+temp1];
 	        }
            /*else {
                  if (id[(*p)-1]==(*ijj)) {
                      temp1=index_of(g,p,id);
                      total=uhat[((*k)-1)*temp+temp1];   	              	
                  }      
            }*/
       }      
       /*printf("%d %d %d %lf \n", *k, *j, l, total);*/
       return(total);
}



double nest_of_sigma_jj_num_k(long *n, long *g, long *p, long *k, long *j, long *jj, long *ijj, long *id, double *uhat, double *pi, double *x, double *mu, double *sigma){
       double total;
       int l;
       long jnext, temp1, temp;
       fflush(stdin);
       temp=pow((*g),(*p));

       l=0-1;       
       total=0;
       
       if ((*j)<=(*p)){
       	        if ((*j)==(*jj)){
                    for(l=0;l<(*g);l++){
                    	if (l==(*ijj)){
               	            id[(*j)-1]=l;
                            jnext=(*j)+1;
                            total=total+nest_of_sigma_jj_num_k(n, g, p, k, &jnext, jj, ijj,id, uhat, pi, x, mu, sigma);   	
                    	/*printf("%d %d %d %d %d %lf \n", *k, *j, l,  *jj, *ijj, total);*/
						}  
                    	else{
                            total=total+0;   
                        /*printf("%d %d %d %d %d %lf \n", *k, *j, l, *jj, *ijj, total);*/
	                    	
	                    }
                        
                    }
                }
                else{
                	for(l=0;l<(*g);l++){
         	            id[(*j)-1]=l;
                        jnext=(*j)+1;
                        total=total+nest_of_sigma_jj_num_k(n, g, p, k, &jnext, jj, ijj,id, uhat, pi, x, mu, sigma);   	 
                	    /*printf("%d %d %d %d %d %lf \n", *k, *j, l, *jj, *ijj, total);*/
					}
				}
       }
       else {
             if ((*j)!=(*jj)){   
		        temp1=index_of(g,p,id);
                total=uhat[((*k)-1)*temp+temp1]*(x[((*k)-1)*(*p)+(*jj)-1]-mu[((*jj)-1)*(*g)+(*ijj)])*(x[((*k)-1)*(*p)+(*jj)-1]-mu[((*jj)-1)*(*g)+(*ijj)]);
             }
             /*else {
                  if (id[(*p)-1]==(*ijj)) {
                      temp1=index_of(g,p,id);
                      total=uhat[((*k)-1)*temp+temp1]*(x[((*k)-1)*(*p)+(*jj)-1]-mu[((*jj)-1)*(*g)+(*ijj)])*(x[((*k)-1)*(*p)+(*jj)-1]-mu[((*jj)-1)*(*g)+(*ijj)]);
 	              }      
             }*/
 	       
       }      
    /*printf("%d %d %d %lf \n", *k, *j, l, total);*/
       return(total);
}

void em_pdim_pcd_exchange(long *n, long *g, long *p, long *j, long *id,  double *uhat, double *pi, double *alpha, double *x, double *mu, double *sigma, double *loglik,double *tol, long *restriction, long *constrain, long *iteration, long *convergence){
     long h,l,jj,ij,jjj,ijj,idx,*id1,*id2,power;
     int i, ii,flag, iter;
     double *Lg,*logLg,loglik2,sumpi,temppi_den, tempmu_num_sum, tempmu_den_sum, tempsigma_num_sum, *tempmu_num_jj,*tempmu_den_jj,*tempsigma_num_jj, *fpcd;
     power=pow((*g),(*p));
	 id1 = (long *) malloc((*p) * sizeof(long));
     id2 = (long *) malloc((*p) * sizeof(long));
     tempmu_num_jj= (double *) malloc((*n) * sizeof(double));
     tempmu_den_jj = (double *) malloc((*n)* sizeof(double));
     tempsigma_num_jj = (double *) malloc((*n) * sizeof(double));
     fpcd = (double *) malloc((*n) * sizeof(double));
     logLg = (double *) malloc((*n) * sizeof(double));


/*theoretical restriction on the null*/
if((*restriction)>0){    
    /*initialization*/
    flag=1;
    (*loglik) = 0;
    	for(h=1;h<=(*n);h++){
          	    for(i=0;i<(*p);i++){
                    id[i]=0;
                }
	            fpcd[h-1]=fpcdpart1_k(n, g, p, &h, pi, alpha, x, mu, sigma)+nest_of_fpcdpart2_k(n, g, p, &h, j, id, x, mu, sigma);
	            logLg[h-1] = log(fpcd[h-1]); 
	            (*loglik)= (*loglik) + logLg[h-1];
//                printf("%d %lf %lf \n", h, fpcd[h-1], logLg[h-1]);
   	    }
           
        
        printf("%lf \n", *loglik);
    /*iteration*/
    iter = 0;
    while(flag>0 & iter<(*iteration)){
    	printf("%d \n",iter);
        /*E step for z*/
      	    	  for(h=1;h<=(*n);h++){
                      for(i=0;i<(*p);i++){
                          id[i]=0;
                      }
                      nest_ez_k(n, g, p, &h, j, id, uhat, pi,alpha, x, mu, sigma);
             	  }
                  for (l=0;l<power*(*n);l++){
//	                   printf("%d %lf", l, uhat[l]);
                  }
                  for(h=1;h<=(*n);h++){
                      for (l=0;l<(*g);l++){
                  	       for (jj=0;jj<(*p);jj++){
   	                  	       id[jj]=l;
   	                       } 
   	                       uhat[(h-1)*power+index_of(g,p,id)]=pi[l]*uhat[(h-1)*power+index_of(g,p,id)]/(*alpha);
							  //prod_of(n, g, p, &h, p, id, x, mu, sigma)/fpcd[h-1];
                       }
                  }
	
        /*M step for pi*/
                 temppi_den=0;
    	          for(l=0;l<power*(*n);l++){      	    	
   	                  temppi_den=temppi_den+uhat[l];
          	      }

                  for(l=0;l<(*g);l++){
  	                  pi[l]=0;
  	                  for(jj=0;jj<(*p);jj++){
                  	      id[jj]=l;
                      }
	                  for(h=0;h<(*n);h++){
                      	pi[l]=pi[l]+uhat[h*power+index_of(g,p,id)];
	        
	                  }
                      pi[l]=pi[l]/temppi_den;
                      printf("%d %lf \n", l ,pi[l]);
        	      }  

        	      
        /*M step for alpha*/
                  sumpi=0;
                  for(l=0;l<(*g);l++){
                  	 sumpi=sumpi+pi[l];
                  }
                  printf("%lf \n", sumpi);
                  
		          (*alpha)=(1.0-sumpi)/(power-(*g));
		          printf("%lf \n", *alpha);
     
	    /*M step for mu*/		

                  for(jj=1;jj<=(*p);jj++){
                      for(ij=0;ij<(*g);ij++){
                          fflush(stdin);
		                  tempmu_num_sum=0;
               	          tempmu_den_sum=0;	
	                      for(h=1;h<=(*n);h++){
		                      for(ii=0;ii<(*p);ii++){
	                              id1[ii]=0;
	                              id2[ii]=0;
                              }      
                              tempmu_num_jj[h-1]=nest_of_mu_jj_num_k(n, g, p, &h, j, &jj, &ij, id1, uhat, pi, x, mu, sigma);
                              tempmu_den_jj[h-1]=nest_of_mu_jj_den_k(n, g, p, &h, j, &jj, &ij, id2, uhat, pi, x, mu, sigma);
                              tempmu_num_sum=tempmu_num_sum+tempmu_num_jj[h-1];
                              tempmu_den_sum=tempmu_den_sum+tempmu_den_jj[h-1];
//                              printf("%d %d %d %lf %lf\n", jj, ij, h, tempmu_num_jj[h-1],tempmu_den_jj[h-1]);
        
	                      } 
	                      mu[(jj-1)*(*g)+ij]=tempmu_num_sum/tempmu_den_sum;
     	                  //printf("%d %d %lf %lf %lf \n",jj,ij,tempmu_num_sum,tempmu_den_sum, mu[(jj-1)*(*g)+ij]);
                      }
                  }

		  
		/*M step for sigma*/	
                  for(jj=1;jj<=(*p);jj++){
                      for(ij=0;ij<(*g);ij++){
                		  tempsigma_num_sum=0;
		                  tempmu_den_sum=0;
	                      for(h=1;h<=(*n);h++){
                      	      for(i=0;i<(*p);i++){
	                              id1[i]=0;
	                              id2[i]=0;
                              }      
                              tempsigma_num_jj[h-1]=nest_of_sigma_jj_num_k(n, g, p, &h, j, &jj, &ij, id1, uhat, pi, x, mu, sigma);
                              tempmu_den_jj[h-1]=nest_of_mu_jj_den_k(n, g, p, &h, j, &jj, &ij, id2, uhat, pi, x, mu, sigma);
							  tempsigma_num_sum=tempsigma_num_sum+tempsigma_num_jj[h-1];
                              tempmu_den_sum=tempmu_den_sum+tempmu_den_jj[h-1];
//                              printf("%d %d %d %lf %lf\n", jj, ij, h, tempsigma_num_jj[h-1],tempmu_den_jj[h-1]);
        
	                      } 
	                      sigma[(jj-1)*(*g)+ij]=tempsigma_num_sum/tempmu_den_sum;
                          //printf("%d %d %lf %lf %lf \n",jj,ij,tempsigma_num_sum,tempmu_den_sum,sigma[(jj-1)*(*g)+ij]);
                      }
                  }	                
                
        /*constrain*/
        for(jj=0;jj<(*p);jj++){
        for(ij=0;ij<(*g);ij++){
            if(constrain[jj*(*g)+ij]==0-1){
                if(mu[jj*(*g)+ij]>0.0){
                    mu[jj*(*g)+ij]=0.0;
                }  
            }
            if(constrain[jj*(*g)+ij]==0){
                mu[jj*(*g)+ij] = 0.0;
                sigma[jj*(*g)+ij] = 1.0;
            }
            if(constrain[jj*(*g)+ij]==0+1){
                if(mu[jj*(*g)+ij]<0.0){
                    mu[jj*(*g)+ij]=0.0;
                } 
            }
        } 
        }
        
        /*check convergence*/
       loglik2 = (*loglik);
        (*loglik) = 0;
                  for(h=1;h<=(*n);h++){
          	          for(i=0;i<(*p);i++){
                          id[i]=0;
                      }
                      fpcd[h-1]=fpcdpart1_k(n, g, p, &h, pi,alpha, x, mu, sigma)+nest_of_fpcdpart2_k(n, g, p, &h, j, id, x, mu, sigma);
	                  logLg[h-1] = log(fpcd[h-1]); 
	                  (*loglik)= (*loglik) + logLg[h-1];
//                      printf("%d %lf %lf \n", h, fpcd[h-1], logLg[h-1]);
                  }
        if(fabs((*loglik)-loglik2)<(*tol)){
            flag = 0;
        }
        
        iter = iter + 1;        
    }
    (*convergence) = 0;
    if(iter<(*iteration)){
        (*convergence) = 1;
    }
   
    for (l=0;l<(*g);l++){
	       printf("%d %lf \n", l, pi[l]);
    }
    printf("%lf \n", *alpha);
    for (l=0;l<(*g)*(*p);l++){
	       printf("%d %lf %lf \n", l, mu[l],sigma[l]);
    }
}
/*no restriction on the null*/
else{
   /*initialization*/
    flag=1;
    (*loglik) = 0;
        for(h=1;h<=(*n);h++){
          	    for(i=0;i<(*p);i++){
                    id[i]=0;
                }
	             fpcd[h-1]=fpcdpart1_k(n, g, p, &h, pi,alpha, x, mu, sigma)+nest_of_fpcdpart2_k(n, g, p, &h, j, id, x, mu, sigma);
	             logLg[h-1] = log(fpcd[h-1]); 
	            (*loglik)= (*loglik) + logLg[h-1];
                printf("%d %lf %lf \n", h, fpcd[h-1], logLg[h-1]);
  	     } 
        /*printf("%lf \n", *loglik);/
    /*iteration*/
    iter = 0;
    while(flag>0 & iter<(*iteration)){
       	 /*E step for z*/
      	    	  for(h=1;h<=(*n);h++){
                      for(i=0;i<(*p);i++){
                          id[i]=0;
                      }
                      nest_ez_k(n, g, p, &h, j, id, uhat, pi, alpha, x, mu, sigma);
             	  }
                  for (l=0;l<power*(*n);l++){
//	                   printf("%d %lf", l, uhat[l]);
                  }
                  for(h=1;h<=(*n);h++){
                      for (l=0;l<(*g);l++){
                  	       for (jj=0;jj<(*p);jj++){
   	                  	       id[jj]=l;
   	                       } 
   	                       uhat[(h-1)*power+index_of(g,p,id)]=pi[l]*prod_of(n, g, p, &h, p, id, x, mu, sigma)/fpcd[h-1];
                       }
                  }
    	
        /*M step for pi*/
                 temppi_den=0;
    	          for(l=0;l<power*(*n);l++){      	    	
   	                  temppi_den=temppi_den+uhat[l];
          	      }

                  for(l=0;l<(*g);l++){
  	                  pi[l]=0;
  	                  for(jj=0;jj<(*p);jj++){
                  	      id[jj]=l;
                      }
	                  for(h=0;h<(*n);h++){
                      	pi[l]=pi[l]+uhat[h*power+index_of(g,p,id)];
	        
	                  }
                      pi[l]=pi[l]/temppi_den;
                      printf("%d %lf \n", l ,pi[l]);
        	      }  
        	      
        	      
        /*M step for alpha*/
                  sumpi=0;
                  for(l=0;l<(*g);l++){
                  	  sumpi=sumpi+pi[l];
                  }
		          (*alpha)=(1.0-sumpi)/(power-(*g));
	
                     
                  /*M step for mu*/		

                  for(jj=1;jj<=(*p);jj++){
                      for(ij=0;ij<(*g);ij++){
                          fflush(stdin);
		                  tempmu_num_sum=0;
               	          tempmu_den_sum=0;	
	                      for(h=1;h<=(*n);h++){
		                      for(ii=0;ii<(*p);ii++){
	                              id1[ii]=0;
	                              id2[ii]=0;
                              }      
                              tempmu_num_jj[h-1]=nest_of_mu_jj_num_k(n, g, p, &h, j, &jj, &ij, id1, uhat, pi, x, mu, sigma);
                              tempmu_den_jj[h-1]=nest_of_mu_jj_den_k(n, g, p, &h, j, &jj, &ij, id2, uhat, pi, x, mu, sigma);
                              tempmu_num_sum=tempmu_num_sum+tempmu_num_jj[h-1];
                              tempmu_den_sum=tempmu_den_sum+tempmu_den_jj[h-1];
                              printf("%d %d %d %lf %lf\n", jj, ij, h, tempmu_num_jj[h-1],tempmu_den_jj[h-1]);
        
	                      } 
	                      mu[(jj-1)*(*g)+ij]=tempmu_num_sum/tempmu_den_sum;
     	                  printf("%d %d %lf %lf %lf \n",jj,ij,tempmu_num_sum,tempmu_den_sum, mu[(jj-1)*(*g)+ij]);
                      }
                  }
                      
                  /*M step for sigma*/	
                  for(jj=1;jj<=(*p);jj++){
                      for(ij=0;ij<(*g);ij++){
                		  tempsigma_num_sum=0;
		                  tempmu_den_sum=0;
	                      for(h=1;h<=(*n);h++){
                      	      for(i=0;i<(*p);i++){
	                              id1[i]=0;
	                              id2[i]=0;
                              }      
                              tempsigma_num_jj[h-1]=nest_of_sigma_jj_num_k(n, g, p, &h, j, &jj, &ij, id1, uhat, pi, x, mu, sigma);
                              tempmu_den_jj[h-1]=nest_of_mu_jj_den_k(n, g, p, &h, j, &jj, &ij, id2, uhat, pi, x, mu, sigma);  
							  tempsigma_num_sum=tempsigma_num_sum+tempsigma_num_jj[h-1];
                              tempmu_den_sum=tempmu_den_sum+tempmu_den_jj[h-1];
                              printf("%d %d %d %lf %lf\n", jj, ij, h, tempsigma_num_jj[h-1],tempmu_den_jj[h-1]);
        
	                      } 
	                      sigma[(jj-1)*(*g)+ij]=tempsigma_num_sum/tempmu_den_sum;
                          printf("%d %d %lf %lf %lf \n",jj,ij,tempsigma_num_sum,tempmu_den_sum,sigma[(jj-1)*(*g)+ij]);
                      }
                  }
 
        /*check convergence*/
       loglik2 = (*loglik);
        (*loglik) = 0;
        for(h=1;h<=(*n);h++){
          	          for(i=0;i<(*p);i++){
                          id[i]=0;
                      }
                      fpcd[h-1]=fpcdpart1_k(n, g, p, &h, pi,alpha, x, mu, sigma)+nest_of_fpcdpart2_k(n, g, p, &h, j, id, x, mu, sigma);
	                  logLg[h-1] = log(fpcd[h-1]); 
	                  (*loglik)= (*loglik) + logLg[h-1];
                      printf("%d %lf %lf \n", h, fpcd[h-1], logLg[h-1]);
                  }
                  
        if(fabs((*loglik)-loglik2)<(*tol)){
            flag = 0;
        }
        
        iter = iter + 1;        
    }
    (*convergence) = 0;
    if(iter<(*iteration)){
        (*convergence) = 1;
    }
   
    for (l=0;l<(*g);l++){
	       printf("%d %lf \n", l, pi[l]);
    }
    printf("%lf \n", *alpha);
    for (l=0;l<(*g)*(*p);l++){
	       printf("%d %lf %lf \n", l, mu[l],sigma[l]);
    }
   
    
}
free(id1);
free(id2);
free(tempmu_num_jj);
free(tempmu_den_jj);
free(tempsigma_num_jj);
free(fpcd);
free(logLg);
}  




void em_pdim_pcd_exchange_me(long *n, long *g, long *p, long *j, long *id,  double *uhat, double *pi, double *alpha, double *x, double *mu, double *sigma, double *loglik,double *tol, long *restriction, long *constrain, long *iteration, long *convergence){
     long h,l,jj,ij,jjj,ijj,idx,*id1,*id2,power;
     int i, ii,flag, iter;
     double *Lg,*logLg,loglik2,sumpi,temppi_den, tempmu_num_sum, tempmu_den_sum, tempsigma_num_sum, *tempmu_num_jj,*tempmu_den_jj,*tempsigma_num_jj, *fpcd;
     power=pow((*g),(*p));
	 id1 = (long *) malloc((*p) * sizeof(long));
     id2 = (long *) malloc((*p) * sizeof(long));
     tempmu_num_jj= (double *) malloc((*n) * sizeof(double));
     tempmu_den_jj = (double *) malloc((*n)* sizeof(double));
     tempsigma_num_jj = (double *) malloc((*n) * sizeof(double));
     fpcd = (double *) malloc((*n) * sizeof(double));
     logLg = (double *) malloc((*n) * sizeof(double));


/*theoretical restriction on the null*/
if((*restriction)>0){    
    /*initialization*/
    flag=1;
    (*loglik) = 0;
    	for(h=1;h<=(*n);h++){
          	    for(i=0;i<(*p);i++){
                    id[i]=0;
                }
	            fpcd[h-1]=fpcdpart1_k(n, g, p, &h, pi, alpha, x, mu, sigma)+nest_of_fpcdpart2_k(n, g, p, &h, j, id, x, mu, sigma);
	            logLg[h-1] = log(fpcd[h-1]); 
	            (*loglik)= (*loglik) + logLg[h-1];
//                printf("%d %lf %lf \n", h, fpcd[h-1], logLg[h-1]);
   	    }
           
        
        printf("%lf \n", *loglik);
   
    /*iteration*/
    iter = 0;
    while(flag>0 & iter<(*iteration)){
    	printf("%d \n",iter);
    	        /*M step for pi*/
                 temppi_den=0;
    	          for(l=0;l<power*(*n);l++){      	    	
   	                  temppi_den=temppi_den+uhat[l];
          	      }

                  for(l=0;l<(*g);l++){
  	                  pi[l]=0;
  	                  for(jj=0;jj<(*p);jj++){
                  	      id[jj]=l;
                      }
	                  for(h=0;h<(*n);h++){
                      	pi[l]=pi[l]+uhat[h*power+index_of(g,p,id)];
	        
	                  }
                      pi[l]=pi[l]/temppi_den;
                      printf("%d %lf \n", l ,pi[l]);
        	      }  
  
        /*M step for alpha*/
                  sumpi=0;
                  for(l=0;l<(*g);l++){
                  	 sumpi=sumpi+pi[l];
                  }
                  printf("%lf \n", sumpi);
                  
		          (*alpha)=(1.0-sumpi)/(power-(*g));
		          printf("%lf \n", *alpha);
     
	    /*M step for mu*/		

                  for(jj=1;jj<=(*p);jj++){
                      for(ij=0;ij<(*g);ij++){
                          fflush(stdin);
		                  tempmu_num_sum=0;
               	          tempmu_den_sum=0;	
	                      for(h=1;h<=(*n);h++){
		                      for(ii=0;ii<(*p);ii++){
	                              id1[ii]=0;
	                              id2[ii]=0;
                              }      
                              tempmu_num_jj[h-1]=nest_of_mu_jj_num_k(n, g, p, &h, j, &jj, &ij, id1, uhat, pi, x, mu, sigma);
                              tempmu_den_jj[h-1]=nest_of_mu_jj_den_k(n, g, p, &h, j, &jj, &ij, id2, uhat, pi, x, mu, sigma);
                              tempmu_num_sum=tempmu_num_sum+tempmu_num_jj[h-1];
                              tempmu_den_sum=tempmu_den_sum+tempmu_den_jj[h-1];
//                              printf("%d %d %d %lf %lf\n", jj, ij, h, tempmu_num_jj[h-1],tempmu_den_jj[h-1]);
        
	                      } 
	                      mu[(jj-1)*(*g)+ij]=tempmu_num_sum/tempmu_den_sum;
     	                  //printf("%d %d %lf %lf %lf \n",jj,ij,tempmu_num_sum,tempmu_den_sum, mu[(jj-1)*(*g)+ij]);
                      }
                  }

		  
		/*M step for sigma*/	
                  for(jj=1;jj<=(*p);jj++){
                      for(ij=0;ij<(*g);ij++){
                		  tempsigma_num_sum=0;
		                  tempmu_den_sum=0;
	                      for(h=1;h<=(*n);h++){
                      	      for(i=0;i<(*p);i++){
	                              id1[i]=0;
	                              id2[i]=0;
                              }      
                              tempsigma_num_jj[h-1]=nest_of_sigma_jj_num_k(n, g, p, &h, j, &jj, &ij, id1, uhat, pi, x, mu, sigma);
                              tempmu_den_jj[h-1]=nest_of_mu_jj_den_k(n, g, p, &h, j, &jj, &ij, id2, uhat, pi, x, mu, sigma);
							  tempsigma_num_sum=tempsigma_num_sum+tempsigma_num_jj[h-1];
                              tempmu_den_sum=tempmu_den_sum+tempmu_den_jj[h-1];
//                              printf("%d %d %d %lf %lf\n", jj, ij, h, tempsigma_num_jj[h-1],tempmu_den_jj[h-1]);
        
	                      } 
	                      sigma[(jj-1)*(*g)+ij]=tempsigma_num_sum/tempmu_den_sum;
                          //printf("%d %d %lf %lf %lf \n",jj,ij,tempsigma_num_sum,tempmu_den_sum,sigma[(jj-1)*(*g)+ij]);
                      }
                  }	
        /*E step for z*/
      	    	  for(h=1;h<=(*n);h++){
                      for(i=0;i<(*p);i++){
                          id[i]=0;
                      }
                      nest_ez_k(n, g, p, &h, j, id, uhat, pi,alpha, x, mu, sigma);
             	  }
                  for (l=0;l<power*(*n);l++){
//	                   printf("%d %lf", l, uhat[l]);
                  }
                  for(h=1;h<=(*n);h++){
                      for (l=0;l<(*g);l++){
                  	       for (jj=0;jj<(*p);jj++){
   	                  	       id[jj]=l;
   	                       } 
   	                       uhat[(h-1)*power+index_of(g,p,id)]=pi[l]*uhat[(h-1)*power+index_of(g,p,id)]/(*alpha);
							  //prod_of(n, g, p, &h, p, id, x, mu, sigma)/fpcd[h-1];
                       }
                  }
	
                
                
        /*constrain*/
        for(jj=0;jj<(*p);jj++){
        for(ij=0;ij<(*g);ij++){
            if(constrain[jj*(*g)+ij]==0-1){
                if(mu[jj*(*g)+ij]>0.0){
                    mu[jj*(*g)+ij]=0.0;
                }  
            }
            if(constrain[jj*(*g)+ij]==0){
                mu[jj*(*g)+ij] = 0.0;
                sigma[jj*(*g)+ij] = 1.0;
            }
            if(constrain[jj*(*g)+ij]==0+1){
                if(mu[jj*(*g)+ij]<0.0){
                    mu[jj*(*g)+ij]=0.0;
                } 
            }
        } 
        }
        
        /*check convergence*/
       loglik2 = (*loglik);
        (*loglik) = 0;
                  for(h=1;h<=(*n);h++){
          	          for(i=0;i<(*p);i++){
                          id[i]=0;
                      }
                      fpcd[h-1]=fpcdpart1_k(n, g, p, &h, pi,alpha, x, mu, sigma)+nest_of_fpcdpart2_k(n, g, p, &h, j, id, x, mu, sigma);
	                  logLg[h-1] = log(fpcd[h-1]); 
	                  (*loglik)= (*loglik) + logLg[h-1];
//                      printf("%d %lf %lf \n", h, fpcd[h-1], logLg[h-1]);
                  }
        if(fabs((*loglik)-loglik2)<(*tol)){
            flag = 0;
        }
        
        iter = iter + 1;        
    }
    (*convergence) = 0;
    if(iter<(*iteration)){
        (*convergence) = 1;
    }
   
    for (l=0;l<(*g);l++){
	       printf("%d %lf \n", l, pi[l]);
    }
    printf("%lf \n", *alpha);
    for (l=0;l<(*g)*(*p);l++){
	       printf("%d %lf %lf \n", l, mu[l],sigma[l]);
    }
}
/*no restriction on the null*/
else{
   /*initialization*/
    flag=1;
    (*loglik) = 0;
        for(h=1;h<=(*n);h++){
          	    for(i=0;i<(*p);i++){
                    id[i]=0;
                }
	             fpcd[h-1]=fpcdpart1_k(n, g, p, &h, pi,alpha, x, mu, sigma)+nest_of_fpcdpart2_k(n, g, p, &h, j, id, x, mu, sigma);
	             logLg[h-1] = log(fpcd[h-1]); 
	            (*loglik)= (*loglik) + logLg[h-1];
                printf("%d %lf %lf \n", h, fpcd[h-1], logLg[h-1]);
  	     } 
        /*printf("%lf \n", *loglik);/
    /*iteration*/
    iter = 0;
    while(flag>0 & iter<(*iteration)){
    	
    	
    	
    	    	
        /*M step for pi*/
                 temppi_den=0;
    	          for(l=0;l<power*(*n);l++){      	    	
   	                  temppi_den=temppi_den+uhat[l];
          	      }

                  for(l=0;l<(*g);l++){
  	                  pi[l]=0;
  	                  for(jj=0;jj<(*p);jj++){
                  	      id[jj]=l;
                      }
	                  for(h=0;h<(*n);h++){
                      	pi[l]=pi[l]+uhat[h*power+index_of(g,p,id)];
	        
	                  }
                      pi[l]=pi[l]/temppi_den;
                      printf("%d %lf \n", l ,pi[l]);
        	      }  
        	      
        	      
        /*M step for alpha*/
                  sumpi=0;
                  for(l=0;l<(*g);l++){
                  	  sumpi=sumpi+pi[l];
                  }
		          (*alpha)=(1.0-sumpi)/(power-(*g));
	
                     
                  /*M step for mu*/		

                  for(jj=1;jj<=(*p);jj++){
                      for(ij=0;ij<(*g);ij++){
                          fflush(stdin);
		                  tempmu_num_sum=0;
               	          tempmu_den_sum=0;	
	                      for(h=1;h<=(*n);h++){
		                      for(ii=0;ii<(*p);ii++){
	                              id1[ii]=0;
	                              id2[ii]=0;
                              }      
                              tempmu_num_jj[h-1]=nest_of_mu_jj_num_k(n, g, p, &h, j, &jj, &ij, id1, uhat, pi, x, mu, sigma);
                              tempmu_den_jj[h-1]=nest_of_mu_jj_den_k(n, g, p, &h, j, &jj, &ij, id2, uhat, pi, x, mu, sigma);
                              tempmu_num_sum=tempmu_num_sum+tempmu_num_jj[h-1];
                              tempmu_den_sum=tempmu_den_sum+tempmu_den_jj[h-1];
                              printf("%d %d %d %lf %lf\n", jj, ij, h, tempmu_num_jj[h-1],tempmu_den_jj[h-1]);
        
	                      } 
	                      mu[(jj-1)*(*g)+ij]=tempmu_num_sum/tempmu_den_sum;
     	                  printf("%d %d %lf %lf %lf \n",jj,ij,tempmu_num_sum,tempmu_den_sum, mu[(jj-1)*(*g)+ij]);
                      }
                  }
                      
                  /*M step for sigma*/	
                  for(jj=1;jj<=(*p);jj++){
                      for(ij=0;ij<(*g);ij++){
                		  tempsigma_num_sum=0;
		                  tempmu_den_sum=0;
	                      for(h=1;h<=(*n);h++){
                      	      for(i=0;i<(*p);i++){
	                              id1[i]=0;
	                              id2[i]=0;
                              }      
                              tempsigma_num_jj[h-1]=nest_of_sigma_jj_num_k(n, g, p, &h, j, &jj, &ij, id1, uhat, pi, x, mu, sigma);
                              tempmu_den_jj[h-1]=nest_of_mu_jj_den_k(n, g, p, &h, j, &jj, &ij, id2, uhat, pi, x, mu, sigma);  
							  tempsigma_num_sum=tempsigma_num_sum+tempsigma_num_jj[h-1];
                              tempmu_den_sum=tempmu_den_sum+tempmu_den_jj[h-1];
                              printf("%d %d %d %lf %lf\n", jj, ij, h, tempsigma_num_jj[h-1],tempmu_den_jj[h-1]);
        
	                      } 
	                      sigma[(jj-1)*(*g)+ij]=tempsigma_num_sum/tempmu_den_sum;
                          printf("%d %d %lf %lf %lf \n",jj,ij,tempsigma_num_sum,tempmu_den_sum,sigma[(jj-1)*(*g)+ij]);
                      }
                  }
       	 /*E step for z*/
      	    	  for(h=1;h<=(*n);h++){
                      for(i=0;i<(*p);i++){
                          id[i]=0;
                      }
                      nest_ez_k(n, g, p, &h, j, id, uhat, pi, alpha, x, mu, sigma);
             	  }
                  for (l=0;l<power*(*n);l++){
//	                   printf("%d %lf", l, uhat[l]);
                  }
                  for(h=1;h<=(*n);h++){
                      for (l=0;l<(*g);l++){
                  	       for (jj=0;jj<(*p);jj++){
   	                  	       id[jj]=l;
   	                       } 
   	                       uhat[(h-1)*power+index_of(g,p,id)]=pi[l]*prod_of(n, g, p, &h, p, id, x, mu, sigma)/fpcd[h-1];
                       }
                  }

 
        /*check convergence*/
       loglik2 = (*loglik);
        (*loglik) = 0;
        for(h=1;h<=(*n);h++){
          	          for(i=0;i<(*p);i++){
                          id[i]=0;
                      }
                      fpcd[h-1]=fpcdpart1_k(n, g, p, &h, pi,alpha, x, mu, sigma)+nest_of_fpcdpart2_k(n, g, p, &h, j, id, x, mu, sigma);
	                  logLg[h-1] = log(fpcd[h-1]); 
	                  (*loglik)= (*loglik) + logLg[h-1];
                      printf("%d %lf %lf \n", h, fpcd[h-1], logLg[h-1]);
                  }
                  
        if(fabs((*loglik)-loglik2)<(*tol)){
            flag = 0;
        }
        
        iter = iter + 1;        
    }
    (*convergence) = 0;
    if(iter<(*iteration)){
        (*convergence) = 1;
    }
   
    for (l=0;l<(*g);l++){
	       printf("%d %lf \n", l, pi[l]);
    }
    printf("%lf \n", *alpha);
    for (l=0;l<(*g)*(*p);l++){
	       printf("%d %lf %lf \n", l, mu[l],sigma[l]);
    }
   
    
}
free(id1);
free(id2);
free(tempmu_num_jj);
free(tempmu_den_jj);
free(tempsigma_num_jj);
free(fpcd);
free(logLg);
}