      program gibbs_dirichlet
*IN THIS PROGRAM WE IMPLEMENT THE conditional_mcmc_local METHOD, WITH LABEL SWTICHINGS, BUT ALSO
*WITH ADDED STEPS WHERE PARAMETERS OF DEAD CLUSTERS ARE UPDATED AFTER EACH UPDATE OF THE ALLOCATION
*VARIABLES.  WHEN A CURRENTLY EMPTY CLUSTER IS VISITED z(1,) AND z(2,)  ARE DRAWN FROM THE PRIOR. 
*WE ALSO GIVE THE POSSIBILITY TO INTRODUCE AUXILIARY STATES WHICH ARE CONSIDERED IN THE
*GIBBS SAMPLING; m_aug IS SPECIFIED AND max_k IS DEFINED AS max{k}+m_aug 
*
*WE DO NOT DO PARAMETER UPDATES
*.....parameters.................................................
      integer n,iter
      integer size,thin,burnin,density_samples
*metropolis parameters
      double precision scale_a
*size is the largest possible dimension of the variable-sized vectors
*thin is the thinning of the markov chain
      parameter(n=1000,iter=2000000,size=10000,thin=1,burnin=0)
      parameter(density_samples=5)
*
      parameter(scale_a=0.08d0)
*.....local scalars..............................................
*VARIABLES
*NOTATION: y        is the data, y(i), i=1,n
*          k        is the classification variables, k(i), i=1,n
*          z        is the set of position of random masses, z(j),j=1,2,...
*                   IN THIS PROGRAM z(j) IS TWO-DIMENSIONAL, z(1,j) is the mean, and z(2,j) the variance
*          v        is the set of gamma variables used to built the p, v(j),j=1,2,...
*          p        is the set of random masses, p(j),j=1,2,...
*          pprop    is the set of proposal probabilities - UN-NORMALISED!!!!!  
*          m        is the set of m(j), j=1,2,..., where m(j)=#{k(i)=j}
*          sumy     sumy(j)=sum{y(i):k(i)=j},j=1,2,...
*          sumysq   sumysq(j)=sum{(y(i)-z(1,j))^2:k(i)=j},j=1,2,...
*          al       is the set of indicators of alive points, i.e all js s.t m(j)=/=0
*          cl       is the size of al, i.e the number of distinct clusters
*          m_aug    the number of auxiliary states to be considered in the Gibbs sampling
*          max_k    is  max{k}+m_aug, where max{k} isdefined in the paper as the largest index of components 
*                   used to construct the proposal distribution 
*                   NOTICE that it might change at the update of each k(i)
*          sump     is the sum of the probabilities p(j) up to the furthest visited point
*                   we need it when wish to search beyond the alive points
*          sumpprop is the sum of proposal probablities for each data point used
*                   to be compared against a uniform in order to propose a new
*                   allocation variable
*          index    is a vector of size n containing the re-ordered indices of the 
*                   original data. This is due to the permutation we do in the updating
*                   of the allocation variables
*          m_norm   is the proportion constant used in the proposal of the allocation variables, 
*                   M in the paper, it controls the probability generating a new cluster from
*                   the prior. See the paper for a choice which depends on the Drichlet
*                   process
*          dev      is the model deviance, see Green and Richardson, we use it as a summary 
*                   of all updated components
****
*          mu    is the mean of z(1,)s
*          sz    is the standard deviation of z(1,)s
*          gamma is the shape parameter of z(2,)s
*          beta  is the scale parameter of z(2,)s
*          alpha is the concentration parameter of the DP
*
****REMARKS***
*   1) the initial allocation of data into clusters is done manually further down.
*      Some care is needed
*   2) Priors for parameters are defined directly in their subroutines, where all
*      hyperparameters and the functional form of the log-prior have to be specified
*
      double precision y(n)
      integer k(n),index(n),max_k,m_aug
      double precision z(2,size),v(size),p(size),sumy(size),sumysq(size)
      integer m(size),al(size),cl
      double precision sump,sumpprop
      double precision m_norm
      double precision dev
****
      double precision mu,gamma,beta,alpha,sz
