PROGRAM spect_mix_logit_rj_sim2
  USE lin_eig_self_int
  USE show_int
  USE wrrrn_int
  USE rnset_int
  IMPLICIT none
  INTERFACE
   FUNCTION log_spect_ar2(nfreq,comp)
    IMPLICIT none
    INTEGER, INTENT(IN) :: nfreq, comp
    DOUBLE PRECISION, DIMENSION(nfreq) :: log_spect_ar2
   END FUNCTION
  END INTERFACE 
  INTEGER, PARAMETER :: nobs=64, nbasis=10
  INTEGER, PARAMETER :: nchunks=16
  INTEGER, PARAMETER :: nsim=100
  INTEGER, PARAMETER :: nfreq=nobs/2 + 1
  INTEGER, PARAMETER :: ny=nchunks*nfreq
  INTEGER, PARAMETER :: nbeta=nbasis+1
  INTEGER, PARAMETER :: nloop=100000,nwarmup=20000
  INTEGER :: i, j, k,iseed
  DOUBLE PRECISION, allocatable :: tp(:),omega(:,:),x(:,:) 
  DOUBLE PRECISION, allocatable :: xx(:,:),y(:,:),yy(:,:),xdelta(:,:)
  DOUBLE PRECISION, allocatable :: d(:), v_s(:,:)
  DOUBLE PRECISION :: sigalpha,sigdelta
  DOUBLE PRECISION, DIMENSION(nfreq,nchunks) :: fit
  DOUBLE PRECISION, DIMENSION(nfreq,nchunks) :: log_spect_true
  INTEGER, DIMENSION(nchunks,nsim) :: true_ind
  DOUBLE PRECISION :: mse
  DOUBLE PRECISION, DIMENSION(nchunks-1) :: clevel
  INTEGER, DIMENSION(nchunks-1) :: iclson,icrson

  
  open(10, file = 'log_per_ar10_100sim')
  open(12, file = 'iter_ar10_100sim22')
  open(15, file = 'fitted_ar10_100sim22')
  open(20, file = 'lcl_ar10_100sim22')
  open(22, file = 'ucl_ar10_100sim22')
  iseed=24378
  !iseed=54674
  call rnset(iseed)
  ALLOCATE(tp(nfreq))
  ALLOCATE(omega(nfreq,nfreq))
  DO i = 1, nfreq
     !tp(i)=(dble(i)-1.0d0)/(2.0d0*(dble(nfreq)-1.0d0))
     tp(i)=(2.0d0*dble(i)-1.0d0)/(2.0d0*dble(nfreq))
  END DO 
  DO i=1,nfreq
     DO j=i,nfreq
        omega(i,j)=tp(i)**2*(tp(j)-tp(i)/3.0d0)/2.0d0
     END DO 
  END DO 
  DO i=2,nfreq
     DO j=1,i-1
        omega(i,j)=omega(j, i)
     END DO 
  END DO 
  !call wrrrn('omega',omega,nra=nfreq)
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  ALLOCATE(d(nfreq),v_s(nfreq,nfreq))
  !Find the eigenvalue decomposition
  call lin_eig_self(omega,d,v=v_s)
  ALLOCATE(x(nfreq,nbasis))
  DO i=1,nfreq
     DO j=1,nbasis
        x(i,j)=(-1)**(j+1)*v_s(i,j)*sqrt(d(j))
     END DO
  END DO
  !Note: multiply by (-1)**(j+1) to match matlab's eigenvectors
  !do i=1,nfreq
  !   print*, x(i,3)
  !end do
  !construct the matrix xx in which the first column is 1's
  !and the rest of the columns are the matrix x. This combines
  !the Z and X matrices. Z is just a column of 1's since alpha_1
  !is constrained to be zero
  ALLOCATE(xx(nfreq,nbeta))
  xx(:,1)=1.0d0
  DO j=2,nbeta
     xx(:,j)=x(:,nbasis-j+2)
  END DO
  !DO j=1,nbeta
  !   DO i=1,nfreq
  !     read(8,*) xx(i,j)
  !   END DO
  !END DO
  DEALLOCATE(x)
  ALLOCATE(xdelta(nchunks,2))
  xdelta(:,1)=1.0d0
  xdelta(:,2)=dble((/(k,k=1,nchunks)/))/dble(nchunks)
  sigalpha=100.0d0
  sigdelta=4.0d0
  ALLOCATE(yy(nchunks*nfreq,nsim),y(nfreq,nchunks))
  DO j=1,nsim
     DO i=1,nchunks*nfreq
       read(10,*) yy(i,j)
     END DO
  END DO
!  DO i=1,nchunks
!     IF(i.LE.8) true_ind(i,:)=1
!     IF(i.GE.9 .AND. i.LE.12) true_ind(i,:)=2
!     IF(i.GE.13 .AND. i.LE.16) true_ind(i,:)=3
!  END DO
  DO i=1,nchunks
     IF(i.LE.3) true_ind(i,:)=1
     IF(i.GE.4 .AND. i.LE.16) true_ind(i,:)=2
  END DO
  DO k=98,98
     print*,k
     y=reshape(yy(:,k),(/nfreq,nchunks/))
     call revjump(iseed,nloop,nwarmup,nfreq,nchunks,nbeta,&
                  nbasis,sigalpha,sigdelta,xx,xdelta,y,fit,&
                  clevel,iclson,icrson)
     !DO j=1,nchunks
     !   log_spect_true(:,j)=log_spect_ar2(nfreq,true_ind(j,k))
     !END DO
     !mse=sum((fit-log_spect_true)**2)/dble(nchunks*nfreq)
     DO j=1,nchunks
        DO i=1,nfreq
           write(15,*) fit(i,j)
        END DO
     END DO
     !write(20,*) mse
  END DO
