automix.c

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

    used. Note that if the parameters are unavailable or inconsistent with
    the user supplied functions or unavailable go straight to section 5.2
    where initial within-model RWM are performed */
   
   if(mode>2){
      printf("\nInvalid mode entered. Mode must be 0,1,2");
      return -100;
   }
   else if(mode==1){
      if((check=fscanf(fpmix,"%d",&k1))==EOF){
         printf("\nEnd of file encountered before parameters read:");
         printf("\nContinuing using RWM to estimate parameters");
         mode=0;
         goto RWMSTART;
      }
      if(k1!=kmax){
         printf("\nFile kmax contradicts getkmax function:");
         printf("\nContinuing using RWM to estimate parameters");
         mode=0;
         goto RWMSTART;
      }
      for(k1=0;k1<kmax;k1++){
         if((check=fscanf(fpmix,"%d",&j1))==EOF){
            printf("\nEnd of file encountered before parameters read:");
            printf("\nContinuing using RWM to estimate parameters");
            mode=0;
            goto RWMSTART;
         }
         if(j1!=nk[k1]){
            printf("\nFile kmax contradicts getnk function:");
            printf("\nContinuing using RWM to estimate parameters");
            mode=0;
            goto RWMSTART;
         }
      }
      for(k1=0;k1<kmax;k1++){
         nkk=nk[k1];
         for(j1=0;j1<nkk;j1++){
            if((check=fscanf(fpmix,"%lf",&(sig[k1][j1])))==EOF){
               printf("\nEnd of file encountered before parameters read:");
               printf("\nContinuing using RWM to estimate parameters");
               mode=0;
               goto RWMSTART;
            }
         }
         if((check=fscanf(fpmix,"%d",&(Lk[k1])))==EOF){
            printf("\nEnd of file encountered before parameters read:");
            printf("\nContinuing using RWM to estimate parameters");
            mode=0;
            goto RWMSTART;
         }
         Lkk=Lk[k1];
         for(l1=0;l1<Lkk;l1++){
            if((check=fscanf(fpmix,"%lf",&(lambda[k1][l1])))==EOF){
               printf("\nEnd of file encountered before parameters read:");
               printf("\nContinuing using RWM to estimate parameters");
               mode=0;
               goto RWMSTART;
            }
            for(j1=0;j1<nkk;j1++){
               if((check=fscanf(fpmix,"%lf",&(mu[k1][l1][j1])))==EOF){
                  printf("\nEnd of file encountered before parameters read:");
                  printf("\nContinuing using RWM to estimate parameters");
                  mode=0;
                  goto RWMSTART;
               }
            }
            for(j1=0;j1<nkk;j1++){
               for(j2=0;j2<=j1;j2++){
                  if((check=fscanf(fpmix,"%lf",&(B[k1][l1][j1][j2])))==EOF){
                     printf("\nEnd of file encountered before parameters read:");
                     printf("\nContinuing using RWM to estimate parameters");
                     mode=0;
                     goto RWMSTART;
                  }
               }
            }
         }
         sumlambda=0.0;
         for(l1=0;l1<Lkk;l1++){
            sumlambda+=lambda[k1][l1];
         }
         if(sumlambda<0.99999||sumlambda>1.00001){
            printf("\nComponents weights read do not sum to one for k=%d:",k1);
            printf("\nContinuing using RWM to estimate parameters");
            mode=0;
            goto RWMSTART;
         }
         if(sumlambda<1.0||sumlambda>1.0){
            for(l1=0;l1<Lkk;l1++){
               lambda[k1][l1]/=sumlambda;
            }
         }
      }
      if(!(fpmix==NULL)){
         fclose(fpmix);
      }
      
   }
   else if(mode==0||mode==2){
      
  /* --- Section 5.2 - Within-model runs if mixture parameters unavailable -*/
      
      RWMSTART:
         
         for(k1=0;k1<kmax;k1++){
     /* --- Section 5.2.1 - RWM Within Model (Stage 1) -------*/
            
            nkk=nk[k1];
            nparams=nkk+(nkk*(nkk+1))/2;
            lendata=1000*nkk;
            nsweepr=max(nsweep2,10000*nkk);
            nburn=nsweepr/10;
            nsweepr+=nburn;
            data=(double**) malloc(lendata*sizeof(double));
            for(i1=0;i1<lendata;i1++){
               data[i1]=(double*)malloc(nkk*sizeof(double));
            }
            rwm=(double*)malloc(nkk*sizeof(double));
            rwmn=(double*)malloc(nkk*sizeof(double));
            nacc=(int*)malloc(nkk*sizeof(int));
            ntry=(int*)malloc(nkk*sizeof(int));
            Znkk=(double*)malloc(nkk*sizeof(double));
            init=(int*)malloc(nkk*sizeof(int));
            
            printf("\nRWM for Model %d",k1+1);
            fprintf(fpl,"\nRWM for Model %d",k1+1);
            fprintf(fpcf,"RWM for Model %d\n",k1+1);
            fprintf(fpad,"RWM for Model %d\n",k1+1);
            fflush(NULL);
            getic(k1,nkk,rwm);
            for(j1=0;j1<nkk;j1++){
               rwmn[j1]=rwm[j1];
               sig[k1][j1]=10.0;
               nacc[j1]=0;
               ntry[j1]=0;
            }
            lp=lpost(k1,nkk,rwm,&llh);
            
            i2=0;
            remain=nsweepr;
            for(sweep=1;sweep<=nsweepr;sweep++){
               remain--;
               if((sweep>=nburn)&&(fmod((sweep-nburn),((nsweepr-nburn)/10))<tol)){
                  printf("\nNo. of iterations remaining: %d",remain);
                  fflush(NULL);
               }
               u=sdrand();
               if(sweep>nburn&&u<0.1){
                  if(dof>0){
                     rt(Znkk,nkk,dof);
                  }
                  else{
                     gauss(Znkk,nkk);
                  }
                  for(j1=0;j1<nkk;j1++){
                     rwmn[j1]=rwm[j1]+sig[k1][j1]*Znkk[j1];
                  }
                  lpn=lpost(k1,nkk,rwmn,&llhn);
                  if(sdrand()<exp(max(-30.0,min(0.0,lpn-lp)))){
                     for(j1=0;j1<nkk;j1++){
                        rwm[j1]=rwmn[j1];
                     }
                     lp=lpn;
                     llh=llhn;
                  }
               }
               else{
                  gamma=10.0*pow(1.0/(sweep+1),2.0/3.0);
                  for(j1=0;j1<nkk;j1++){
                     rwmn[j1]=rwm[j1];
                  }
                  for(j1=0;j1<nkk;j1++){
                     if(dof>0){
                        rt(Z,1,dof);
                     }
                     else{
                        gauss(Z,1);
                     }
                     rwmn[j1]=rwm[j1]+sig[k1][j1]*Z[0];
                     lpn=lpost(k1,nkk,rwmn,&llhn);
                     accept=min(1,exp(max(-30.0,min(0.0,lpn-lp))));
                     if(sdrand()<accept){
                        (nacc[j1])++;
                        (ntry[j1])++;
                        rwm[j1]=rwmn[j1];
                        lp=lpn;
                        llh=llhn;
                        sig[k1][j1]=max(0,sig[k1][j1]-gamma*(alphastar-1));
                     }
                     else{
                        (ntry[j1])++;
                        rwmn[j1]=rwm[j1];
                        sig[k1][j1]=max(0,sig[k1][j1]-gamma*(alphastar));
                     }
                  }
               }
               if(remain<(10000*nkk)&&fmod(remain,10.0)<0.05){
                  for(j1=0;j1<nkk;j1++){
                     data[i2][j1]=rwm[j1];
                  }
                  i2++;
               }
               if(fmod(sweep,100.0)<0.05){
                  for(j1=0;j1<nkk;j1++){
                     fprintf(fpad,"%lf %lf ",sig[k1][j1],
                     (double)nacc[j1]/(double)ntry[j1]);
                  }
                  fprintf(fpad,"\n");
               }
               
            }
            free(init);
            free(rwm);
            free(rwmn);
            free(nacc);
            free(ntry);
            free(Znkk);
            
     /* --- Section 5.2.2 - Fit Mixture to within-model sample, (stage 2)- */
     /* Note only done if mode 0 (m=0) if mode m=2, go to section 5.2.3*/
     /* Mixture fitting done component wise EM algorithm described in
    Figueiredo and Jain, 2002 (see thesis for full reference) */
            
            printf("\nMixture Fitting: Model %d",k1+1);
            if(mode==0){
               Lkk=Lkmaxmax;
               init=(int*)malloc(Lkk*sizeof(int));
               l1=0;
               while(l1<Lkk){
                  indic=0;
                  u=sdrand();
                  init[l1]=(int)floor(lendata*u);
                  if(l1>0){
                     for(l2=0;l2<l1;l2++){
                        if(init[l2]==init[l1]){
                           indic=1;
                           break;
                        }
                     }
                  }
                  if(indic==0){
                     l1++;
                  }
               }
               
               datamean=(double*)malloc(nkk*sizeof(double));
               M1=(double**)malloc(nkk*sizeof(double));
               for(j1=0;j1<nkk;j1++){
                  M1[j1]=(double*)malloc(nkk*sizeof(double));
               }
               for(j1=0;j1<nkk;j1++){
                  datamean[j1]=0.0;
                  for(i1=0;i1<lendata;i1++){
                     datamean[j1]+=data[i1][j1];
                  }
                  datamean[j1]/=((double)lendata);
               }
               for(j1=0;j1<nkk;j1++){
                  for(j2=0;j2<nkk;j2++){
                     M1[j1][j2]=0;
                     for(i1=0;i1<lendata;i1++){
                        M1[j1][j2]+=(data[i1][j1]-datamean[j1])*
                        (data[i1][j2]-datamean[j2]);
                     }
                     M1[j1][j2]/=((double)lendata);
                  }
               }
               sigma=0.0;
               for(j1=0;j1<nkk;j1++){
                  sigma+=M1[j1][j1];
               }
               sigma/=(10.0*nkk);
               
               for(l1=0;l1<Lkk;l1++){
                  for(j1=0;j1<nkk;j1++){
                     mu[k1][l1][j1]=data[init[l1]][j1];
                     BBT[k1][l1][j1][j1]=sigma;
                     B[k1][l1][j1][j1]=BBT[k1][l1][j1][j1];
                     for(j2=0;j2<j1;j2++){
                        BBT[k1][l1][j1][j2]=0.0;
                        B[k1][l1][j1][j2]=BBT[k1][l1][j1][j2];
                     }
                  }
                  chol(nkk,B[k1][l1]);
                  lambda[k1][l1]=1.0/Lkk;
               }
               
               w=(double**)malloc(lendata*sizeof(double));
               logw=(double*)malloc(Lkk*sizeof(double));
               lpdatagivenl=(double**)malloc(lendata*sizeof(double));
               for(i1=0;i1<lendata;i1++){
                  w[i1]=(double*)malloc(Lkk*sizeof(double));
                  lpdatagivenl[i1]=(double*)malloc(Lkk*sizeof(double));
               }
               
               for(i1=0;i1<lendata;i1++){
                  sum=0.0;
                  for(l1=0;l1<Lkk;l1++){
                     lpdatagivenl[i1][l1]=lnormprob(k1,nkk,l1,mu,B,data[i1]);
                     logw[l1]=log(lambda[k1][l1])+lpdatagivenl[i1][l1];
                     w[i1][l1]=exp(logw[l1]);
                     sum+=w[i1][l1];
                  }
                  for(l1=0;l1<Lkk;l1++){
                     w[i1][l1]/=sum;
                  }
               }
               
               sumw=(double*)malloc(Lkk*sizeof(double));
               
               stop=0;
               count=0;
               
               while(!stop){
                  count++;
                  l1=0;
                  natann=0;
                  forceann=0;
                  while(l1<Lkk){
                     sumwnew=0.0;
                     for(l2=0;l2<Lkk;l2++){
                        sumw[l2]=0.0;
                        for(i1=0;i1<lendata;i1++){
                           sumw[l2]+=w[i1][l2];
                        }
                        wnew=max(0.0,(sumw[l2]-nparams/2.0));
                        if(l2==l1){
                           wnewl1=wnew;
                        }
                        sumwnew+=wnew;
                     }
                     lambda[k1][l1]=wnewl1/sumwnew;
                     sumlambda=0.0;
                     for(l2=0;l2<Lkk;l2++){
                        sumlambda+=lambda[k1][l2];
                     }
                     for(l2=0;l2<Lkk;l2++){
                        lambda[k1][l2]/=sumlambda;
                     }
                     
                     if(lambda[k1][l1]>0.005){
        /*changed to 0.005 from 0.0 -renormalise else */
                        for(j1=0;j1<nkk;j1++){
                           mu[k1][l1][j1]=0.0;
                           for(i1=0;i1<lendata;i1++){
                              mu[k1][l1][j1]+=data[i1][j1]*w[i1][l1];
                           }
                           mu[k1][l1][j1]/=sumw[l1];
                           
                           for(j2=0;j2<=j1;j2++){
                              BBT[k1][l1][j1][j2]=0.0;
                              for(i1=0;i1<lendata;i1++){
                                 BBT[k1][l1][j1][j2]+=(data[i1][j1]-mu[k1][l1][j1])*
                                 (data[i1][j2]-mu[k1][l1][j2])*w[i1][l1];
                              }
                              BBT[k1][l1][j1][j2]/=sumw[l1];
                              B[k1][l1][j1][j2]=BBT[k1][l1][j1][j2];
                           }
                        }
                        
                        chol(nkk,B[k1][l1]);
                        
                        for(i1=0;i1<lendata;i1++){
                           lpdatagivenl[i1][l1]=lnormprob(k1,nkk,l1,mu,B,data[i1]);
                        }
                        l1++;
                        
                     }
                     else{
                        if(fmod(Lkk,5)<0.05){
                           printf("\n");
                        }
                        printf("%d(%d-n) ",Lkk,count);
                        natann=1;
                        if(l1<(Lkk-1)){
                           for(l2=l1;l2<(Lkk-1);l2++){
                              lambda[k1][l2]=lambda[k1][l2+1];
                              for(j1=0;j1<nkk;j1++){
                                 mu[k1][l2][j1]=mu[k1][l2+1][j1];
                                 for(j2=0;j2<=j1;j2++){
                                    BBT[k1][l2][j1][j2]=BBT[k1][l2+1][j1][j2];
                                    B[k1][l2][j1][j2]=B[k1][l2+1][j1][j2];
                                 }
                              }
                              for(i1=0;i1<lendata;i1++){
                                 lpdatagivenl[i1][l2]=lpdatagivenl[i1][l2+1];
                              }
                           }
                        }
                        Lkk--;
                        sumlambda=0.0;
                        for(l2=0;l2<Lkk;l2++){
                           sumlambda+=lambda[k1][l2];
                        }
                        for(l2=0;l2<Lkk;l2++){
                           lambda[k1][l2]/=sumlambda;
                        }
                        
                     }