****  
*.....dummy arguments..............................................
      double precision gam1,gam2,dumgamma(1)
      integer xmax,xmin
      double precision dat_range,prop_m_norm,datmin
      integer ifail,count,count_den,temp_m
      integer al_ind(n)
      double precision temp_z,dir
*temp_z,al_ind,dir  are variables for the label switching step
      double precision mean,var,summ,const,sumpinc,den
      double precision u
      integer propind,temp_k
      double precision log_prob,log_den,pi,prop_const
*
      double precision devi,maxp
      double precision grid(280),density(280,density_samples)
      integer dens_est,step
      double precision priorprop
      double precision c_time
*
*
*.....external functions...........................................
      double precision g05ddf,g05dbf,g05caf,g05dff
      integer g05dyf
      external g05ddf,g05dbf,g05caf,g05dff,g05dyf
*.....external routines............................................
      external g05ccf,g05fff,g05cbf,g01alf,g05ehf,g05ejf
*.....common.......................................................
*.....commands..................................................... 
      open(10,file='bimod-5000.dat')
      open(11,file='f_res.dat')
      open(12,file='f_res-density.dat')
      open(13,file='f_x1.dat')  
      open(14,file='f_x2.dat')
      open(15,file='prior_proposal_proportion.dat')
      open(16,file='f_random_measure.dat')
      open(17,file='f_random_pars.dat')
*************
*************
*SET THE SEED
************* 
*      call g05cbf(1)
      call g05ccf
*
      write(*,*) 'This is the CONDITIONAL RETROSPECTIVE algorithm'
      write(*,*) 'A LABEL SWITCHING STEP IS ATTEMPTED BETWEEN NEIGHBOURI
     *NG CLUSTERS'
      write(*,*) 'MEANS OF DEAD CLUSTERS ARE UPDATED AFTER ALLOCATION UP
     *DATE'
      write(*,*) 'SPECIFY NUMBER OF AUXILIARY STATES, m_aug=0,1,...'
*      read(*,*) m_aug
      m_aug=0
      write(*,*) 'DATASET SIZE IS ',n 
      write(*,*)
      write(*,*) 'number of iterations', iter
      write(*,*) 'chain is thinned every ', thin
      write(*,*) 'and the burnin is ', burnin
      write(*,*) 'therefore only values after ',burnin,'are saved'
      write(*,*) 'click to proceed'
*      read(*,*)
      write(*,*) 'Do you want to compute density estimates? (0/1)'
*      read(*,*) dens_est
      dens_est=0
      step=(iter-burnin)/density_samples
      if (dens_est.eq.1) then
         write(*,*) 'density will be computed'
         write(*,*) density_samples,' samples'
*find the grid
         grid(1)=-3.0d0
         do 1, ij=2,40
            grid(ij)=grid(ij-1)+0.05d0
 1       continue
         do 2, ij=1,200
            grid(40+ij)=grid(ij+40-1)+0.01d0
 2       continue
         do 3,ij=1,40
            grid(240+ij)=grid(240+ij-1)+0.05d0
 3       continue
      endif
**************
*READ THE DATA 
**************
      do 10, i=1,n
         read(10,*) y(i)
*index set
         index(i)=i
 10   continue
*
**********************
*INITIALIZE PARAMETERS
**********************
      alpha=1.00d0
*USE THE DATA-DEPENDENT PRIORS SUGGESTED BY GREEN & RICHARDSON, p.370
      call range_of_data(n,y,dat_range,datmin)
      mu=datmin+dat_range/2
      sz=dat_range
      gamma=2.0d0
*      gamma=0.10d0
      beta=0.02*(dat_range**2)
*      beta=0.01d0
*
*                   
      write(*,*) 'hyperparameter values:'
      write(*,*) 'alpha: ',alpha
      write(*,*) 'mu: ',mu
      write(*,*) 'sz: ',sz
      write(*,*) 'gamma: ',gamma
      write(*,*) 'beta: ',beta  
      write(*,*) 'click to proceed'
*      read(*,*)
      pi=dacos(-1.0d0)
*
********************************
*PUT 0s IN ALL COUNTERS AND SUMS
********************************
      do 11, j=1,size
         m(j)=0
         sumy(j)=0.0d0
         sumysq(j)=0.0d0
 11   continue
      count=0
      count_den=0
