SUBROUTINE CODE3_LOOP
USE CODE8
USE CODE4
USE CODE7, ONLY : MATINV
USE CODE5, ONLY : Sample_Truncated_Normal,Sample_Truncated_Mixture_Normal,Sample_Multivariate_Normal
IMPLICIT NONE
REAL(8) :: tr,mn_zi,max_y,one,iff(Hh,Hh),iff_o(Hh,Hh),Tf,mf,a1,a2,vf(Hh,Hh), &
        vf_sub(Hh-1,Hh-1),mf_c,tf_c,zero,fall(Ii*Tt,Hh),time,time2, &
        fs(Hh),u,rho,countf,crap
INTEGER :: n,cci,aa,o,j,st,ind(Hh-1),k,k2,k3,k4,hz,initt,finish,hz2,in2,fin2, &
        s,ss,f_accepted(Hh)
LOGICAL :: lstep,lburnin,logc

one=1.0d0
zero=0.0d0
CALL system_clock(count_rate=hz)
CALL system_clock(count=initt)

! Initial values
cf = cf +1
IF (n_p>0) THEN
        DO np = 1, n_p
                DO i = 1, Ii
                        IF (pr(np)%eq%m(i)==1.0d0) THEN
                                logc = pr(np)%d(i)==1.0d0
                                mn_zi = DOT_PRODUCT(pr(np)%eq%X_all(i,:),pr(np)%eq%b) + DOT_PRODUCT(f(i,:),pr(np)%eq%a)
                                pr(np)%eq%Y_all(i) = Sample_Truncated_Normal(mn_zi,1.0d0/pr(np)%eq%t,zero,logc)
                        END IF
                END DO
                WRITE(4,fmt=*) 'DONE INITIALIZING LATENT IN PROBIT',np
        END DO