END PROGRAM spect_mix_logit_rj_sim2
!
!
SUBROUTINE revjump(iseed,nloop,nwarmup,nfreq,nchunks,nbeta,&
                 nbasis,sigalpha,sigdelta,xx,xdelta,y,fit,clevel,iclson,icrson)
  !USE show_int
  USE wrrrn_int
  USE rnunf_int
  USE linear_operators
  USE lin_eig_self_int
  USE eqtil_int
  USE maxb_m, ONLY: maxb
  IMPLICIT none
  INTERFACE
    FUNCTION findloc(m,n,ind,ii)
     INTEGER, INTENT(IN) :: m,n
     INTEGER, DIMENSION(n), INTENT(IN) :: ind
     INTEGER, DIMENSION(m) :: findloc
     INTEGER :: ii,j,k
    END FUNCTION
    FUNCTION mvrnorm(iseed,n,mean,var)
      USE lin_eig_self_int
      USE linear_operators
      IMPLICIT none
      INTEGER, INTENT(IN) :: n,iseed
      DOUBLE PRECISION, DIMENSION(n), INTENT(IN) :: mean
      DOUBLE PRECISION, DIMENSION(n,n), INTENT(IN) :: var
      DOUBLE PRECISION, DIMENSION(n) :: mvrnorm
    END FUNCTION
    FUNCTION mvrnorm1(iseed,n,mean,var)
      USE linear_operators
      USE show_int
      IMPLICIT none
      INTEGER, INTENT(IN) :: n,iseed
      DOUBLE PRECISION, DIMENSION(n), INTENT(IN) :: mean
      DOUBLE PRECISION, DIMENSION(n,n), INTENT(IN) :: var
      DOUBLE PRECISION, DIMENSION(n) :: mvrnorm1
    END FUNCTION
    FUNCTION log_like(nfreq,nchunks,nbeta,nexp,xx,y,beta,prior_ind)
      IMPLICIT none
      INTEGER, INTENT(IN) :: nfreq,nchunks,nbeta,nexp
      DOUBLE PRECISION, DIMENSION(nfreq,nbeta), INTENT(IN) :: xx
      DOUBLE PRECISION, DIMENSION(nfreq,nchunks), INTENT(IN) :: y
      DOUBLE PRECISION, DIMENSION(nexp,nchunks), INTENT(IN) :: prior_ind
      DOUBLE PRECISION, DIMENSION(nexp,nbeta), INTENT(IN) :: beta
      DOUBLE PRECISION :: log_like
    END FUNCTION
    FUNCTION cumsum(n,x)
      IMPLICIT none
      INTEGER, INTENT(IN) :: n
      DOUBLE PRECISION, DIMENSION(n), INTENT(IN) :: x
      DOUBLE PRECISION, DIMENSION(n) :: cumsum
    END FUNCTION
  END INTERFACE
  INTEGER, INTENT(IN) :: nfreq,nchunks,nbeta,nbasis,nloop,&
                         nwarmup,iseed
  DOUBLE PRECISION, INTENT(IN) :: sigalpha,sigdelta
  DOUBLE PRECISION, DIMENSION(nfreq,nbeta), INTENT(IN) :: xx
  DOUBLE PRECISION, DIMENSION(nchunks,2), INTENT(IN) :: xdelta
  DOUBLE PRECISION, DIMENSION(nfreq,nchunks), INTENT(IN) :: y
  DOUBLE PRECISION, DIMENSION(nfreq,nchunks), INTENT(OUT) :: fit
  INTEGER, PARAMETER :: nquan=2, nexp_max=10
  INTEGER :: nexperts
  DOUBLE PRECISION, ALLOCATABLE :: beta(:,:),delta(:,:),tau(:),prior_ind(:,:),&
                       prior_ind_prop(:,:),delta_prop(:,:),tau_prop(:),&
                       beta_prop(:,:)
  DOUBLE PRECISION, DIMENSION(2) :: mean_d
  DOUBLE PRECISION, DIMENSION(2,2) :: var_d
  DOUBLE PRECISION, DIMENSION(nbeta) :: mean_b
  DOUBLE PRECISION, DIMENSION(nbeta,nbeta) :: var_b
  DOUBLE PRECISION, ALLOCATABLE :: h(:,:)
  DOUBLE PRECISION :: u,drnunf
  DOUBLE PRECISION :: u1,u2
  DOUBLE PRECISION :: log_like_curr,log_like_new,met_rat_jump,epsilon_jump
  INTEGER :: ndelta,nproposed,counter
  INTEGER :: i,j,k,iter
  INTEGER, DIMENSION(nchunks) :: ind0
  INTEGER, DIMENSION(nchunks) :: ind,kk
  INTEGER, ALLOCATABLE :: kk_ind(:)
  INTEGER :: n_chunk_ind
  DOUBLE PRECISION, ALLOCATABLE :: hh(:)
  DOUBLE PRECISION, ALLOCATABLE :: post_ind(:,:)
  DOUBLE PRECISION, DIMENSION(nbeta) :: param
  DOUBLE PRECISION, DIMENSION(nfreq) :: fit_chunk
  DOUBLE PRECISION, DIMENSION(nfreq,nchunks) :: fit_ch
  DOUBLE PRECISION, DIMENSION(nchunks) :: sum_fit_ch
  DOUBLE PRECISION, DIMENSION(nfreq,nchunks) :: fhat
  DOUBLE PRECISION, DIMENSION(nfreq,nchunks) :: cum_fhat
  DOUBLE PRECISION, DIMENSION(nchunks-1) :: dist
  DOUBLE PRECISION, DIMENSION(nchunks-1) :: clevel
  INTEGER, DIMENSION(nchunks-1) :: iclson,icrson
  DOUBLE PRECISION, DIMENSION(nfreq,3) :: fit_exp,fit_exp_sum
  INTEGER :: count_fit, nmiss
  DOUBLE PRECISION :: qprop(nquan), quan(nquan), xlo(nquan), xhi(nquan)
  DOUBLE PRECISION :: lcl(nfreq,nchunks), ucl(nfreq,nchunks)
  DOUBLE PRECISION :: fit_it(nfreq,nloop-nwarmup,nchunks)
  !DOUBLE PRECISION :: prior_ind_hat(4,nchunks)
  INTEGER :: acc_death, acc_birth
  !
  nexperts=1
  ndelta=2*(nexperts-1)
  ALLOCATE(tau(nexperts),beta(nexperts,nbeta),delta(2,nexperts),&
           prior_ind(nexperts,nchunks),h(nexperts,nchunks))


  tau(1)=10000.0d0
  !Note that if nexperts<2 the loop is not executed
  DO j=2,nexperts
     tau(j)=tau(j-1)/10.0d0
  END DO
  delta=0.0d0
  mean_d=0.0d0
  mean_b=0.0d0
  var_d=sigdelta*eye(2)
  !call em(iseed,nfreq,nchunks,nbeta,ndelta,nexperts,nbasis,sigalpha,&
  !            sigdelta,xx,xdelta,y,ind0,beta,prior_ind)
  !ind=ind0
  DO j=1,nexperts
     call drnun(nchunks,h(j,:))
  END DO
  ind=maxloc(h,dim=1)
  DEALLOCATE(h)
  DO j=1,nexperts
     kk=0
     where(ind.eq.j) kk=1
     n_chunk_ind=sum(kk)
     ALLOCATE(kk_ind(n_chunk_ind),hh(n_chunk_ind))
     !find the chunks belonging to expert j
     kk_ind=findloc(n_chunk_ind,nchunks,ind,j)
     hh=1.0d0
     param=0.0d0
     call maxb(param,tau(j),hh,sigalpha,nfreq,n_chunk_ind,nbeta,&
               nexperts,nbasis,xx,y(:,kk_ind),beta(j,:))
     DEALLOCATE(kk_ind,hh)
  END DO
  fit=0.0d0
  !!!!!!!!!!!!!!!!
  count_fit=0
  fit_exp_sum=0.0d0
  !!!!!!!!!!!!!!!
  acc_birth=0
  acc_death=0
  DO iter=1,nloop
     !print*,'iter',iter
     !print*,'nexperts',nexperts
     write(12,*) iter,nexperts
     !dimension jumping step
     u=rnunf()
     IF(u>0.5d0 .OR. nexperts==1) THEN
       !BIRTH STEP
       nproposed=min(nexp_max,nexperts+1)
       ALLOCATE(beta_prop(nproposed,nbeta),delta_prop(2,nproposed),&
                prior_ind_prop(nproposed,nchunks),tau_prop(nproposed))
       !generating new parameters beta, delta and tau
       DO j=1,nexperts
          delta_prop(:,j)=delta(:,j)
          beta_prop(j,:)=beta(j,:)
          tau_prop(j)=tau(j)
       END DO
       delta_prop(:,nproposed)=mvrnorm1(iseed,2,mean_d,var_d)
       tau_prop(nproposed)=drnunf()*tau(nexperts)
       var_b=tau_prop(nproposed)*eye(nbeta)
       var_b(1,1)=sigalpha
       beta_prop(nproposed,:)=mvrnorm1(iseed,nbeta,mean_b,var_b)
       DO j=1,nproposed
          prior_ind_prop(j,:)=exp(matmul(xdelta,delta_prop(:,j)))/&
             sum(exp(matmul(xdelta,delta_prop)),2)
       END DO
       DO j=1,nexperts
          prior_ind(j,:)=exp(matmul(xdelta,delta(:,j)))/&
             sum(exp(matmul(xdelta,delta)),2)
       END DO
       !Checking dominance condition
       counter=0
       DO j=1,nproposed
          DO k=1,nproposed
             IF(maxval(prior_ind_prop(j,:)-prior_ind_prop(k,:))<0.0d0) &
             counter=counter+1
          END DO
       END DO
       log_like_curr=log_like(nfreq,nchunks,nbeta,nexperts,xx,y,beta,prior_ind)
       log_like_new=log_like(nfreq,nchunks,nbeta,nproposed,xx,y,beta_prop,&
                             prior_ind_prop)
       met_rat_jump=log_like_new-log_like_curr
       IF(nexperts==1) THEN
           epsilon_jump=min(1.0d0,0.5d0*exp(met_rat_jump))
       ELSE IF (nexperts==nexp_max-1) THEN
           epsilon_jump=min(1.0d0,2.0d0*exp(met_rat_jump))
       ELSE
           epsilon_jump=min(1.0d0,exp(met_rat_jump))
       END IF
       !print*,'epsilon_birth',epsilon_jump,' counter',counter
       u=rnunf()
       !write(26,100) "birth", iter,nexperts, epsilon_jump,counter
       !write(*,100) "birth", iter,nexperts, epsilon_jump,counter
       !100 FORMAT(a5,1x,i5,1x,i2,1x,f10.2,1x,i2)
       IF(epsilon_jump>u .AND. counter==0) THEN
          acc_birth=acc_birth+1
          nexperts=nproposed
          ndelta=2*(nexperts-1)
          DEALLOCATE(beta,delta,tau,prior_ind)
          ALLOCATE(beta(nproposed,nbeta),delta(2,nproposed),tau(nproposed),&
                   prior_ind(nproposed,nchunks)) 
          beta=beta_prop
          delta=delta_prop
          tau=tau_prop
          !CALCULATE NEW PRIOR_IND
          DO j=1,nexperts
             prior_ind(j,:)=exp(matmul(xdelta,delta(:,j)))/&
                sum(exp(matmul(xdelta,delta)),2)
          END DO
       END IF
       DEALLOCATE(tau_prop,beta_prop,delta_prop,prior_ind_prop)
     ELSE
       !DEATH STEP
       nproposed=nexperts-1
       ALLOCATE(beta_prop(nproposed,nbeta),delta_prop(2,nproposed),&
                prior_ind_prop(nproposed,nchunks),tau_prop(nproposed))
       !delta(1,:)=(/0.0, -0.64304866686812/)
       !delta(2,:)=(/0.0,0.02764115636602/)
       !beta(1,:)=(/3.32283832639074, -73.01996795651158, -19.05586827637581,&
       !           -34.81691967298258, -46.57188017605456,-8.91738098595595,&
       !           59.12663838240293,30.88905387793989,10.01531221750966,&
       !           -20.12078557915168, -12.60812700514549/)
       !beta(2,:)=(/2.91184728926154,66.09989179687619,-26.36129212167623,&
       !           -13.96764243259529,3.66334777738740,-16.79573458936665,&
       !           -23.42244397765204,13.10281341818712,22.74223121584767,&
       !           -8.74595277994963,  -8.67254310645845/)
       !tau=(/578.022163287414,566.172917759205/)
       DO j=1,nproposed
          delta_prop(:,j)=delta(:,j)   
          beta_prop(j,:)=beta(j,:)
          tau_prop(j)=tau(j)
       END DO
       !call wrrrn('beta',beta,nra=nexperts)
       !call wrrrn('beta_prop',beta_prop,nra=nproposed)
       !call wrrrn('delta',delta,nra=2)
       !call wrrrn('delta_prop',delta_prop,nra=2)
       !print*,'tau',tau
       !print*,'tau_prop',tau_prop
       DO j=1,nproposed
          prior_ind_prop(j,:)=exp(matmul(xdelta,delta_prop(:,j)))/&
             sum(exp(matmul(xdelta,delta_prop)),2)
       END DO
       DO j=1,nexperts
          prior_ind(j,:)=exp(matmul(xdelta,delta(:,j)))/&
             sum(exp(matmul(xdelta,delta)),2)
       END DO
       !Checking dominance condition
       counter=0
       DO j=1,nproposed
          DO k=1,nproposed
             IF(maxval(prior_ind_prop(j,:)-prior_ind_prop(k,:))<0.0d0) &
             counter=counter+1
          END DO
       END DO
       !print*,'counter',counter
       log_like_curr=log_like(nfreq,nchunks,nbeta,nexperts,xx,y,beta,prior_ind)
       log_like_new=log_like(nfreq,nchunks,nbeta,nproposed,xx,y,beta_prop,&
                             prior_ind_prop)
       !print*,'log_like_curr', log_like_curr
       !print*,'log_like_new', log_like_new
       met_rat_jump=log_like_new-log_like_curr
       IF(nexperts==2) THEN
          epsilon_jump=min(1.0d0,2.0d0*exp(met_rat_jump))
       ELSE IF(nexperts==nexp_max) THEN
          epsilon_jump=min(1.0d0,0.5d0*exp(met_rat_jump))
       ELSE
          epsilon_jump=min(1.0d0,exp(met_rat_jump))
       END IF
       !print*,'epsilon_death',epsilon_jump
       u=rnunf()
       !write(26,105)  "death",iter,nexperts, epsilon_jump, counter
       !105 FORMAT(a5,1x,i5,1x,i2,1x,f10.2,1x,i2)
       IF(epsilon_jump>u .AND. counter==0) THEN
       !IF(epsilon_jump>u) THEN
          acc_death=acc_death+1
          nexperts=nproposed
          ndelta=2*(nexperts-1)
          DEALLOCATE(beta,delta,tau,prior_ind)
          ALLOCATE(beta(nproposed,nbeta),delta(2,nproposed),tau(nproposed),&
                   prior_ind(nproposed,nchunks)) 
          beta=beta_prop
          delta=delta_prop
          tau=tau_prop
          !CALCULATE NEW PRIOR_IND
          DO j=1,nexperts
             prior_ind(j,:)=exp(matmul(xdelta,delta(:,j)))/&
                sum(exp(matmul(xdelta,delta)),2)
          END DO
       END IF
       DEALLOCATE(tau_prop,beta_prop,delta_prop,prior_ind_prop)
     END IF
     !print*,'iteration: ',iter,' nexperts: ',nexperts
     IF(nexperts>1) THEN
       !Drawing the indicators
       ALLOCATE(post_ind(nexperts,nchunks))
       call draw_ind(iseed,nfreq,nchunks,nbeta,nexperts,xx,y,beta,&
                     prior_ind,post_ind,ind)
       !do i=1,nexperts
       !   do j=1,nchunks
       !      write(15,*) post_ind(i,j)
       !   end do
       !end do
       !if(nexperts==4) print*,ind
       !IF(mod(iter,100)==0 .AND. iter>nwarmup) THEN
       !   DO i=1,nchunks
       !      write(25,*) iter,ind(i)
       !   END DO
       !END IF

       DO j=1,nexperts
          kk=0
          where(ind==j) kk=1
          IF(sum(kk)==0) ind(j)=j
       END DO
       !Drawing the deltas 
       call draw_delta(iseed,iter,nchunks,ndelta,nexperts,sigdelta,xdelta,&
                       ind,delta)
     END IF
     !Drawing the betas
     call draw_beta(iseed,nfreq,nchunks,nbeta,nexperts,nbasis,sigalpha,xx,y,&
                    ind,tau,beta)
     call draw_tau(iseed,nbeta,nexperts,nbasis,beta,tau)
     !print*,'tau',tau
     !prior_ind_hat=0.0
     IF(iter>nwarmup) THEN
     DO k=1,nchunks
        fit_chunk=0.0d0
        IF(nexperts>1) THEN
           DO j=1,nexperts
              fit_chunk=fit_chunk+matmul(xx,beta(j,:))*post_ind(j,k)
              !fit_chunk=fit_chunk+matmul(xx,beta(j,:))*prior_ind(j,k)
              !prior_ind_hat=prior_ind_hat+prior_ind/dble(nloop-nwarmup)
           END DO
           fit_ch(:,k)=fit_chunk
        ELSE
           fit_chunk=matmul(xx,beta(1,:))
           fit_ch(:,k)=fit_chunk
        END IF
        fit(:,k)=fit(:,k)+fit_chunk/dble(nloop-nwarmup)
        fit_it(:,iter-nwarmup,k)=fit_chunk
     END DO
     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     !IF(nexperts==3) THEN
       !count_fit=count_fit+1
       !DO j=1,3
       !   fit_exp(:,j)=matmul(xx,beta(j,:))
       !END DO
       !fit_exp_sum=fit_exp_sum+fit_exp
       !IF(mod(iter,10)==0) THEN
       !   DO j=1,nexperts
       !      DO k=1,nchunks
       !         write(45,*) post_ind(j,k)
       !      END DO
       !   END DO
       !END IF
     !ENDIF
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
     fit_ch=exp(fit_ch)
     sum_fit_ch=sum(fit_ch,1)
     DO i=1,nfreq
        fit_ch(i,:)=fit_ch(i,:)/sum_fit_ch
     END DO
     DO k=1,nchunks
        cum_fhat(:,k)=cumsum(nfreq,fit_ch(:,k))
     END DO
     DO i=1,nchunks-1
        !dist(i)=sqrt(sum((cum_fhat(:,i)-cum_fhat(:,i+1))**2)/dble(nfreq))
        dist(i)=maxval(abs(cum_fhat(:,i)-cum_fhat(:,i+1)))
     END DO
     !IF(iter>98000) THEN
     !   DO k=1,nexperts
     !      DO j=1,nchunks
     !         write(24,*) prior_ind(k,j)
     !      END DO
     !   END DO 
     !END IF
   END IF
   IF(nexperts>1) DEALLOCATE(post_ind)
  END DO !end main loop
  !write(*,*) 'acc_birth: ',acc_birth,' acc_death: ', acc_death
  qprop(1)=0.025d0
  qprop(2)=0.975d0
  DO i=1,nfreq
       DO k=1,nchunks
          call eqtil(fit_it(i,:,k),nquan,qprop,quan,xlo,xhi,nmiss)
          lcl(i,k)=quan(1)
          ucl(i,k)=quan(2)
       END DO
  END DO
  DO j=1,nchunks
        DO i=1,nfreq
           write(20,*) lcl(i,j)
           write(22,*) ucl(i,j)
        END DO
  END DO
  !DO k=1,nexperts
  !   DO j=1,nchunks
  !      write(24,*) prior_ind_hat(k,j)
  !   END DO
  !END DO

  !!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !fit_exp_sum=fit_exp_sum/dble(count_fit)
  !DO j=1,3
  !   DO i=1,nfreq
  !      write(40,*) fit_exp_sum(i,j)
  !   END DO
  !END DO
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!
  
  DEALLOCATE(tau,beta,delta,prior_ind)