*
*************************************************************************
*INITIAL NUMBER OF CLUSTERS                                             *
*************************************************************************
*  
      cl=1
      al(1)=1
*      do 23,j=1,cl                                                      
*         al(j)=j                                                        
* 23   continue                                                   
************************************************                             
*DECIDE INITIAL ALLOCATION OF DATA INTO CLUSTERS                                            
************************************************
*important cluster summaries:
      do 21,i=1,n                                                       
         k(i)=1                                                         
         m(k(i))=m(k(i))+1                                              
         sumy(k(i))=sumy(k(i))+y(i)                                                                             
         sumysq(k(i))=sumysq(k(i))+(y(i)-z(1,k(i)))**2 
 21   continue                                                     
**********************************                                      
*GIVE BIRTH TO THE REQUIRED POINTS                                      
**********************************                                                       
      call integer_extremes(n,k,xmax,xmin)
      max_k=xmax+m_aug                                           
*simulate first v(1),p(1),z(1,1),z(2,1)                                          
      ifail=0                                                           
      call g05fff(1.0d0,1.0d0,1,dumgamma,ifail)                            
      gam1=dumgamma(1)                                                     
      call g05fff(alpha,1.0d0,1,dumgamma,ifail)                      
      gam2=dumgamma(1)                                                     
      v(1)=gam1/(gam1+gam2)                                             
      p(1)=v(1)  
      z(1,1)=g05ddf(mu,sz)                                          
      ifail=0                                                           
      call g05fff(gamma,1.0d0/beta,1,dumgamma,ifail)   
      z(2,1)=1.0d0/dumgamma(1)
      sump=p(1)                                 
      do 20,j=2,max_k                                                    
*simulate p(j) from the prior                                           
         ifail=0                                                        
         call g05fff(1.0d0,1.0d0,1,dumgamma,ifail)                         
         gam1=dumgamma(1)                                                  
         call g05fff(alpha,1.0d0,1,dumgamma,ifail)                   
         gam2=dumgamma(1)                                                  
         v(j)=gam1/(gam1+gam2)                                                
         p(j)=v(j)*(1.0d0-sump)
         sump=sump+p(j)
*simulate z(1,j) and z(2,j) from the prior                                           
         z(1,j)=g05ddf(mu,sz)                                                            
         call g05fff(gamma,1.0d0/beta,1,dumgamma,ifail)   
         z(2,j)=1.0d0/dumgamma(1) 
 20   continue    
*                               
      write(*,*) 'initial number of clusters', cl
      write(*,*) 'cluster parameters', m(k(1)),sumy(k(1)),sumysq(k(1))                                        
***********************
*CALCULATE THE DEVIANCE
***********************
      dev=0.0d0
      do 24, i=1,n
         devi=0.0d0
         do 25, j=1,cl
            devi=devi+((m(al(j))*1.0d0)/(n*sqrt(2*pi*z(2,al(j)))))*exp(-
     *      (0.5d0/z(2,al(j)))*(y(i)-z(1,al(j)))**2)
 25      continue
        if (devi.le.0.0d0) write(*,*) 'problem in deviance'
         dev=dev+log(devi)
 24   continue
      dev=-2*dev   
      write(11,*) cl,dev   
      write(*,*) 'initial deviance', dev                
      write(*,*) 'start Gibbs samping - click to continue'  
*     
*      read(*,*) 
*************************************************************************
***********************************
*START THE GIBBS SAMPLING ALGORITHM
***********************************
*
      do 30,l=1,iter
*TO MEASURE CPU TIME IN SEC
         if (l.eq.50000) then
            call CPU_Time(c_time)
            write(*,*) c_time
         endif
*   
*calculate max_k                                                   
         call integer_extremes(n,k,xmax,xmin)
         max_k=xmax+m_aug 
*************************************************
*STEP 1: UPDATE THE RANDOM MEASURE 
*************************************************
*update data summaries 
         do 42,i=1,n                                  
            sumysq(k(i))=0.0d0
            sumy(k(i))=0.0d0
 42      continue
         do 41,i=1,n        
            sumy(k(i))=sumy(k(i))+y(i)
            sumysq(k(i))=sumysq(k(i))+(y(i)-z(1,k(i)))**2 
 41      continue
         summ=0.0d0
         sump=0.0d0
         maxp=0.0d0
