automix.c

he AutoMix package is a C program for Unix-like systems, implementing the automatic reversible jum

                     lpn=0.0;
                     for(i1=0;i1<lendata;i1++){
                        sum=0.0;
                        for(l2=0;l2<Lkk;l2++){
                           logw[l2]=log(lambda[k1][l2])+lpdatagivenl[i1][l2];
                           w[i1][l2]=exp(logw[l2]);
                           sum+=w[i1][l2];
                        }
                        if(sum>0){
                           for(l2=0;l2<Lkk;l2++){
                              w[i1][l2]/=sum;
                           }
                           lpn+=log(sum);
                        }
                        else{
      /* if no component fits point well make equally likely */
                           for(l2=0;l2<Lkk;l2++){
                              w[i1][l2]=1.0/Lkk;
                           }
                           lpn+=(-500.0);
                        }
                     }
                  }
                  
                  sum=0.0;
                  for(l1=0;l1<Lkk;l1++){
                     sum+=log(lendata*lambda[k1][l1]/12.0);
                  }
                  costfnnew=(nparams/2.0)*sum+(Lkk/2.0)*log(lendata/12.0)+
                  Lkk*(nparams+1)/2.0-lpn;
                  
                  if(count==1){
                     costfn=costfnnew;
                  }
                  if(count==1||costfnnew<costfnmin){
                     Lkkmin=Lkk;
                     costfnmin=costfnnew;
                     for(l1=0;l1<Lkk;l1++){
                        lambdamin[k1][l1]=lambda[k1][l1];
                        for(j1=0;j1<nkk;j1++){
                           mumin[k1][l1][j1]=mu[k1][l1][j1];
                           for(j2=0;j2<=j1;j2++){
                              Bmin[k1][l1][j1][j2]=B[k1][l1][j1][j2];
                           }
                        }
                     }
                  }
                  if((fabs(costfn-costfnnew)<min(tol*fabs(costfn),0.01))&&(count>1)){
                     if(Lkk==1){
                        stop=1;
                     }
                     else{
                        if(fmod(Lkk,5)<0.05){
                           printf("\n");
                        }
                        printf("%d(%d-f) ",Lkk,count);
                        forceann=2;
                        minlambda=lambda[k1][0];
                        ldel=0;
                        for(l1=1;l1<Lkk;l1++){
                           if(minlambda>lambda[k1][l1]){
                              minlambda=lambda[k1][l1];
                              ldel=l1;
                           }
                        }
                        if(ldel<(Lkk-1)){
                           for(l1=ldel;l1<(Lkk-1);l1++){
                              lambda[k1][l1]=lambda[k1][l1+1];
                              for(j1=0;j1<nkk;j1++){
                                 mu[k1][l1][j1]=mu[k1][l1+1][j1];
                                 for(j2=0;j2<=j1;j2++){
                                    BBT[k1][l1][j1][j2]=BBT[k1][l1+1][j1][j2];
                                    B[k1][l1][j1][j2]=B[k1][l1+1][j1][j2];
                                 }
                              }
                              for(i1=0;i1<lendata;i1++){
                                 lpdatagivenl[i1][l1]=lpdatagivenl[i1][l1+1];
                              }
                           }
                        }
                        Lkk--;
                        sumlambda=0.0;
                        for(l1=0;l1<Lkk;l1++){
                           sumlambda+=lambda[k1][l1];
                        }
                        for(l1=0;l1<Lkk;l1++){
                           lambda[k1][l1]/=sumlambda;
                        }
                        
                        lpn=0.0;
                        for(i1=0;i1<lendata;i1++){
                           sum=0.0;
                           for(l2=0;l2<Lkk;l2++){
                              logw[l2]=log(lambda[k1][l2])+lpdatagivenl[i1][l2];
                              w[i1][l2]=exp(logw[l2]);
                              sum+=w[i1][l2];
                           }
                           if(sum>0){
                              for(l2=0;l2<Lkk;l2++){
                                 w[i1][l2]/=sum;
                              }
                              lpn+=log(sum);
                           }
                           else{
        /* if no component fits point well make equally likely */
                              for(l2=0;l2<Lkk;l2++){
                                 w[i1][l2]=1.0/Lkk;
                              }
                              lpn+=(-500.0);
                           }
                           
                        }
                        
                        sum=0.0;
                        for(l1=0;l1<Lkk;l1++){
                           sum+=log(lendata*lambda[k1][l1]/12.0);
                        }
                        costfnnew=(nparams/2.0)*sum+(Lkk/2.0)*log(lendata/12.0)+
                        Lkk*(nparams+1)/2.0-lpn;
                     }
                  }
                  if(count>5000){
                     stop=1;
                  }
                  costfn=costfnnew;
                  fprintf(fpcf,"%d %lf %lf %d\n",Lkk,lpn,costfnnew,(natann+forceann));
                  fflush(NULL);
               }
               
               for(j1=0;j1<nkk;j1++){
                  free(M1[j1]);
               }
               free(M1);
               free(datamean);
               for(i1=0;i1<lendata;i1++){
                  free(data[i1]);
                  free(w[i1]);
                  free(lpdatagivenl[i1]);
               }
               free(w);
               free(lpdatagivenl);
               free(data);
               free(logw);
               free(sumw);
               free(init);
               Lk[k1]=Lkkmin;
               for(l1=0;l1<Lkkmin;l1++){
                  lambda[k1][l1]=lambdamin[k1][l1];
                  for(j1=0;j1<nkk;j1++){
                     mu[k1][l1][j1]=mumin[k1][l1][j1];
                  }
                  for(j1=0;j1<nkk;j1++){
                     for(j2=0;j2<=j1;j2++){
                        B[k1][l1][j1][j2]=Bmin[k1][l1][j1][j2];
                     }
                  }
               }
            }
            else if(mode==2){
   /* --- Section 5.2.3 - Fit AutoRJ single mu vector and B matrix --*/
   /* Note only done if mode 2 (m=2).*/
               Lk[k1]=1;
               lambda[k1][0]=1.0;
               for(j1=0;j1<nkk;j1++){
                  mu[k1][0][j1]=0.0;
                  for(i1=0;i1<lendata;i1++){
                     mu[k1][0][j1]+=data[i1][j1];
                  }
                  mu[k1][0][j1]/=((double)lendata);
               }
               for(j1=0;j1<nkk;j1++){
                  for(j2=0;j2<=j1;j2++){
                     B[k1][0][j1][j2]=0.0;
                     for(i1=0;i1<lendata;i1++){
                        B[k1][0][j1][j2]+=(data[i1][j1]-mu[k1][0][j1])*
                        (data[i1][j2]-mu[k1][0][j2]);
                     }
                     B[k1][0][j1][j2]/=((double)(lendata-1));
                  }
               }
               chol(nkk,B[k1][0]);
               
               for(i1=0;i1<lendata;i1++){
                  free(data[i1]);
               }
               free(data);
            }
         }
   }
   
  /* Print mixture parameters to file (log and mix files) for reference
    and use in future runs. */
   
   sprintf(fname1,fname);
   strcat(fname1,"_mix.data");
   fpmix = fopen(fname1,"w");
   fprintf(fpmix,"%d\n",kmax);
   for(k1=0;k1<kmax;k1++){
      fprintf(fpmix,"%d\n",nk[k1]);
   }
   
   for(k1=0;k1<kmax;k1++){
      fprintf(fpl,"\nModel:%d\n",k1+1);
      Lkk=Lk[k1];
      nkk=nk[k1];
      fprintf(fpl,"\nARW params:\n");
      for(j1=0;j1<nkk;j1++){
         fprintf(fpl,"%lf ",sig[k1][j1]);
      }
      fprintf(fpl,"\n");
      fprintf(fpl,"\nLkk:%d\n",Lkk);
      for(j1=0;j1<nkk;j1++){
         fprintf(fpmix,"%lf\n",sig[k1][j1]);
      }
      fprintf(fpmix,"%d\n",Lkk);
      for(l1=0;l1<Lkk;l1++){
         fprintf(fpl,"\nComponent:%d\n",l1+1);
         fprintf(fpl,"lambda:%lf\n",lambda[k1][l1]);
         fprintf(fpmix,"%lf\n",lambda[k1][l1]);
         fprintf(fpl,"mu:\n");
         for(j1=0;j1<nkk;j1++){
            fprintf(fpl,"%lf ",mu[k1][l1][j1]);
            fprintf(fpmix,"%lf\n",mu[k1][l1][j1]);
         }
         fprintf(fpl,"\nB:\n");
         for(j1=0;j1<nkk;j1++){
            for(j2=0;j2<=j1;j2++){
               fprintf(fpl,"%lf ",B[k1][l1][j1][j2]);
               fprintf(fpmix,"%lf\n",B[k1][l1][j1][j2]);
            }
            fprintf(fpl,"\n");
         }
         detB[k1][l1]=det(k1,nkk,l1,B);
      }
   }
   fflush(NULL);
   
  /* --Section 6 - Secondary file handling -------------*/
   
   sprintf(fname1,fname);
   strcat(fname1,"_k");
   strcat(fname1,".data");
   fpk=fopen(fname1,"w");
   sprintf(fname1,fname);
   strcat(fname1,"_lp");
   strcat(fname1,".data");
   fplp=fopen(fname1,"w");
   
   for(k1=0;k1<kmax;k1++){
      sprintf(fname1,fname);
      sprintf(kno,"%d",k1+1);
      strcat(fname1,"_theta");
      strcat(fname1,kno);
      strcat(fname1,".data");
      fpt[k1]=fopen(fname1,"w");
   }
   
  /* --Section 7 - Final initialisation of variables ----*/
   
   naccrwmb=0;
   ntryrwmb=0;
   naccrwms=0;
   ntryrwms=0;
   nacctd=0;
   ntrytd=0;
   
   constt=100000.0;
   Lkmax=Lk[0];
   for(k1=1;k1<kmax;k1++){
      Lkmax=max(Lkmax,Lk[k1]);
   }
   k=(int)floor(kmax*sdrand());
   nkk=nk[k];
   Lkk=Lk[k];
   
   theta=(double*)malloc(nkmax*sizeof(double));
   thetan=(double*)malloc(nkmax*sizeof(double));
   work=(double*)malloc(nkmax*sizeof(double));
   palloc=(double*)malloc(Lkmax*sizeof(double));
   pallocn=(double*)malloc(Lkmax*sizeof(double));
   Znkk=(double*)malloc(nkmax*sizeof(double));
   propk=(double*)malloc(kmax*sizeof(double));
   pk=(double*)malloc(kmax*sizeof(double));
   
   getic(k,nkk,theta);
   
   lp=lpost(k,nkk,theta,&llh);
   
   for(k1=0;k1<kmax;k1++){
      pk[k1]=1.0/kmax;
      if(k1==k){
         propk[k1]=1.0;
      }
      else{
         propk[k1]=0.0;
      }
   }
   nreinit=1;
   reinit=0;
   pkllim=1.0/10.0;
   
   tol=0.5/nsweep;
   nburn=max(10000,(int)(nsweep/10));
   
   nsokal=1;
   nkeep=nsweep/(2*nsokal);
   nkeep=(int)pow(2.0,min(15,(int)(log(nkeep)/log(2.0)+0.001)));
   keep=nburn+(nsweep-nkeep*nsokal);
   xr=(double*)malloc(nkeep*sizeof(double));
   
   
   
 /* -----Start of main loop ----------------*/
   for(sweep=1;sweep<=(nburn+nsweep);sweep++){
      
   /* --Section 8 - RWM within-model moves ---*/
      
      
   /* Every 10 sweeps to block RWM */
      if(fmod(sweep,10)<0.05){
         ntryrwmb++;
         if(dof>0){
            rt(Znkk,nkk,dof);
         }
         else{
            gauss(Znkk,nkk);
         }
         for(j1=0;j1<nkk;j1++){
            thetan[j1]=theta[j1]+sig[k][j1]*Znkk[j1];
         }
         lpn=lpost(k,nkk,thetan,&llhn);
         if(sdrand()<exp(max(-30.0,min(0.0,lpn-lp)))){
            naccrwmb++;
            for(j1=0;j1<nkk;j1++){
               theta[j1]=thetan[j1];
            }
            lp=lpn;
            llh=llhn;
         }
      }
      else{
     /* else do component-wise RWM */
         for(j1=0;j1<nkk;j1++){
            thetan[j1]=theta[j1];
         }
         for(j1=0;j1<nkk;j1++){
            ntryrwms++;
            if(dof>0){
               rt(Z,1,dof);
            }
            else{
               gauss(Z,1);
            }
            thetan[j1]=theta[j1]+sig[k][j1]*Z[0];
            lpn=lpost(k,nkk,thetan,&llhn);
            if(sdrand()<exp(max(-30.0,min(0.0,lpn-lp)))){
               naccrwms++;
               theta[j1]=thetan[j1];
               lp=lpn;
               llh=llhn;
            }
            else{
               thetan[j1]=theta[j1];
            }
         }
      }
      
   /* --- Section 9 Reversible Jump Moves -------------------*/
      
   /* --Section 9.1 - Allocate current position to a component --*/
      
      ntrytd++;
      if(Lkk>1){
         sum=0.0;
         for(l1=0;l1<Lkk;l1++){
            palloc[l1]=log(lambda[k][l1])+lnormprob(k,nkk,l1,mu,B,theta);
            palloc[l1]=exp(palloc[l1]);
            sum+=palloc[l1];
         }
         if(sum>0){
            for(l1=0;l1<Lkk;l1++){
               palloc[l1]/=sum;
            }
         }
         else{
            for(l1=0;l1<Lkk;l1++){
               palloc[l1]=1.0/Lkk;
            }
         }
         u=sdrand();
         thresh=0.0;
         for(l1=0;l1<Lkk;l1++){
            thresh+=palloc[l1];
            if(u<thresh){
               l=l1;
               break;
            }
         }
      }
      else{
         l=0;