END SUBROUTINE revjump
!
!
FUNCTION log_spect_ar2(nfreq,comp)
  IMPLICIT none
  INTEGER, INTENT(IN) :: nfreq, comp
  DOUBLE PRECISION, DIMENSION(nfreq) :: log_spect_ar2
  INTEGER :: i
  DOUBLE PRECISION :: phi1,phi2,a,b,c,pi
  DOUBLE PRECISION, DIMENSION(nfreq) :: tp
  pi=4.0d0*atan(1.0d0)
  DO i = 1, nfreq
     tp(i)=(dble(i)-1.0d0)/(2.0d0*(dble(nfreq)-1.0d0))
  END DO 
  IF(comp==1) THEN
    phi1=0.9d0
    phi2=0.0d0
!  ELSE IF (comp==2) THEN
!    phi1=1.69d0
!    phi2=-0.81d0
  ELSE
!    phi1=1.32d0
!    phi2=-0.81d0
    phi1=-0.9d0
    phi2=0.0d0
  END IF
  a=1.0d0+phi1**2+phi2**2
  b=phi1*phi2-phi1
  c=phi2
  log_spect_ar2=-log(a+2.0d0*b*cos(2.0d0*pi*tp)-2.0d0*c*cos(4.0d0*pi*tp))
END FUNCTION log_spect_ar2
!
!
FUNCTION cumsum(n,x)
  IMPLICIT none
  INTEGER, INTENT(IN) :: n
  DOUBLE PRECISION, DIMENSION(n), INTENT(IN) :: x
  DOUBLE PRECISION, DIMENSION(n) :: cumsum
  INTEGER :: i
  cumsum(1)=x(1)
  DO i=2,n
     cumsum(i)=cumsum(i-1)+x(i)
  END DO
END FUNCTION cumsum
!
!
SUBROUTINE draw_ind(iseed,nfreq,nchunks,nbeta,nexperts,xx,y,beta,prior_ind,&
                    post_ind,ind)
  IMPLICIT none
  INTEGER, INTENT(IN) :: iseed,nfreq,nchunks,nbeta,nexperts
  DOUBLE PRECISION, DIMENSION(nfreq,nbeta), INTENT(IN) :: xx
  DOUBLE PRECISION, DIMENSION(nfreq,nchunks), INTENT(IN) :: y
  DOUBLE PRECISION, DIMENSION(nexperts,nchunks), INTENT(IN) :: prior_ind
  DOUBLE PRECISION, DIMENSION(nexperts,nbeta),INTENT(IN) :: beta
  INTEGER, DIMENSION(nchunks), INTENT(OUT) :: ind
  DOUBLE PRECISION, DIMENSION(nexperts,nchunks),INTENT(OUT) :: post_ind
  DOUBLE PRECISION, DIMENSION(nfreq,nexperts) :: ft
  DOUBLE PRECISION, DIMENSION(nfreq) :: ydev
  DOUBLE PRECISION, DIMENSION(nexperts,nchunks) :: loglike
  DOUBLE PRECISION :: u,drnunf
  DOUBLE PRECISION :: sum1,sum2
  INTEGER :: i,j,k
  !Drawing the indicators
  DO k=1,nchunks
    DO j=1,nexperts
        ft(:,j)=matmul(xx,beta(j,:))
        ydev=y(:,k)-ft(:,j)
        loglike(j,k)=sum(ydev(2:nfreq-1)-exp(ydev(2:nfreq-1)))+&
            0.5d0*ydev(1)-0.5d0*exp(ydev(1))+0.5d0*ydev(nfreq)-&
            0.5d0*exp(ydev(nfreq))
        post_ind(j,k)=loglike(j,k)+log(prior_ind(j,k))
    END DO
    sum1=sum(exp(post_ind(:,k)))
    DO j=1,nexperts
        post_ind(j,k)=exp(post_ind(j,k))/sum1
    END DO
    u=drnunf()
    sum2=0.0d0
    DO j=1,nexperts
        IF(u> sum2 .AND. u<post_ind(j,k)+sum2) THEN
           ind(k)=j
        END IF
        sum2=sum2+post_ind(j,k)
    END DO
   END DO !end drawing indicators