*update the zs,ps until max_k
*DO NOT SIMULATE z VALUES FOR DEAD CLUSTERS SINCE THIS ARE SIMULATED AGAIN
*IN THE UPDATING OF THE ALLOCATION VARIABLES
         do 40,j=1,max_k
*update zjs - if the cluster is alive
            if (m(j).eq.0) goto 43
*z(2,j) first, conditionally on z(1,j)                                                         
            call g05fff(gamma+m(j)*1.0d0/2.0d0,1.0d0/(beta+sumysq(j)/2.0
     *      d0),1,dumgamma,ifail)   
            z(2,j)=1.0d0/dumgamma(1)
*
*z(1,j) conditionally on z(2,j)
            var=1.0d0/(m(j)/z(2,j)+1.0d0/(sz**2))
            mean=var*(mu/(sz**2)+sumy(j)/z(2,j))
            z(1,j)=g05ddf(mean,sqrt(var))         
*update pjs
 43         summ=summ+m(j)
            ifail=0
            call g05fff((m(j)+1.0d0),1.0d0,1,dumgamma,ifail)
            gam1=dumgamma(1)
            call g05fff((alpha+n-summ),1.0d0,1,dumgamma,ifail)
            gam2=dumgamma(1)   
            v(j)=gam1/(gam1+gam2)
            p(j)=(1.0d0-sump)*v(j)
            sump=sump+p(j)
            if (sump.gt.1.0d0) then
               write(*,*) 'error after 50'
               read(*,*)
            endif
            if (maxp.le.p(j)) maxp=p(j)
 40      continue
*************************************************
*STEP 2.1: UPDATE THE CLASSIFICATION VARIABLES
*************************************************
*do a permutation of indices
*
         ifail=0
         call g05ehf(index,n,ifail)
*
* Use the quasi-independence sampler
         do 60,i=1,n    
            call integer_extremes(n,k,xmax,xmin)
            max_k=xmax+m_aug        
            ip=index(i)
*calculate first the scaling factor M  and constant
*as suggested in the paper we take M=max(f(Yi|Z_j),j<kmax)
            m_norm=0.0d0
            const=0.0d0
            sumpinc=0.0d0
            sumpprop=0.0d0
            sump=0.0d0
*re-calculate also the sum of prior probabilities, the fact that we change random measure according to max_k
*can create inconsistencies in the sump counter, which might have included from previous iteration a currently
*incorrect number (since some clusters have in the meantime become dead)
            do 65,j=1,max_k
               sumpinc=sumpinc+p(j)
               sump=sump+p(j)
**********************************************************************
*FOR ANY CURRENTLY DEAD CLUSTER SIMULATE THE PARAMETERS FROM THE PRIOR
**********************************************************************
               if (m(j).eq.0) then
                  z(1,j)=g05ddf(mu,sz)                                                           
                  call g05fff(gamma,1.0d0/beta,1,dumgamma,ifail)   
                  z(2,j)=1.0d0/dumgamma(1)
               endif
*the density of data
               den=(1.0d0/sqrt(2*pi*z(2,j)))*exp(-(0.5d0/z(2,j))*((y(ip)
     *         -z(1,j))**2))
               if (den.ge.m_norm) m_norm=den
               const=const+p(j)*den
 65         continue
            if (m_norm.eq.0.0d0) then
               write(*,*) 'problem in the proposal, m_norm'
               read(*,*)
            endif
            const=const+(1.0d0-sumpinc)*m_norm
*the probability of proposing from the prior is (1-sumpinc)*m_norm/max_k 
*This implies a probability of proposing from the prior equal to
            priorprop=(1.0d0-sumpinc)*m_norm/const
            if (priorprop.lt.0.0d0) then
               write(*,*) 'problem in mnorm'
               write(*,*) max_k,sumpinc,sump,ip
               write(*,*) (p(ij),ij=1,max_k)
               read(*,*)
            endif
            u=g05caf(u)
*start searching for the required cluster - SEARCH FIRST UP TO THE KNOWN ALIVE & GHOST POINTS
            do 62,j=1,max_k