END IF
CALL system_clock(count_rate=hz2)
CALL system_clock(count=in2)
DO j = 1, (step*iterations + burnin) 
        IF (n_p>0) THEN
                DO np = 1, n_p
                        pr(np)%eq%coun=0
                        pr(np)%eq%f = f(pr(np)%eq%i,:)
                        pr(np)%eq%Y = pr(np)%eq%Y_all(pr(np)%eq%i)
                        pr(np)%eq%Yq = pr(np)%eq%Y - MATMUL(pr(np)%eq%x,pr(np)%eq%b)
                        pr(np)%eq%Yq = pr(np)%eq%Yq - MATMUL(pr(np)%eq%f,pr(np)%eq%a)
                        pr(np)%eq%xmix = pr(np)%eq%x
                ! sample a|b
                        IF (pr(np)%eq%nor%l_f) THEN
                        pr(np)%eq%f_f = pr(np)%eq%f(:,pr(np)%eq%nor%i_f)
                                pr(np)%eq%Yq = pr(np)%eq%Y - MATMUL(pr(np)%eq%f,pr(np)%eq%nor%zo_f*pr(np)%eq%a)
                                pr(np)%eq%Yq = pr(np)%eq%Yq - MATMUL(pr(np)%eq%xmix,pr(np)%eq%b) 
                                CALL SAMPLE_ALPHA(pr(np)%eq,prior_var)                  
                        END IF
                        IF (pr(np)%eq%nor%l_s) THEN
                                pr(np)%eq%f_s = pr(np)%eq%f(:,pr(np)%eq%nor%i_s)
                                pr(np)%eq%Yq = pr(np)%eq%Y - MATMUL(pr(np)%eq%f,pr(np)%eq%nor%zo_s*pr(np)%eq%a)
                                pr(np)%eq%Yq = pr(np)%eq%Yq - MATMUL(pr(np)%eq%xmix,pr(np)%eq%b) 
                                CALL SAMPLE_ALPHA_TRUNCATED(pr(np)%eq,prior_var)
                        END IF  
                        ! sample b|a
                        pr(np)%eq%Yq = pr(np)%eq%Y - MATMUL(pr(np)%eq%f,pr(np)%eq%a)
                        CALL SAMPLE_BETA(pr(np)%eq)
                                st=0
                                WRITE(pr(np)%eq%file,'(40F16.8)') pr(np)%eq%b,pr(np)%eq%a,pr(np)%eq%t
                        IF (j==1) WRITE(4,fmt=*) 'DONE WITH PROBIT',np,' ITERATION 1' 
                END DO
        END IF  
        IF (n_e>0) THEN
                DO ne = 1, n_e
                        eq(ne)%coun=0
                        eq(ne)%f = f(eq(ne)%i,:)
                        eq(ne)%Y = eq(ne)%Y_all(eq(ne)%i)
                        eq(ne)%Yq = eq(ne)%Y - MATMUL(eq(ne)%x,eq(ne)%b)
                        eq(ne)%Yq = eq(ne)%Yq - MATMUL(eq(ne)%f,eq(ne)%a)
						eq(ne)%xmix = eq(ne)%x 
                ! sample a|b
                        IF (eq(ne)%nor%l_f) THEN
                        eq(ne)%f_f = eq(ne)%f(:,eq(ne)%nor%i_f)
                                eq(ne)%Yq = eq(ne)%Y - MATMUL(eq(ne)%f,eq(ne)%nor%zo_f*eq(ne)%a)
                                eq(ne)%Yq = eq(ne)%Yq - MATMUL(eq(ne)%xmix,eq(ne)%b)
                                CALL SAMPLE_ALPHA(eq(ne),prior_var)                     
                        END IF
                        IF (eq(ne)%nor%l_s) THEN
                                eq(ne)%f_s = eq(ne)%f(:,eq(ne)%nor%i_s)
                                eq(ne)%Yq = eq(ne)%Y - MATMUL(eq(ne)%f,eq(ne)%nor%zo_s*eq(ne)%a)
                                eq(ne)%Yq = eq(ne)%Yq - MATMUL(eq(ne)%xmix,eq(ne)%b)
                                CALL SAMPLE_ALPHA_TRUNCATED(eq(ne),prior_var)
                        END IF  
                        ! sample b|a
                        eq(ne)%Yq = eq(ne)%Y - MATMUL(eq(ne)%f,eq(ne)%a)
                        CALL SAMPLE_BETA(eq(ne))
                        ! Sample precision
                        eq(ne)%t = TAU(g1,g2,eq(ne)%Y,eq(ne)%x,eq(ne)%f,eq(ne)%b,eq(ne)%a)
                                st=0
                                WRITE(eq(ne)%file,'(40F16.8)') eq(ne)%b,eq(ne)%a,eq(ne)%t
                        IF (j==1) WRITE(4,fmt=*) 'DONE WITH LINEAR EQUATION',ne,' ITERATION 1' 
                END DO
        END IF  
        DO h = 1, Hh
        ind = 0
                m = 0
                DO k = 1, Hh
                        IF (h.ne.k)     THEN
                        m = m + 1
                        ind(m) = k
                        END IF
                END DO
                countf=0.0d0
                DO i = 1, Ii
                    IF (h.le.sample_factor(i)) THEN
                        Tf = 0.0d0
                        mf = 0.0d0                  
                        IF (n_p>0) THEN
                                DO np = 1, n_p
                                        IF (pr(np)%eq%m(i)==1.0d0) THEN
                                                IF (h==1) THEN
                                                        logc = pr(np)%d(i)==1.0d0
                                                        mn_zi = DOT_PRODUCT(pr(np)%eq%X_all(i,:),pr(np)%eq%b) + &
                                                        DOT_PRODUCT(f(i,:),pr(np)%eq%a)
                                                        pr(np)%eq%Y_all(i) = Sample_Truncated_Normal(mn_zi, &
                                                        1.0d0/pr(np)%eq%t,zero,logc)
                                                END IF
                                                Tf = Tf + pr(np)%eq%t*(pr(np)%eq%a(h)**2)
                                                IF (Hh==1) THEN
                                                        mf = mf + pr(np)%eq%t*pr(np)%eq%a(h)*(pr(np)%eq%Y_all(i) - &
                                                        DOT_PRODUCT(pr(np)%eq%x_all(i,:),pr(np)%eq%b)) 
                                                ELSE
                                                        mf = mf + pr(np)%eq%t*pr(np)%eq%a(h)*(pr(np)%eq%Y_all(i) - &
                                                        DOT_PRODUCT(pr(np)%eq%x_all(i,:),pr(np)%eq%b) &
                                                        - DOT_PRODUCT(f(i,ind),pr(np)%eq%a(ind))) 
                                                END IF
                                        END IF
                                END DO
                        END IF                                    
                        IF (n_e>0) THEN
                                DO ne = 1, n_e
                                        IF (eq(ne)%m(i)==1.0d0) THEN
                                                Tf = Tf + eq(ne)%t*(eq(ne)%a(h)**2)
                                                IF (Hh==1) THEN
                                                        mf = mf + eq(ne)%t*eq(ne)%a(h)*(eq(ne)%Y_all(i) - &
                                                        DOT_PRODUCT(eq(ne)%x_all(i,:),eq(ne)%b)) 
                                                ELSE
                                                        mf = mf + eq(ne)%t*eq(ne)%a(h)*(eq(ne)%Y_all(i) - &
                                                        DOT_PRODUCT(eq(ne)%x_all(i,:),eq(ne)%b) &
                                                        - DOT_PRODUCT(f(i,ind),eq(ne)%a(ind))) 
                                                END IF
                                        END IF
                                END DO
                        END IF                  
                        mf_c = mf/Tf
                        tf_c = Tf
                                f(i,h) = SAMPLE_F(fmix(h),mf_c,tf_c)
                    END IF
                END DO
                IF (j==1) WRITE(4,fmt=*) 'SAMPLED FACTOR',h,' ITERATION 1' 
                CALL update_Uni_normal_mixture(fmix(h),f(:,h))
                IF (j==1) WRITE(4,fmt=*) 'UPDATED MIXTURE FACTOR',h,' ITERATION 1' 
                        st=0
                        CALL write_output_Uni_nor_mixture(fmix(h))
        END DO
                OPEN(cf,file='fout.out')
                DO i = 1, Ii
                        WRITE(cf,'(40F16.8)') f(i,:)
                END DO
                CLOSE(cf)
        
        lburnin = j .ge. burnin
        if (lburnin) st=st+1
        lstep=st==step

        IF (j<100) WRITE(*,*) j
        IF (j==100) THEN
                WRITE(*,*) 'REACHED 100 ITERATIONS'
                CALL system_clock(count=fin2)
                time2=dble(fin2-in2)/dble(hz2)
                WRITE(*,*) time2
        END IF

        IF (j==burnin) WRITE(*,*) "DONE WITH BURNIN" 
        IF (j==burnin+((step*iterations)/2)) WRITE(*,*) "DONE WITH HALF",step*iterations/2," AFTER BURNIN DONE" 
        CLOSE(4)
END DO

CALL system_clock(count=finish)
time=dble(finish-initt)/dble(hz)
WRITE(*,*) time

END SUBROUTINE CODE3_LOOP