END SUBROUTINE draw_ind
!
!
SUBROUTINE draw_delta(iseed,iter,nchunks,ndelta,nexperts,sigdelta,xdelta,ind,&
                      delta)
  USE show_int
  USE rnunf_int
  USE linear_operators
  !USE maxdel_m, ONLY:maxdel
  IMPLICIT none
  INTERFACE
   FUNCTION log_target_delta(ndelta,nchunks,nexperts,xdelta,z,delta,sigdelta)
     USE show_int
     IMPLICIT none
     INTEGER, INTENT(IN) :: ndelta,nchunks,nexperts
     DOUBLE PRECISION, DIMENSION(nchunks,nexperts), INTENT(IN) :: z
     DOUBLE PRECISION, DIMENSION(nchunks,2), INTENT(IN) :: xdelta
     DOUBLE PRECISION, DIMENSION(ndelta), INTENT(IN) :: delta
     DOUBLE PRECISION, INTENT(IN) :: sigdelta
     DOUBLE PRECISION :: log_target_delta
   END FUNCTION
   FUNCTION log_prop(n,x,mean,var_inv)
     IMPLICIT none
     INTEGER, INTENT(IN) :: n
     DOUBLE PRECISION, DIMENSION(n), INTENT(IN) :: x, mean
     DOUBLE PRECISION, DIMENSION(n,n), INTENT(IN) :: var_inv
     DOUBLE PRECISION :: log_prop
   END FUNCTION
   FUNCTION mvrnorm(iseed,n,mean,var)
     USE lin_eig_self_int
     USE linear_operators
     IMPLICIT none
     INTEGER, INTENT(IN) :: n,iseed
     DOUBLE PRECISION, DIMENSION(n), INTENT(IN) :: mean
     DOUBLE PRECISION, DIMENSION(n,n), INTENT(IN) :: var
     DOUBLE PRECISION, DIMENSION(n) :: mvrnorm
    END FUNCTION
    FUNCTION mvrnorm1(iseed,n,mean,var)
      USE linear_operators
      USE show_int
      IMPLICIT none
      INTEGER, INTENT(IN) :: n,iseed
      DOUBLE PRECISION, DIMENSION(n), INTENT(IN) :: mean
      DOUBLE PRECISION, DIMENSION(n,n), INTENT(IN) :: var
      DOUBLE PRECISION, DIMENSION(n) :: mvrnorm1
    END FUNCTION
    FUNCTION mvrnorm2(iseed,n,mean,var) RESULT(r)
      !n is the number of draws from the multivariate normal
      USE rnmvn_int
      USE chfac_int
      IMPLICIT none
      INTEGER, PARAMETER :: DP = kind(1.d0)
      INTEGER :: iseed, n
      REAL(DP) :: mean(:), var(:,:)
      REAL (DP) :: rsig(size(mean),size(mean)), r(n,size(mean))
      INTEGER :: irank
    END FUNCTION
  END INTERFACE

  INTEGER, INTENT(IN) :: nchunks,ndelta,nexperts,iseed
  DOUBLE PRECISION, INTENT(IN) :: sigdelta
  DOUBLE PRECISION, DIMENSION(nchunks,2), INTENT(IN) :: xdelta
  INTEGER, DIMENSION(nchunks), INTENT(IN) :: ind
  INTEGER, INTENT(IN) :: iter
  DOUBLE PRECISION, DIMENSION(2,nexperts), INTENT(INOUT) :: delta
  DOUBLE PRECISION, DIMENSION(2,ndelta) :: temp_delta
  DOUBLE PRECISION, DIMENSION(ndelta) :: temp1
  DOUBLE PRECISION, DIMENSION(ndelta) :: delta_mean
  DOUBLE PRECISION, DIMENSION(ndelta,ndelta) :: delta_var
  DOUBLE PRECISION, DIMENSION(ndelta) :: sold,sold1
  DOUBLE PRECISION, DIMENSION(ndelta) :: derd
  DOUBLE PRECISION, DIMENSION(ndelta,ndelta) :: hessd
  DOUBLE PRECISION, DIMENSION(nchunks,nexperts) :: z
  DOUBLE PRECISION :: met_rat_delta
  INTEGER :: reject
  DOUBLE PRECISION :: tolf,tolx,tol,dmach,u
  DOUBLE PRECISION :: log_target_delta_old,log_target_delta_new
  !DOUBLE PRECISION :: drnunf,epsilond
  DOUBLE PRECISION :: epsilond
  DOUBLE PRECISION :: log_prop_delta_old,log_prop_delta_new
  DOUBLE PRECISION, DIMENSION(ndelta) :: pard
  INTEGER, DIMENSION(nchunks) :: kk
  INTEGER :: irank,nr
  INTEGER :: i,j,k,jj
  INTEGER :: count, counter
  DOUBLE PRECISION :: rnorm(1,ndelta),prior_ind_prop(nexperts,nchunks)
  DOUBLE PRECISION :: delta_prop(2,nexperts)

  !Drawing the deltas
  z=0.0d0
  DO i=1,nchunks
     DO j=1,nexperts
        IF(ind(i)==j) z(i,j)=1.0d0
     END DO
  END DO
  count=0
  DO j=2,nexperts
      DO k=1,2
         count=count+1
         temp_delta(1,count)=delta(k,j)
      END DO
  END DO
  !IF(iter<100) THEN
      !call drnun(ndelta,pard)
      !call drnnor(ndelta,pard)
      !pard=0.1d0*pard
     pard=0.0d0
  !ELSE
  !   pard=temp_delta(1,:)
  !END IF
  !call maxdel(ndelta,pard,ind,nchunks,nexperts,xdelta,sold1)
  !print*,'opt',sold1
  tolf=1.0d0
  tolx=1.0d0
  jj=0
  DO WHILE(tolf>.01d0 .AND. tolx>.01d0 .AND. jj<=20)
     jj=jj+1
     call derivsd(z,ndelta,nexperts,nchunks,pard,xdelta,sigdelta,derd,hessd)
     call dlslsf(ndelta,hessd,ndelta,derd,sold)
     pard=pard-sold
     tolf=sum(abs(derd))
     tolx=sum(abs(-hessd.ix.derd))
  END DO
  delta_mean=pard
  !print*,'nr',delta_mean
  delta_var=-.i.hessd
  rnorm=mvrnorm2(iseed,1,delta_mean,delta_var)
  temp_delta(2,:)=rnorm(1,:)
  !temp_delta(2,:)=mvrnorm1(iseed,ndelta,delta_mean,delta_var)
  log_prop_delta_old=log_prop(ndelta,temp_delta(1,:),delta_mean,-hessd)
  log_prop_delta_new=log_prop(ndelta,temp_delta(2,:),delta_mean,-hessd)
  log_target_delta_old=log_target_delta(ndelta,nchunks,nexperts,xdelta,z,&
                       temp_delta(1,:),sigdelta)
  !print*,'log_target_delta_old',log_target_delta_old
  log_target_delta_new=log_target_delta(ndelta,nchunks,nexperts,xdelta,z,&
                       temp_delta(2,:),sigdelta)
  met_rat_delta=log_target_delta_new-log_target_delta_old+log_prop_delta_old-&
                log_prop_delta_new
  !
  count=0
  delta_prop=0.0
  DO j=2,nexperts
      DO k=1,2
         count=count+1
         delta_prop(k,j)=temp_delta(2,count)
      END DO
  END DO
  DO j=1,nexperts
     prior_ind_prop(j,:)=exp(matmul(xdelta,delta_prop(:,j)))/&
        sum(exp(matmul(xdelta,delta_prop)),2)
  END DO
  !Checking dominance condition
  counter=0
  DO j=1,nexperts
     DO k=1,nexperts
        IF(maxval(prior_ind_prop(j,:)-prior_ind_prop(k,:))<0.0d0) &
        counter=counter+1
     END DO
  END DO
  u=rnunf()
  epsilond=min(1.0d0,exp(met_rat_delta))
  !print*,'epsilond',epsilond,'counter',counter
  !write(26,105)  "gating",epsilond,counter
  !105 FORMAT(a5,1x,f10.2,1x,i2)
  IF(epsilond > u .AND. counter==0) THEN
  count=0
  DO j=2,nexperts
      DO k=1,2
         count=count+1
         delta(k,j)=temp_delta(2,count)
      END DO
  END DO
  END IF
  !IF(u>epsilond) THEN
  !   temp_delta(2,:)=temp_delta(1,:)
  !END IF 