*the proposal probability is p(j)*den/const, so:
*the density of data
               den=(1.0d0/sqrt(2*pi*z(2,j)))*exp(-(0.5d0/z(2,j))*((y(ip)
     *            -z(1,j))**2))
               sumpprop=sumpprop+p(j)*den/const
               if (sumpprop.gt.u) goto 63
 62         continue
*
 63         if (sumpprop.gt.u) then
               propind=j
            else
*continue searching in unvisited clusters (of current size 0)
*reset j, which is already increased by 1 in the last loop
               j=max_k
 64            j=j+1
               if (j.ge.size) then
                  write(*,*) 'exceeded limits for vector lengths'
                  write(*,*) '****POSSIBLY PROBLEM****'
                  write(*,*) size,sumpprop
                  write(*,*) 'iteration',l
                  write(*,*) 'max_k',max_k
                  write(*,*) 'data',i
                  read(*,*)
               endif
*simulate parameters for the newly-visited cluster: FROM THE PRIOR (since they correspond to the
*previous step of the algorithm
               ifail=0
               call g05fff(1.0d0,1.0d0,1,dumgamma,ifail)
               gam1=dumgamma(1)
               call g05fff(alpha,1.0d0,1,dumgamma,ifail)
               gam2=dumgamma(1)
               v(j)=gam1/(gam1+dumgamma(1))
               p(j)=(1.0d0-sump)*v(j)
               if (maxp.le.p(j)) maxp=p(j)
               sump=sump+p(j)
               sumpprop=sumpprop+m_norm*p(j)/const
               if (sump.gt.1.0d0) then
                  write(*,*) 'error after 64'
                  read(*,*)
               endif 
 2000          if (sumpprop.le.0.0d0) write(*,*) 'problem in 2000'
               if (log(sumpprop).gt.10.d-15) then
                  write(*,*) 'error after 64 - pprop',j,log(sumpprop)
                  write(*,*) u,v(j),sumpprop-m_norm*p(j)/const
                  write(*,*) sump-p(j)
                  write(*,*) nast,sump,sumpprop
                  write(*,*) (p(ij),ij=1,j)
                  write(*,*) m_norm,const,ip,l
                  read(*,*)
               endif
*simulate z(1,j) and z(2,j) from the prior                                           
               z(1,j)=g05ddf(mu,sz)                                                            
               call g05fff(gamma,1.0d0/beta,1,dumgamma,ifail)   
               z(2,j)=1.0d0/dumgamma(1)
               if (sumpprop.lt.u) goto 64
*
*at this point it has exited and propind=j
               propind=j 
            endif
*At this point the proposed value propind has been generated. The acceptance probability
*has to be computed and all relevant parameters to be updated.
*The acceptance probability will depend on following scenarios:
            if (propind.le.max_k) then
*
*find the maximum value of the vector k, where ip-element is replaced by propind
*if this changes by the proposed move detailed balance requires an appropriate A/R
*
***************TEMPORARILY!!********
               temp_k=k(ip)
               k(ip)=propind
               call integer_extremes(n,k,xmax,xmin)
               xmax=xmax+m_aug
               k(ip)=temp_k
************************************
*SCENARIO 1
               if (xmax.eq.max_k) then
*accept with probability 1
                  goto 70
               endif
*the acceptance probability will depend on whether k(ip)>xmax or not
*SCENARIO 2
               if ((xmax.gt.max_k).or.((xmax.lt.max_k).and.(k(ip).le.xma
     *         x))) then
*moves k(ip)-->propind and propind-->k(ip) are done with Gibbs, but the normalising constant differs
*this move changes the proposal distribution and m_norm
*NOTICE: this move might involve new clusters if xmax>max_k
*the probabilities and parameters of the new clusters are so far unknown
******* these are the clusters from max_k+1 till xmax, we first generate these values
*notice that this will not be excuted if xmax<max_k
                  do 66,ij=max_k+1,xmax
                     call g05fff(1.0d0,1.0d0,1,dumgamma,ifail)
                     gam1=dumgamma(1)
                     call g05fff(alpha,1.0d0,1,dumgamma,ifail)
                     gam2=dumgamma(1)
                     v(ij)=gam1/(gam1+dumgamma(1))
                     p(ij)=(1.0d0-sump)*v(ij)
                     if (maxp.le.p(ij)) maxp=p(ij)
                     sump=sump+p(ij)
*simulate z(1,j) and z(2,j) from the prior                                           
                     z(1,ij)=g05ddf(mu,sz)                                                            
                     call g05fff(gamma,1.0d0/beta,1,dumgamma,ifail)   
                     z(2,ij)=1.0d0/dumgamma(1)
 66               continue
                  prop_m_norm=0.0d0
                  prop_const=0.0d0
                  sumpinc=0.0d0
                  do 81,j=1,xmax
                     sumpinc=sumpinc+p(j)
*the density of data
                     den=(1.0d0/sqrt(2*pi*z(2,j)))*exp(-(0.5d0/z(2,j))*(
     *               (y(ip)-z(1,j))**2))
                     if (den.ge.prop_m_norm) prop_m_norm=den
                     prop_const=prop_const+p(j)*den
 81               continue
                  if (prop_m_norm.eq.0.0d0) then
                     write(*,*) 'problem in scenario 2, m_norm'
                     read(*,*)
                  endif
                  prop_const=prop_const+(1.0d0-sumpinc)*prop_m_norm
                  log_prob=log(const)-log(prop_const)
                  u=g05caf(u)
                  if (log(u).le.log_prob) then
*accept the move
                     goto 70
                  else
*reject the move- no changes needed
                     goto 71
                  endif
               endif
               if (xmax.lt.k(ip)) then
*move k(ip)-->propind is with Gibbs but move from propind-->k(ip) is from the prior
                  prop_m_norm=0.0d0
                  prop_const=0.0d0
                  sumpinc=0.0d0
                  do 82,j=1,xmax
*the density of data
                     den=(1.0d0/sqrt(2*pi*z(2,j)))*exp(-(0.5d0/z(2,j))*(
     *            (y(ip)-z(1,j))**2))
                     if (den.ge.prop_m_norm) prop_m_norm=den
                     prop_const=prop_const+p(j)*den
                     sumpinc=sumpinc+p(j)
 82               continue
                  if (prop_m_norm.eq.0.0d0) then
                     write(*,*) 'problem in scenario 2, m_norm'
                     read(*,*)
                  endif
                  prop_const=prop_const+(1.0d0-sumpinc)*prop_m_norm
*the log-density of data evaluated at the current data point: 
calucalte log directly is mORE stable, 0's have been noticed
*when density is evaluated - log works ok
                  log_den=-0.5d0*log(2*pi*z(2,k(ip)))-(0.5d0/z(2,k(ip)))
     *          *((y(ip)-z(1,k(ip)))**2)
                  log_prob=log(const)+log(prop_m_norm)-log(prop_const)-l
     *            og_den
                  u=g05caf(u)
                  if (log(u).le.log_prob) then
*accept the move
                     goto 70
                  else
*reject the move - no changes needed
                     goto 71
                  endif
               endif
*
            else
***************TEMPORARILY!!********
               temp_k=k(ip)
               k(ip)=propind
               call integer_extremes(n,k,xmax,xmin)
               xmax=xmax+m_aug
               k(ip)=temp_k
************************************ 
*SCENARIO 3
*the normalising contant is different at the proposed configuration of the
*classification variables. Needs to be computed*also m_norm is different
*NOTICE: this move will involve new clusters whose probabilities and parameters are so far unknown
******* these are the clusters from propind+1 till xmax, we first generate these values
                  do 67,ij=propind+1,xmax
                     call g05fff(1.0d0,1.0d0,1,dumgamma,ifail)
                     gam1=dumgamma(1)
                     call g05fff(alpha,1.0d0,1,dumgamma,ifail)
                     gam2=dumgamma(1)
                     v(ij)=gam1/(gam1+dumgamma(1))
                     p(ij)=(1.0d0-sump)*v(ij)
                     if (maxp.le.p(ij)) maxp=p(ij)
                     sump=sump+p(ij)
*simulate z(1,j) and z(2,j) from the prior                                           
                     z(1,ij)=g05ddf(mu,sz)                                                            
                     call g05fff(gamma,1.0d0/beta,1,dumgamma,ifail)   
                     z(2,ij)=1.0d0/dumgamma(1)
 67               continue
               prop_m_norm=0.0d0
               prop_const=0.0d0
               sumpinc=0.0d0
               do 91,j=1,xmax
*the density of data
                  den=(1.0d0/sqrt(2*pi*z(2,j)))*exp(-(0.5d0/z(2,j))*((y(
     *            ip)-z(1,j))**2))
                  if (den.ge.prop_m_norm) prop_m_norm=den
                  prop_const=prop_const+p(j)*den
                  sumpinc=sumpinc+p(j)
 91            continue
               if (prop_m_norm.eq.0.0d0) then
                  write(*,*) 'problem in scenario 2, m_norm'
                  read(*,*)
               endif
               prop_const=prop_const+(1.0d0-sumpinc)*prop_m_norm
*the log-density of data: calucalte log directly is mORE stable, 0's have been noticed
*when density is evaluated - log works ok
*               den=(1.0d0/sqrt(2*pi*z(2,propind)))*exp(-(0.5d0/z(2,propi
*     *         nd))*((y(ip)-z(1,propind))**2))
               log_den=-0.5d0*log(2*pi*z(2,propind))-(0.5d0/z(2,propind)
     *          )*((y(ip)-z(1,propind))**2)
               log_prob=log(const)-log(m_norm)-log(prop_const)+log_den
               u=g05caf(u)
               if (log(u).le.log_prob) then
*accept the move
                  goto 70
               else
*reject the move - no changes needed
                  goto 71
               endif
            endif
*******************************
*move ip ---> propind is accepted
 70         continue
*update parameters of previous cluster
            m(k(ip))=m(k(ip))-1
            if ((m(k(ip)).eq.0).and.(propind.ne.k(ip))) then
*we need to update al and cl
               call remove_cluster(al,cl,k(ip))
*update sumsq to ensure that it is 0 - THIS IS DONE SINCE AT THE UPDATE OF RANDOM MEASURE WE UPDATE CLUSTERS 1:nstar SINCE 
*THE FORMULAS ARE GENERAL AND DO NOT DISTINGUISH BETWEEN ALIVE AND DEAD. THE INITIALISATION FOR sumysq IN THE UPDATE OF THE
*RANDOM MEASURE HAPPENS ONLY FOR ALIVE CLUSTERS
            endif
*update the summaries of the new cluster 
            m(propind)=m(propind)+1
            if ((m(propind).eq.1).and.(propind.ne.k(ip))) then
*a new cluster is born - it needs to be added to alive clusters
*accept the move - update al,cl
               call add_cluster(al,cl,propind)
            endif 
            k(ip)=propind
            goto 71
*******************************
 71         continue
*continue updating the rest of allocation variables
 60      continue                        
*************************************************
*STEP 2.2: DO A LABEL SWITCHING STEP 
*************************************************
*here we implement the following scheme: choose two ALIVE clusters at 
*random and propose to exchange their indices
         if (cl.le.1) goto 94
         ifail=0
         call g05ejf(al,cl,al_ind,2,ifail)
*propose an exchange 
         u=g05caf(u)
         if (log(u).le.(m(al_ind(2))-m(al_ind(1)))*(log(p(al_ind(1)))-lo
     *      g(p(al_ind(2))))) then
            temp_z=z(1,al_ind(1))
            z(1,al_ind(1))=z(1,al_ind(2))
            z(1,al_ind(2))=temp_z  
            temp_z=z(2,al_ind(1))
            z(2,al_ind(1))=z(2,al_ind(2))
            z(2,al_ind(2))=temp_z
            temp_z=sumy(al_ind(1))
            sumy(al_ind(1))=sumy(al_ind(2))
            sumy(al_ind(2))=temp_z
*
            temp_m=m(al_ind(1))
            m(al_ind(1))=m(al_ind(2))
            m(al_ind(2))=temp_m
*change the allocation indices:
            do 72, i=1,n
               if (k(i).eq.al_ind(1)) then
                  k(i)=al_ind(2)
               elseif (k(i).eq.al_ind(2)) then
                  k(i)=al_ind(1)
               endif
 72         continue      
         endif
***************************
*SAVE INTERESTING SUMMARIES
***************************
 94      continue
         if ((((l-burnin)/thin-count).ne.0).and.(l.ge.burnin)) then 
*write iteration number on screen         
*            write(*,*) l
            count=count+1
*
*CALCULATE THE DEVIANCE
*
           dev=0.0d0
           do 95, i=1,n
              devi=0.0d0
              do 96, j=1,cl
                 devi=devi+(m(al(j))*1.0d0)/(n*sqrt(2*pi*z(2,al(j))))*ex
     *           p(-(0.5d0/z(2,al(j)))*(y(i)-z(1,al(j)))**2)
 96           continue
              dev=dev+log(devi)
 95        continue
           dev=-2*dev
           write(11,*) cl,dev
*the mean of y(3): 
           write(13,*) z(1,k(3))         
           write(14,*) z(1,k(2))
           write(15,*) priorprop
           if ((l.ge.50000).and.(l.le.80000)) then
              write(16,99) (p(j),j=1,4)
              write(17,98) (z(1,j),j=1,4)
           endif
        endif
        if ((((l-burnin)/step-count_den).ne.0).and.(l.ge.burnin).and.(de
     *  ns_est.eq.1)) then
           count_den=count_den+1
*
*CALCULATE THE DENSITY ESTIMATE
           do 1000, ij=1,280             
              devi=0.0d0
              do 1010,j=1,cl
                 devi=devi+(m(al(j))*1.0d0)/(n*sqrt(2*pi*z(2,al(j))))*ex
     *           p(-(0.5d0/z(2,al(j)))*(grid(ij)-z(1,al(j)))**2)
 1010         continue
              density(ij,count_den)=devi
 1000      continue
        endif
 30   continue
      if (dens_est.eq.1) then
         do 1100,ij=1,280
            write(12,99) grid(ij),(density(ij,j),j=1,5)
 1100    continue
      endif
 98   format(4(1x,f10.4))
 99   format(4(1X,f8.6))
      end
********************************************************************************
******SUBROUTINE****************************************************************
      subroutine remove_cluster(al,cl,ind)
*this removes clustr ind from the vector of existing clusters al, reduces the
*number of clusters cl by one, and reorders the elements of al to be in 
*increasing order
*.....local scalars................................................
      integer al(*),cl,ind,j
*.....commands.....................................................
      do 100,j=1,cl-1
         if (al(j).lt.ind) then
            goto 100
         else
            al(j)=al(j+1)
         endif
 100  continue
      cl=cl-1
      return
      end 
********************************************************************************
******SUBROUTINE****************************************************************
      subroutine add_cluster(al,cl,ind)
*this adds cluster ind to the vector of existing clusters al, increases the
*number of clusters cl by one, and reorders the elements of al to be in 
*increasing order
*.....local scalars................................................
      integer al(*),cl,ind,j,k
*.....commands.....................................................
      do 200,j=1,cl
         if (al(j).gt.ind) goto 210
 200  continue
 210  continue
      do 220,k=cl+1,j+1,-1
         al(k)=al(k-1)
 220  continue
      al(j)=ind
      cl=cl+1
      return
      end 
********************************************************************************
******SUBROUTINE****************************************************************
      subroutine integer_extremes(n,x,xmax,xmin)
*it finds the maximum and minimum elements of an integer vector
*.....local scalars................................................
      integer n,x(n),xmax,xmin
*.....commands.....................................................
      xmax=x(1)
      xmin=x(1)
      do 300,j=1,n
         if (x(j).le.xmin) xmin=x(j)
         if (x(j).ge.xmax) xmax=x(j)
 300  continue
      return
      end 
********************************************************************************
******SUBROUTINE****************************************************************
      subroutine range_of_data(n,x,dat_range,datmin)
*it finds the range of the data x
*.....local scalars................................................
      integer n
      double precision x(n),xmax,xmin,dat_range,datmin
*.....commands.....................................................
      xmax=x(1)
      xmin=x(1)
      do 400,j=1,n
         if (x(j).le.xmin) xmin=x(j)
         if (x(j).ge.xmax) xmax=x(j)
 400  continue
      dat_range=xmax-xmin
      datmin=xmin
      return
      end