END SUBROUTINE draw_delta
!
!
SUBROUTINE draw_beta(iseed,nfreq,nchunks,nbeta,nexperts,nbasis,sigalpha,xx,y,&
                    ind,tau,beta)
  USE show_int
  USE wrrrn_int
  USE rnunf_int
  USE lin_eig_self_int
  USE linear_operators
  USE maxb_m, ONLY: maxb
  IMPLICIT none
  INTERFACE
    FUNCTION findloc(m,n,ind,ii)
     INTEGER, INTENT(IN) :: m,n
     INTEGER, DIMENSION(n), INTENT(IN) :: ind
     INTEGER, DIMENSION(m) :: findloc
     INTEGER :: ii,j,k
    END FUNCTION
    FUNCTION log_prop(n,x,mean,var_inv)
      IMPLICIT none
      INTEGER, INTENT(IN) :: n
      DOUBLE PRECISION, DIMENSION(n), INTENT(IN) :: x, mean
      DOUBLE PRECISION, DIMENSION(n,n), INTENT(IN) :: var_inv
      DOUBLE PRECISION :: log_prop
   END FUNCTION
   FUNCTION log_target(nfreq,ncoly,nbeta,x,y,par,sigalpha,tau)
     IMPLICIT none
     INTEGER, INTENT(IN) :: nfreq,ncoly,nbeta
     DOUBLE PRECISION, DIMENSION(nfreq,nbeta), INTENT(IN) :: x
     DOUBLE PRECISION, DIMENSION(nfreq,ncoly), INTENT(IN) :: y
     DOUBLE PRECISION, DIMENSION(nbeta), INTENT(IN) :: par
     DOUBLE PRECISION, INTENT(IN) :: sigalpha, tau
     DOUBLE PRECISION :: log_target
   END FUNCTION
   FUNCTION mvrnorm1(iseed,n,mean,var)
     USE linear_operators
     IMPLICIT none
     INTEGER, INTENT(IN) :: n,iseed
     DOUBLE PRECISION, DIMENSION(n), INTENT(IN) :: mean
     DOUBLE PRECISION, DIMENSION(n,n), INTENT(IN) :: var
     DOUBLE PRECISION, DIMENSION(n) :: mvrnorm1
   END FUNCTION
   FUNCTION mvrnorm(iseed,n,mean,var)
     USE lin_eig_self_int
     USE linear_operators
     IMPLICIT none
     INTEGER, INTENT(IN) :: n,iseed
     DOUBLE PRECISION, DIMENSION(n), INTENT(IN) :: mean
     DOUBLE PRECISION, DIMENSION(n,n), INTENT(IN) :: var
     DOUBLE PRECISION, DIMENSION(n) :: mvrnorm
   END FUNCTION
   FUNCTION mvrnorm2(iseed,n,mean,var) RESULT(r)
     !n is the number of draws from the multivariate normal
     USE rnmvn_int
     USE chfac_int
     IMPLICIT none
     INTEGER, PARAMETER :: DP = kind(1.d0)
     INTEGER :: iseed, n
     REAL(DP) :: mean(:), var(:,:)
     REAL (DP) :: rsig(size(mean),size(mean)), r(n,size(mean))
     INTEGER :: irank
   END FUNCTION
  END INTERFACE
  INTEGER, INTENT(IN) :: nfreq,nchunks,nbeta,nexperts,nbasis,iseed
  DOUBLE PRECISION, INTENT(IN) :: sigalpha
  DOUBLE PRECISION, DIMENSION(nfreq,nbeta), INTENT(IN) :: xx
  DOUBLE PRECISION, DIMENSION(nfreq,nchunks), INTENT(IN) :: y
  DOUBLE PRECISION, DIMENSION(nexperts),INTENT(IN) :: tau
  INTEGER, DIMENSION(nchunks),INTENT(IN) :: ind
  DOUBLE PRECISION, DIMENSION(nexperts,nbeta), INTENT(INOUT) :: beta
  DOUBLE PRECISION, DIMENSION(nexperts,nbeta,2) :: temp_beta
  DOUBLE PRECISION, DIMENSION(nfreq,nchunks) :: fit
  DOUBLE PRECISION, DIMENSION(nfreq) :: fit_chunk
  DOUBLE PRECISION, DIMENSION(nbeta) :: derb
  DOUBLE PRECISION, DIMENSION(nbeta) :: param
  DOUBLE PRECISION, DIMENSION(nbeta) :: temp
  DOUBLE PRECISION, DIMENSION(nbeta,nbeta) :: hessb_inv
  DOUBLE PRECISION, DIMENSION(nbeta) :: beta_mean
  DOUBLE PRECISION, DIMENSION(nbeta,nbeta) :: beta_var
  DOUBLE PRECISION, DIMENSION(nbeta) :: solb
  DOUBLE PRECISION, DIMENSION(nbeta) :: sol
  DOUBLE PRECISION, DIMENSION(nbeta) :: rnor
  DOUBLE PRECISION, DIMENSION(nbeta,nbeta) :: hessb,beta_var_inv
  DOUBLE PRECISION, DIMENSION(nbeta,nbeta) :: rsig
  DOUBLE PRECISION, DIMENSION(nfreq,nexperts) :: ft_hat,ft
  DOUBLE PRECISION, DIMENSION(nfreq,nchunks) :: fhat
  DOUBLE PRECISION, DIMENSION(nfreq) :: ydev
  DOUBLE PRECISION, DIMENSION(nexperts,nchunks) :: loglike,post_ind,&
                                                   prior_ind
  DOUBLE PRECISION, DIMENSION(nexperts) :: met_rat
  INTEGER :: reject
  DOUBLE PRECISION, DIMENSION(nexperts,nbeta) :: beta0
  DOUBLE PRECISION :: tolf,tolx,tol,dmach,log_prop_old,log_prop_new,u
  !DOUBLE PRECISION :: log_target_old,log_target_new,epsilonb,drnunf
  DOUBLE PRECISION :: log_target_old,log_target_new,epsilonb
  INTEGER, DIMENSION(nchunks) :: kk
  INTEGER, DIMENSION(nbeta) :: nev
  INTEGER, ALLOCATABLE :: kk_ind(:)
  DOUBLE PRECISION, ALLOCATABLE :: hh(:)
  INTEGER :: i,j,k,iter,ii,jj
  INTEGER :: n_chunk_ind,n_exp
  DOUBLE PRECISION, DIMENSION(nbeta,nbeta) :: eigenvec
  DOUBLE PRECISION, DIMENSION(nbeta) :: eigenval
  DOUBLE PRECISION :: rnorm(1,nbeta)
  temp_beta(:,:,1)=beta
  DO j=1,nexperts
     kk=0
     where(ind.eq.j) kk=1
     n_chunk_ind=sum(kk)
     ALLOCATE(kk_ind(n_chunk_ind))
     !find the chunks belonging to expert j
     kk_ind=findloc(n_chunk_ind,nchunks,ind,j)
     !generate beta
     !find the posterior mode of beta
     !param=temp_beta(j,:,1)
     param=0.0d0
     tolf=1.0d0
     tolx=1.0d0
     ii=0
     DO WHILE(tolf>.001d0 .AND. tolx>.00001d0)
        ii=ii+1
        nev=0
        call derivs(nbeta,nfreq,n_chunk_ind,param,xx,y(:,kk_ind),&
                 sigalpha,tau(j),derb,hessb)
        call lin_eig_self(hessb,eigenval,v=eigenvec)
        WHERE(eigenval.LE.0.0d0) nev=1
        IF(SUM(nev)>0) EXIT
        call dlslsf(nbeta,hessb,nbeta,derb,solb)
        param=param-solb
        tolf=sum(abs(derb))
        tolx=sum(abs(-hessb.ix.derb))
     END DO
     IF(SUM(nev)>0) THEN
        param=temp_beta(j,:,1)
        ALLOCATE(hh(n_chunk_ind))
        hh=1.0d0
        call maxb(param,tau(j),hh,sigalpha,nfreq,n_chunk_ind,nbeta,&
                  nexperts,nbasis,xx,y(:,kk_ind),sol)
        call derivs(nbeta,nfreq,n_chunk_ind,sol,xx,y(:,kk_ind),&
                    sigalpha,tau(j),derb,hessb)
        DEALLOCATE(hh)
     END IF 
     beta_mean=param
     beta_var_inv=-hessb/1.2d0
     beta_var=-1.2d0*.i.hessb
     !temp_beta(j,:,2)=mvrnorm1(iseed,nbeta,beta_mean,beta_var)
     rnorm=mvrnorm2(iseed,1,beta_mean,beta_var)
     !print*,'rnorm',rnorm
     temp_beta(j,:,2)=rnorm(1,:)
     log_prop_old=log_prop(nbeta,temp_beta(j,:,1),beta_mean,beta_var_inv)
     !beta(1,:)=(/2.95326589079431, -52.28456352975171, -59.68971488606051,&
     !-14.03614907332720, -15.24975538988489, -39.33304649805422,&
     ! 47.47076878242155,  43.85075118448089,  16.11726555949620,&
     ! -15.46952991800888, -11.25003783572547/)
     log_prop_new=log_prop(nbeta,temp_beta(j,:,2),beta_mean,beta_var_inv)
     log_target_old=log_target(nfreq,n_chunk_ind,nbeta,xx,y(:,kk_ind),&
                    temp_beta(j,:,1),sigalpha,tau(j))
     log_target_new=log_target(nfreq,n_chunk_ind,nbeta,xx,y(:,kk_ind),&
                    temp_beta(j,:,2),sigalpha,tau(j))
     !Metropolis ratio
     met_rat(j)=log_target_new-log_target_old+log_prop_old-log_prop_new
     u=rnunf()
     epsilonb=min(1.0d0,exp(met_rat(j)))
     IF(u>epsilonb) THEN
         temp_beta(j,:,2)=temp_beta(j,:,1)
         !IF(j==2) reject=reject+1
     END IF 
     beta(j,:)=temp_beta(j,:,2)
     DEALLOCATE(kk_ind)
  END DO !end expert loop
END SUBROUTINE draw_beta
!
!
SUBROUTINE draw_tau(iseed,nbeta,nexperts,nbasis,beta,tau)
  USE show_int
  USE gamin_int
  USE gamdf_int
  USE rngam_int
  USE rnunf_int
  IMPLICIT none
  INTEGER, INTENT(IN) :: nbeta,nexperts,nbasis,iseed
  DOUBLE PRECISION, DIMENSION(nexperts,nbeta), INTENT(IN) :: beta
  DOUBLE PRECISION, DIMENSION(nexperts),INTENT(INOUT) :: tau
  !DOUBLE PRECISION :: taua,taub,gamdev,const1,pinv,dgamin,dgamdf
  DOUBLE PRECISION :: taua,taub,gamdev(1),const1,pinv
  !DOUBLE PRECISION :: u,drnunf
  DOUBLE PRECISION :: u
  INTEGER :: j
  !Drawing tau
  DO j=1,nexperts
     taua=0.5d0*dble(nbasis)
     taub=0.5d0*sum(beta(j,2:nbeta)**2)
     IF (j==1) THEN
         !call drngam(1,taua,gamdev)
         call rngam(taua,gamdev)
         tau(j)=taub/gamdev(1)
     !Drawing Tauj constrained to be less than tauj-1
     ELSE
         !const1=dgamdf(taub/tau(j-1),taua)
         const1=gamdf(taub/tau(j-1),taua)
         !u=drnunf()
         u=rnunf()
         !pinv=u*(1.0d0-const1)+const1
         !tau(j)=taub/dgamin(pinv,taua)
         pinv=1.0d0-u*(1.0d0-const1)
         IF(pinv>0.999d0) THEN
              tau(j)=tau(j-1)-1.0d0
         ELSE
              !tau(j)=taub/dgamin(pinv,taua)
              tau(j)=taub/gamin(pinv,taua)
         END IF
     END IF
  END DO
END SUBROUTINE draw_tau
!
!
SUBROUTINE em(iseed,nfreq,nchunks,nbeta,ndelta,nexperts,nbasis,sigalpha,&
              sigdelta,xx,xdelta,y,ind0,beta0,prior_ind)
  USE show_int
  USE maxb_m, ONLY: maxb
  USE maxd_m, ONLY: maxd
  IMPLICIT none
  INTEGER, INTENT(IN) :: iseed,nfreq,nchunks,nbeta,ndelta,nexperts,nbasis
  DOUBLE PRECISION, DIMENSION(nfreq,nbeta), INTENT(IN) :: xx
  DOUBLE PRECISION, DIMENSION(nchunks,2), INTENT(IN) :: xdelta
  DOUBLE PRECISION, DIMENSION(nfreq,nchunks), INTENT(IN) :: y
  DOUBLE PRECISION, INTENT(IN) :: sigalpha,sigdelta
  DOUBLE PRECISION, DIMENSION(nexperts,nchunks) :: h
  DOUBLE PRECISION, DIMENSION(nchunks) :: h1
  INTEGER, DIMENSION(nchunks), INTENT(OUT) :: ind0
  DOUBLE PRECISION, DIMENSION(nexperts,nbeta), INTENT(OUT) :: beta0
  DOUBLE PRECISION, DIMENSION(nexperts,nchunks), INTENT(OUT) :: prior_ind
  DOUBLE PRECISION, ALLOCATABLE :: tau(:),beta_curr(:,:),delta(:,:)
  DOUBLE PRECISION, ALLOCATABLE :: loglike(:,:)
  DOUBLE PRECISION, ALLOCATABLE :: ydev(:),ft(:,:),denom(:),&
                 h_sum(:),param(:), sol(:),par(:),sold(:)
  INTEGER :: i,j,k,mm,iter,jj,kk
  DOUBLE PRECISION :: f,tolf,fd
  DOUBLE PRECISION :: crit
  ALLOCATE(beta_curr(2,nbeta*nexperts),ft(nfreq,nexperts),ydev(nfreq), &
  loglike(nexperts,nchunks),denom(nchunks),delta(2,nexperts),&
  h_sum(nchunks),tau(nexperts),param(nbeta),sol(nbeta),par(ndelta),&
  sold(ndelta))
  IF(nexperts>=2) THEN
  beta_curr=0.0d0
  tau(1)=10000.0d0
  DO j=2,nexperts
     tau(j)=tau(j-1)/10.0d0
  END DO
  delta=0.0d0
  iter=0
  crit=100.0d0
  !DO mm=1,5 !start EM loop
  DO WHILE(crit>.00001d0 .AND. iter <=20) !start EM loop
     iter=iter+1
     !E-step
     DO k=1,nchunks
        DO j=1,nexperts
           ft(:,j)=matmul(xx,beta_curr(1,1+(j-1)*nbeta:j*nbeta))
           ydev=y(:,k)-ft(:,j)
           loglike(j,k)=sum(ydev(2:nfreq-1)-exp(ydev(2:nfreq-1)))+&
                        0.5*(ydev(1)-exp(ydev(1)))+&
                        0.5*(ydev(nfreq)-exp(ydev(nfreq)))
        END DO
     END DO
     denom=0.0d0
     DO j=1,nexperts
        denom=denom+exp(matmul(xdelta,delta(:,j)))
     END DO 
     DO j=1,nexperts
        prior_ind(j,:)=exp(matmul(xdelta,delta(:,j)))/denom
     END DO
     h=prior_ind*exp(loglike)
     h_sum=sum(h,1)
     DO j=1,nexperts
        h(j,:)=h(j,:)/h_sum
     END DO
     IF(iter==1) THEN
        DO j=1,nexperts
           call drnun(nchunks,h(j,:))
        END DO
        h_sum=sum(h,1)
        DO j=1,nexperts
           h(j,:)=h(j,:)/h_sum
        END DO
     END IF
     !M-step
     !wrt the betas
     DO j=1,nexperts
        param=beta_curr(1,1+(j-1)*nbeta:j*nbeta)
        call maxb(param,tau(j),h(j,:),sigalpha,nfreq,nchunks,nbeta,&
                  nexperts,nbasis,xx,y,sol)
        beta_curr(2,1+(j-1)*nbeta:j*nbeta)=sol
     END DO
     !call show(beta_curr)

     !wrt the deltas
     par=0.0d0
     call maxd(ndelta,par,h,sigdelta,nchunks,nexperts,xdelta,sold)
     ! print*,sold
     DO j=1,nexperts-1
        delta(:,j+1)=sold(1+(j-1)*2:j*2)
     END DO
     !wrt the taus
     DO j=1,nexperts
        tau(j)=0.5d0*sum(beta_curr(2,2+(j-1)*nbeta:j*nbeta)**2)
     END DO
     !print*,tau
     crit=maxval(abs(beta_curr(1,:)-beta_curr(2,:))/(beta_curr(1,:)+.001d0))
     beta_curr(1,:)=beta_curr(2,:)     
   END DO !end of EM loop
   DO j=1,nexperts
      beta0(j,:)=beta_curr(1,1+(j-1)*nbeta:j*nbeta)
   END DO
   ind0=maxloc(h,dim=1)

   ELSE
     h1=1.0d0
     call maxb(param,tau(j),h1,sigalpha,nfreq,nchunks,nbeta,&
               nexperts,nbasis,xx,y,sol)
     beta0(1,:)=sol
     ind0=1
     prior_ind(1,:)=1.0d0    
   END IF
   DEALLOCATE(beta_curr,ft,loglike,ydev,denom,delta,h_sum,param,&
              tau,sol,par,sold)
END SUBROUTINE em
!
!
  SUBROUTINE funcb(nbeta,param,f)
  USE maxb_m, ONLY:tau_exp,pp,sigalpha,nfreq,nchunks,nexperts,nbasis,xx,y
  IMPLICIT none
  INTEGER, INTENT(IN) :: nbeta
  DOUBLE PRECISION, DIMENSION(nbeta), INTENT(IN) :: param
  DOUBLE PRECISION, INTENT(OUT) :: f
  DOUBLE PRECISION, ALLOCATABLE :: ydev(:,:),xxparam(:)
  INTEGER :: i,j,k

  ALLOCATE(xxparam(nfreq),ydev(nfreq,nchunks))
  xxparam=matmul(xx,param)
  DO j=1,nchunks
     ydev(:,j)=y(:,j)-xxparam
  END DO
  f=0.0d0
  DO k=1,nchunks
     f=f+pp(k)*(0.5*ydev(1,k)-0.5*exp(ydev(1,k)) &
       +0.5*ydev(nfreq,k)-0.5*exp(ydev(nfreq,k)) &
       +sum(ydev(2:nfreq-1,k)-exp(ydev(2:nfreq-1,k))))
  END DO 
  f=f-0.5*(sum(param(2:nbeta)**2)/tau_exp+param(1)**2/sigalpha);
  f=-f
  DEALLOCATE(xxparam,ydev)
END SUBROUTINE funcb
!
!
SUBROUTINE gradb(nbeta,param,g)
  USE maxb_m, ONLY:tau_exp,pp,sigalpha,nfreq,nchunks,nexperts,xx,y
  IMPLICIT none
  INTEGER, INTENT(IN) :: nbeta
  DOUBLE PRECISION, DIMENSION(nbeta), INTENT(IN) :: param
  DOUBLE PRECISION, DIMENSION(nbeta),INTENT(OUT) :: g
  INTEGER :: i,j,k

  g=0.0d0
  g(1)=-param(1)/sigalpha
  g(2:nbeta)=-param(2:nbeta)/tau_exp
  DO i=1,nfreq
     DO k=1,nchunks
             IF(i>1 .AND. i<nfreq) THEN
                 g=g-pp(k)*xx(i,:)+pp(k)*xx(i,:)*exp(y(i,k)-sum(xx(i,:)*param))
             ELSE
                 g=g-0.5d0*pp(k)*xx(i,:)+0.5*pp(k)*xx(i,:) &
                   *exp(y(i,k)-sum(xx(i,:)*param))
             END IF
     END DO
  END DO 
  g=-g
END SUBROUTINE gradb
!
!
SUBROUTINE hessb(nbeta,param,hess,ldh)
  USE maxb_m, ONLY:tau_exp,pp,sigalpha,nfreq,nchunks,nexperts,nbasis,xx,y
  IMPLICIT none
  INTERFACE
    FUNCTION out_prod(n,x)
     INTEGER, INTENT(IN) :: n
     DOUBLE PRECISION, DIMENSION(n), INTENT(IN) :: x
     DOUBLE PRECISION, DIMENSION(n,n) :: out_prod
    END FUNCTION
  END INTERFACE
  INTEGER, INTENT(IN) :: nbeta,ldh
  DOUBLE PRECISION, DIMENSION(nbeta), INTENT(IN) :: param
  DOUBLE PRECISION, DIMENSION(nbeta,nbeta),INTENT(OUT) :: hess
  INTEGER :: i,j,k
  hess=0.0d0
  hess(1,1)=-1.0/sigalpha
  hess(2:nbeta,2:nbeta)=-1.0d0/tau_exp
  DO i=1,nfreq
     DO k=1,nchunks
             IF(i>1 .AND. i<nfreq) THEN
                 hess=hess-pp(k)*out_prod(nbeta,xx(i,:))*&
                      exp(y(i,k)-sum(xx(i,:)*param))
             ELSE
                 hess=hess-pp(k)*out_prod(nbeta,xx(i,:))*&
                      0.5*exp(y(i,k)-sum(xx(i,:)*param))
             END IF
     END DO
  END DO 
  hess=-hess
END SUBROUTINE hessb
!
!
FUNCTION out_prod(n,x)
  IMPLICIT none
  INTEGER, INTENT(IN) :: n
  DOUBLE PRECISION, DIMENSION(n), INTENT(IN) :: x
  DOUBLE PRECISION, DIMENSION(n,n) :: out_prod
  INTEGER :: i

  DO i=1,n
     out_prod(i,:)=x(i)*x
  END DO
END FUNCTION out_prod
!
!
SUBROUTINE funcd(ndelta,param,f)
  USE maxd_m, ONLY:h,sigdelta,nchunks,nexperts,xdelta
  IMPLICIT none
  INTEGER, INTENT(IN) :: ndelta
  DOUBLE PRECISION, DIMENSION(ndelta), INTENT(IN) :: param
  DOUBLE PRECISION, INTENT(OUT) :: f
  DOUBLE PRECISION, DIMENSION(nexperts*2) :: par1
  DOUBLE PRECISION, DIMENSION(nchunks) :: denom
  DOUBLE PRECISION, DIMENSION(nexperts,nchunks) :: prior_ind
  INTEGER :: i,j,k
  par1=0.0d0
  par1(3:2*nexperts)=param
  denom=0.0d0
  DO j=1,nexperts
     denom=denom+exp(matmul(xdelta,par1(1+(j-1)*2:2*j)))
  END DO 
  DO j=1,nexperts
     prior_ind(j,:)=exp(matmul(xdelta,par1(1+(j-1)*2:2*j)))/denom
  END DO
  f=SUM(h*log(prior_ind))-0.5d0*SUM(param**2)/sigdelta
  f=-f
END SUBROUTINE funcd
!
!
SUBROUTINE gradd(ndelta,param,g)
  USE maxd_m, ONLY:h,sigdelta,nchunks,nexperts,xdelta
  IMPLICIT none
  INTEGER, INTENT(IN) :: ndelta
  DOUBLE PRECISION, DIMENSION(ndelta), INTENT(IN) :: param
  DOUBLE PRECISION, DIMENSION(ndelta),INTENT(OUT) :: g
  DOUBLE PRECISION, DIMENSION(nexperts*2) :: par1
  DOUBLE PRECISION, DIMENSION(nchunks) :: denom
  DOUBLE PRECISION, DIMENSION(nexperts,nchunks) :: prior_ind
  INTEGER :: i,j,k
  par1=0.0d0
  par1(3:2*nexperts)=param
  denom=0.0d0
  DO j=1,nexperts
     denom=denom+exp(matmul(xdelta,par1(1+(j-1)*2:2*j)))
  END DO 
  DO j=1,nexperts
     prior_ind(j,:)=exp(matmul(xdelta,par1(1+(j-1)*2:2*j)))/denom
  END DO
  g=-par1(3:2*nexperts)/sigdelta
  DO k=1,nexperts-1
     DO i=1,nchunks
        g(1+(k-1)*2:k*2)=g(1+(k-1)*2:k*2)+(h(k,i)-prior_ind(k,i))*xdelta(i,:)
     END DO
  END DO
  g=-g
END SUBROUTINE gradd
!
!
SUBROUTINE hessd(ndelta,param,hd,ldh)
  USE maxd_m, ONLY:h,sigdelta,nchunks,nexperts,xdelta
  IMPLICIT none
  INTERFACE
    FUNCTION out_prod(n,x)
     INTEGER, INTENT(IN) :: n
     DOUBLE PRECISION, DIMENSION(n), INTENT(IN) :: x
     DOUBLE PRECISION, DIMENSION(n,n) :: out_prod
    END FUNCTION
  END INTERFACE
  INTEGER, INTENT(IN) :: ndelta,ldh
  DOUBLE PRECISION, DIMENSION(ndelta), INTENT(IN) :: param
  DOUBLE PRECISION, DIMENSION(ndelta,ndelta),INTENT(OUT) :: hd
  DOUBLE PRECISION, DIMENSION(nexperts*2) :: par1
  DOUBLE PRECISION, DIMENSION(nchunks) :: denom
  DOUBLE PRECISION, DIMENSION(nexperts,nchunks) :: prior_ind
  INTEGER :: i,j,k,l,kron
  par1=0.0d0
  par1(3:2*nexperts)=param
  denom=0.0d0
  DO j=1,nexperts
     denom=denom+exp(matmul(xdelta,par1(1+(j-1)*2:2*j)))
  END DO 
  DO j=1,nexperts
     prior_ind(j,:)=exp(matmul(xdelta,par1(1+(j-1)*2:2*j)))/denom
  END DO
  hd=0.0d0
  DO i=1,(nexperts-1)*2
     hd(i,i)=-1.0d0/sigdelta
  END DO
  DO k=1,nexperts-1
     DO l=1,nexperts-1
        IF(k==l) THEN
           kron=1
        ELSE
           kron=0
        END IF
        DO i=1,nchunks
           hd(1+(k-1)*2:2*k,1+(l-1)*2:2*l)=hd(1+(k-1)*2:2*k,1+(l-1)*2:2*l)-&
           prior_ind(k,i)*(kron-prior_ind(l,i))*out_prod(2,xdelta(i,:))
        END DO
     END DO
  END DO
  hd=-hd
END SUBROUTINE hessd
!
!
FUNCTION findloc(m,n,ind,ii)
  !find the locations of all the values of
  !ind equal to ii
  IMPLICIT none
  INTEGER, INTENT(IN) :: m,n
  INTEGER, DIMENSION(n), INTENT(IN) :: ind
  INTEGER, DIMENSION(m) :: findloc
  INTEGER :: ii,j,k
  k=1
  DO j=1,n
     IF(ind(j).eq. ii) THEN
        findloc(k)=j
        k=k+1
     END IF
  END DO
END FUNCTION findloc
!
!
SUBROUTINE derivs(nbeta,nfreq,ncoly,beta,x,y,sigalpha,tau,der,hess)
   USE show_int
   IMPLICIT none
   INTERFACE
    FUNCTION out_prod(n,x)
     INTEGER, INTENT(IN) :: n
     DOUBLE PRECISION, DIMENSION(n), INTENT(IN) :: x
     DOUBLE PRECISION, DIMENSION(n,n) :: out_prod
     END FUNCTION
   END INTERFACE
   INTEGER, INTENT(IN) :: nbeta,nfreq,ncoly
   DOUBLE PRECISION, DIMENSION(nbeta), INTENT(IN) :: beta
   DOUBLE PRECISION, DIMENSION(nfreq,nbeta), INTENT(IN) :: x
   DOUBLE PRECISION, DIMENSION(nfreq,ncoly), INTENT(IN) :: y
   DOUBLE PRECISION, INTENT(IN) :: sigalpha,tau
   DOUBLE PRECISION, DIMENSION(nbeta), INTENT(OUT) :: der
   DOUBLE PRECISION, DIMENSION(nbeta,nbeta), INTENT(OUT) :: hess
   INTEGER :: i,k
   DOUBLE PRECISION :: exp_ys
   der=0.0d0
   der(1)=-beta(1)/sigalpha
   der(2:nbeta)=-beta(2:nbeta)/tau
   hess=0.0d0
   hess(1,1)=-1.0d0/sigalpha
   DO i=2,nbeta
      hess(i,i)=-1.0d0/tau
   END DO
   DO i=1,nfreq
      DO k=1,ncoly
         !exp_ys=exp(y(i,k)-sum(x(i,:)*beta))
         IF(i>1 .and. i<nfreq) THEN
            der=der-x(i,:)+x(i,:)*exp(y(i,k)-sum(x(i,:)*beta))
            hess=hess-out_prod(nbeta,x(i,:))*exp(y(i,k)-sum(x(i,:)*beta))
         ELSE
            der=der-0.5d0*x(i,:)+0.5d0*x(i,:)*exp(y(i,k)-sum(x(i,:)*beta))
            hess=hess-0.5d0*out_prod(nbeta,x(i,:))*&
                           exp(y(i,k)-sum(x(i,:)*beta))
         END IF
      END DO
   END DO
END SUBROUTINE derivs
!
!
SUBROUTINE derivs1(nbeta,nfreq,ncoly,beta,x,y,sigalpha,tau,der,hess)
   USE show_int
   IMPLICIT none
   INTERFACE
    FUNCTION out_prod(n,x)
     INTEGER, INTENT(IN) :: n
     DOUBLE PRECISION, DIMENSION(n), INTENT(IN) :: x
     DOUBLE PRECISION, DIMENSION(n,n) :: out_prod
    END FUNCTION
   END INTERFACE
   INTEGER, INTENT(IN) :: nbeta,nfreq,ncoly
   DOUBLE PRECISION, DIMENSION(nbeta), INTENT(IN) :: beta
   DOUBLE PRECISION, DIMENSION(nfreq,nbeta), INTENT(IN) :: x
   DOUBLE PRECISION, DIMENSION(nfreq,ncoly), INTENT(IN) :: y
   DOUBLE PRECISION, INTENT(IN) :: sigalpha,tau
   DOUBLE PRECISION, DIMENSION(nbeta), INTENT(OUT) :: der
   DOUBLE PRECISION, DIMENSION(nbeta,nbeta), INTENT(OUT) :: hess
   INTEGER :: i,k
   DOUBLE PRECISION :: exp_ys
   der=0.0d0
   der(1)=-beta(1)/sigalpha
   der(2:nbeta)=-beta(2:nbeta)/tau
   hess=0.0d0
   hess(1,1)=-1.0d0/sigalpha
   DO i=2,nbeta
      hess(i,i)=-1.0d0/tau
   END DO
   call show(y)
   DO i=1,nfreq
      DO k=1,ncoly
         IF (mod(i,nfreq)==0.OR.mod(i,nfreq)==1) THEN
            der=der-0.5d0*x(i,:)+0.5d0*x(i,:)*exp(y(i,k)-sum(x(i,:)*beta))
            hess=hess-out_prod(nbeta,x(i,:))*0.5d0*&
                            exp(y(i,k)-sum(x(i,:)*beta))
         ELSE
            der=der-x(i,:)+x(i,:)*exp(y(i,k)-sum(x(i,:)*beta))
            hess=hess-out_prod(nbeta,x(i,:))*exp(y(i,k)-sum(x(i,:)*beta))
         END IF
      END DO
   END DO
END SUBROUTINE derivs1
!
!
SUBROUTINE derivsd(z,ndelta,nexperts,nchunks,param,xdelta,sigdelta,der,hess)
   USE show_int
   USE linear_operators
   IMPLICIT none
   INTERFACE
     FUNCTION out_prod(n,x)
      INTEGER, INTENT(IN) :: n
      DOUBLE PRECISION, DIMENSION(n), INTENT(IN) :: x
      DOUBLE PRECISION, DIMENSION(n,n) :: out_prod
     END FUNCTION
   END INTERFACE
   INTEGER, INTENT(IN) :: nexperts,nchunks,ndelta
   !INTEGER, INTENT(IN),DIMENSION(nchunks) :: ind
   DOUBLE PRECISION, INTENT(IN),DIMENSION(nchunks,nexperts) :: z
   DOUBLE PRECISION, DIMENSION(nchunks,2), INTENT(IN) :: xdelta
   DOUBLE PRECISION, DIMENSION(ndelta), INTENT(IN) :: param
   DOUBLE PRECISION, INTENT(IN) :: sigdelta
   DOUBLE PRECISION, DIMENSION(ndelta), INTENT(OUT) :: der
   DOUBLE PRECISION, DIMENSION(ndelta,ndelta), INTENT(OUT) :: hess
   DOUBLE PRECISION, DIMENSION(nexperts,nchunks) :: prior_ind
   DOUBLE PRECISION, DIMENSION(nchunks) :: denom
   DOUBLE PRECISION, DIMENSION(nexperts*2) :: param1
   DOUBLE PRECISION :: kron
   INTEGER :: i,j,k,l
   param1=0.0d0
   param1(3:nexperts*2)=param
   denom=0.0d0
   DO j=1,nexperts
      denom=denom+exp(matmul(xdelta,param1(1+(j-1)*2:2*j)))
   END DO 
   DO j=1,nexperts
      prior_ind(j,:)=exp(matmul(xdelta,param1(1+(j-1)*2:2*j)))/denom
   END DO
   der=-param/sigdelta
   hess=-1.0d0/sigdelta*eye(ndelta)
   DO k=2,nexperts
      DO i=1,nchunks 
        der(1+(k-2)*2:2*(k-1))=der(1+(k-2)*2:2*(k-1))+(z(i,k)-prior_ind(k,i))*&
        xdelta(i,:)
      END DO
   END DO
   DO k=2,nexperts
      DO l=2,nexperts
         IF(k==l) THEN
           kron=1.0d0
         ELSE
           kron=0.0d0
         END IF
         DO i=1,nchunks
            hess(1+(k-2)*2:2*(k-1),1+(l-2)*2:2*(l-1))=&
            hess(1+(k-2)*2:2*(k-1),1+(l-2)*2:2*(l-1))-prior_ind(k,i)*&
            (kron-prior_ind(l,i))*out_prod(2,xdelta(i,:))
         END DO
      END DO
   END DO
END SUBROUTINE derivsd
!
!
FUNCTION log_prop(n,x,mean,var_inv)
  IMPLICIT none
  INTEGER, INTENT(IN) :: n
  DOUBLE PRECISION, DIMENSION(n), INTENT(IN) :: x, mean
  DOUBLE PRECISION, DIMENSION(n,n), INTENT(IN) :: var_inv
  DOUBLE PRECISION :: log_prop
                                                                                
  DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: xm
  DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: xmtv
                                                                                
  ALLOCATE(xm(n),xmtv(n))
  xm=x-mean
  xmtv=matmul(xm,var_inv)
  log_prop=-0.5d0*sum(xmtv*xm)
  DEALLOCATE(xm,xmtv)
END FUNCTION log_prop

FUNCTION log_target(nfreq,ncoly,nbeta,x,y,par,sigalpha,tau)
  IMPLICIT none
  INTEGER, INTENT(IN) :: nfreq,ncoly,nbeta
  DOUBLE PRECISION, DIMENSION(nfreq,nbeta), INTENT(IN) :: x
  DOUBLE PRECISION, DIMENSION(nfreq,ncoly), INTENT(IN) :: y
  DOUBLE PRECISION, DIMENSION(nbeta), INTENT(IN) :: par
  DOUBLE PRECISION, INTENT(IN) :: sigalpha, tau
  DOUBLE PRECISION :: log_target
                                                                                
  DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: xpar
  DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: ys
  INTEGER :: j
                                                                                
  ALLOCATE(xpar(nfreq),ys(nfreq,ncoly))
                                                                                
  xpar=matmul(x,par)
  DO j=1,ncoly
     ys(:,j)=y(:,j)-xpar
  END DO
  log_target=sum(ys(2:(nfreq-1),:) - exp(ys(2:(nfreq-1),:)))+ &
              0.5d0*sum(ys(1,:)-exp(ys(1,:))) + &
              0.5d0*sum(ys(nfreq,:)-exp(ys(nfreq,:)))- &
              0.5d0*par(1)**2/sigalpha-0.5d0*sum(par(2:nbeta)**2)/tau
   DEALLOCATE(xpar,ys)
END FUNCTION log_target
!
!
FUNCTION log_like(nfreq,nchunks,nbeta,nexp,xx,y,beta,prior_ind)
  USE wrrrn_int
  IMPLICIT none
  INTEGER, INTENT(IN) :: nfreq,nchunks,nbeta,nexp
  DOUBLE PRECISION, DIMENSION(nfreq,nbeta), INTENT(IN) :: xx
  DOUBLE PRECISION, DIMENSION(nfreq,nchunks), INTENT(IN) :: y
  DOUBLE PRECISION, DIMENSION(nexp,nbeta), INTENT(IN) :: beta
  DOUBLE PRECISION, DIMENSION(nexp,nchunks), INTENT(IN) :: prior_ind
  DOUBLE PRECISION :: log_like
  DOUBLE PRECISION, DIMENSION(nfreq,nchunks,nexp) :: ys
  DOUBLE PRECISION, DIMENSION(nfreq,nexp) :: xb
  DOUBLE PRECISION, DIMENSION(nfreq,nchunks,nexp) :: log_like_arr
  DOUBLE PRECISION, DIMENSION(nfreq,nchunks) :: like_mat
  INTEGER :: i,j,k
                                                                                
  !print*,'log_like func'
  !call wrrrn('beta',beta,nra=nexp)
  !call wrrrn('xx',xx,nra=nfreq)
  DO j=1,nexp
     xb(:,j)=matmul(xx,beta(j,:))
  END DO
  DO j=1,nexp
     DO k=1,nchunks
        ys(:,k,j)=y(:,k)-xb(:,j)
     END DO
  END DO
  !IF(nexp==1) THEN
  !   call wrrrn('ys',ys(:,:,1),nra=nfreq)
  !END IF
  DO j=1,nexp
     log_like_arr(2:nfreq-1,:,j)=ys(2:nfreq-1,:,j)-exp(ys(2:nfreq-1,:,j))
     log_like_arr(1,:,j)=0.5d0*(ys(1,:,j)-exp(ys(1,:,j)))
     log_like_arr(nfreq,:,j)=0.5d0*(ys(nfreq,:,j)-exp(ys(nfreq,:,j)))
  END DO
  !call wrrrn('log_like_arr',log_like_arr(:,:,1),nra=nfreq)
  like_mat=0.0d0
  DO k=1,nchunks
     DO j=1,nexp
        like_mat(:,k)=like_mat(:,k)+prior_ind(j,k)*exp(log_like_arr(:,k,j))
     END DO
  END DO
  log_like=sum(log(like_mat))
END FUNCTION log_like
!
!
FUNCTION log_target_delta(ndelta,nchunks,nexperts,xdelta,z,delta,sigdelta)
  USE show_int
  IMPLICIT none
  INTEGER, INTENT(IN) :: ndelta,nchunks,nexperts
  DOUBLE PRECISION, DIMENSION(nchunks,nexperts), INTENT(IN) :: z
  DOUBLE PRECISION, DIMENSION(nchunks,2), INTENT(IN) :: xdelta
  DOUBLE PRECISION, DIMENSION(ndelta), INTENT(IN) :: delta
  DOUBLE PRECISION, INTENT(IN) :: sigdelta
  DOUBLE PRECISION :: log_target_delta
  DOUBLE PRECISION, DIMENSION(nchunks,nexperts) :: prob
  DOUBLE PRECISION, DIMENSION(nchunks,nexperts) :: temp
  INTEGER :: i,j,k
  temp=0.0d0
  DO j=2,nexperts
     temp(:,j)=matmul(xdelta,delta(1+(j-2)*2:2*(j-1)))
  END DO
  DO j=1,nexperts
     prob(:,j)=exp(temp(:,j))/sum(exp(temp),2)
  END DO
  log_target_delta=sum(z*log(prob))-0.5d0*sum(delta**2)/sigdelta
END FUNCTION log_target_delta
!
!
FUNCTION mvrnorm(iseed,n,mean,var)
  USE lin_eig_self_int
  USE linear_operators
  IMPLICIT none
  INTEGER, INTENT(IN) :: n,iseed
  DOUBLE PRECISION, DIMENSION(n), INTENT(IN) :: mean
  DOUBLE PRECISION, DIMENSION(n,n), INTENT(IN) :: var
  DOUBLE PRECISION, DIMENSION(n) :: mvrnorm
  DOUBLE PRECISION, DIMENSION(n) :: ev
  DOUBLE PRECISION, DIMENSION(n) :: z
  DOUBLE PRECISION, DIMENSION(n,n) :: evc
  DOUBLE PRECISION, DIMENSION(n) :: temp
  !find eigen value decomposition of var
  call lin_eig_self(var,ev,v=evc)
  !draw n independent N(0,1)
  call drnnor(n,z)
  temp=evc .x. diag(sqrt(ev)) .x. z
  mvrnorm=mean+temp
END FUNCTION mvrnorm
!
!
FUNCTION mvrnorm1(iseed,n,mean,var)
  USE linear_operators
  USE show_int
  IMPLICIT none
  INTEGER, INTENT(IN) :: n,iseed
  DOUBLE PRECISION, DIMENSION(n), INTENT(IN) :: mean
  DOUBLE PRECISION, DIMENSION(n,n), INTENT(IN) :: var
  DOUBLE PRECISION, DIMENSION(n) :: mvrnorm1
  DOUBLE PRECISION, DIMENSION(n) :: z
  DOUBLE PRECISION, DIMENSION(n) :: temp
  DOUBLE PRECISION, DIMENSION(n,n) :: fac
  DOUBLE PRECISION :: rcond
  INTEGER :: i,j
  !Find Cholesky decomposition
  call dlfcds(n,var,n,fac,n,rcond)
  DO i=2,n
     DO j=1,i-1
        fac(i,j)=0.0d0
     END DO
  END DO
  !draw n independent N(0,1)
  call drnnor(n,z)
  temp=fac.tx.z
  mvrnorm1=mean+temp
END FUNCTION mvrnorm1
!
FUNCTION mvrnorm2(iseed,n,mean,var) RESULT(r)
  !n is the number of draws from the multivariate normal
  USE rnmvn_int
  USE chfac_int
  IMPLICIT none
  INTEGER, PARAMETER :: DP = kind(1.d0)
  INTEGER :: iseed, n
  REAL(DP) :: mean(:), var(:,:)
  REAL (DP) :: rsig(size(mean),size(mean)), r(n,size(mean))
  INTEGER :: irank
  INTEGER :: i
  call chfac(var,irank,rsig)
  call rnmvn(rsig,r)
  DO i=1,n
     r(i,:)=r(i,:)+mean
  END DO
END FUNCTION mvrnorm2
!
!FUNCTION euclid(n,x)
!  IMPLICIT none
!  INTEGER, INTENT(IN) :: n
!  DOUBLE PRECISION, DIMENSION(n), INTENT(IN) :: x
!END FUNCTION euclid