! Rong Chen, 2021-4-22 program EM_mix use random use random2 use samplers implicit none integer, parameter :: i4=selected_int_kind(9) ! i4 range is about 10^-10 ~ 10^10. integer, parameter :: i8=selected_int_kind(15) integer, parameter :: r8=selected_real_kind(15,9) !zz real(kind=r8), parameter :: zero=0.0_r8, & !zz one=1.0_r8, & !zz two=2.0_r8, & real(kind=r8), parameter :: three=3.0_r8, & !zz four=4.0_r8, & !zz five=5.0_r8, & !zz six=6.0_r8, & !zz seven=7.0_r8, & !zz eight=8.0_r8, & !zz nine=9.0_r8, & !zz ten=10.0_r8, & tenth=0.1_r8 ! , & !zz half=.5_r8, & !zz third=one/3, & !zz sixth=one/6, & !zz pi=four*atan(one), & !zz normpdf_factor=one/sqrt(two*pi) ! 1/sqrt(2 pi) integer(kind=i8) :: irn integer(kind=i4) :: i,iter,itermax,nsub,nsub1,nsub2,mi,dim_p,kmix & ,day,hour,minute,istart,nest,istart_begin,istart_end,istart_step ! k,kmax, real(kind=r8) :: mu01,sig01,mu02,sig02,mu0V,sig0V,w01,sig0 & ,mu1,sig1,mu2,sig2,muV,sigV,w1,w2 & ,D,sig,t0,t1,t2,t01,t02,second,ttmp !! ,w1old,w1now,LL,portion,error_old,error_old1,error_old2 real(kind=r8), dimension(:), allocatable :: qt1,qt2,qt,t,V real(kind=r8), allocatable :: LL_iter(:) !! W1_iter(:),Mu1_iter(:),Mu2_iter(:),MuV_iter(:), Sigma1_iter(:),Sigma2_iter(:),SigmaV_iter(:) real(kind=r8), dimension(:,:), allocatable :: e,Y logical :: readYji ! ,result_showed,result_found,convcheck, !z real(kind=r8), dimension(:), allocatable :: x integer(kind=i4) :: m,mgauss ! n, real(kind=r8), dimension(:), allocatable :: valmusigma(:),errormusigma(:),kappa(:) character(len=8), dimension(:), allocatable :: valname(:) open ( 6, file='EM_mix.log' ) call delta_sec ( '_START' ) call delta_sec ( '# pYq_i_detail' ) ! -------------- True Values ----------------- mu01 = 0.3_r8 sig01 = 0.06_r8 mu02 = 0.6_r8 sig02 = 0.06_r8 mu0V = 20.0_r8 sig0V = 2.0_r8 w01 = 0.8_r8 sig0 = 0.1_r8 !------------------------------------------------------- irn = 90027 ! set the base random number seed call setrn(irn) itermax = 50 nsub = 100 nsub1 = int(w01*real(nsub)) nsub2 = nsub-nsub1 mi = 5 ! the size of t array. kmix = 2 dim_p = 2 readYji = .false. ! .true. m = int(10**4) ! # of samples in each iteration mgauss = int(10**3) ! set to 1000 is enough for now. mgauss need to be smaller than m the msample. !zz result_showed = .false. !zz result_found = .false. allocate(LL_iter(0:itermax)) ! '%%%%%%%%%%%% Initial Conditions %%%%%%%%%%%%' mu1 = mu01 sig1 = sig01 mu2 = mu02 sig2 = sig02 muV = mu0V sigV = sig0V w1 = w01 w2 = 1 - w1 sig = sig0 !------------------------------------- mu1 = 2 mu2 = 2 muV = 50 w1 = 0.5 sig1 = mu1/three ! five means 5 sigma ~ mean. so that values are with 5 sigma, this is very reasonable assumption. sig2 = mu2/three sigV = muV/three w2 = 1 - w1 sig = 1.0 call get_datetime write(6,'(''nsub ='',t10,i10)') nsub write(6,'(''nsub1 ='',t10,i10)') nsub1 write(6,'(''nsub2 ='',t10,i10)') nsub2 write(6,'(''----- initial condition -----'')') write(6,'(''w1 ='',t10,f12.7)') w1 write(6,'(''w2 ='',t10,f12.7)') w2 write(6,'(''mu1 ='',t10,f12.7)') mu1 write(6,'(''sig1 ='',t10,f12.7)') sig1 write(6,'(''mu2 ='',t10,f12.7)') mu2 write(6,'(''sig2 ='',t10,f12.7)') sig2 write(6,'(''muV ='',t10,f12.7)') muV write(6,'(''sigV ='',t10,f12.7)') sigV write(6,'(''sigma ='',t10,f12.7)') sig write(6,'(''read Y? '',t10,l10)') readYji write(6,*) !z call delta_sec ( 'SET-UP' ) ! insignificant time : removed ! Generate Observations data. D = 100.0_r8 qt1 = mu01 + sig01*gaussian(nsub1) qt2 = mu02 + sig02*gaussian(nsub2) qt = [qt1, qt2] V = mu0V + sig0V*gaussian(nsub) t = (/ 1.5, 2.0, 3.0, 4.0, 5.5 /) e = reshape( gaussian(nsub*mi),(/mi,nsub/) ) call samplers_init(itermax,m,mgauss,mi,nsub,dim_p,kmix & ,mu1,sig1,mu2,sig2,muV,sigV,w1,w2,sig & ,D,t) if (readYji) then call read_Yji(mi,nsub) else allocate(Y(mi,nsub)) !open(unit=11,file='Y.dat',form='formatted') do i = 1, nsub Y(:,i) = D/V(i)*exp(-t(:)*qt(i))*( 1.0_r8 + sig0*e(:,i) ) !write (11,'(5(f12.8,1x))') Y(:,i) enddo !close(11) call push_Yji(Y) deallocate(Y) endif !do i = 1, nsub ! write(6,'(5(f12.8,1x))') Y(:,i) !enddo !call step write(6,'(''Number of samples in each iteration = '',i10)') m write(6,'(''Gaussian sample (type-1) size in one loop = '',i10)') mgauss write(6,'(''Total iteration = '',i10)') itermax call delta_sec ( 'INITIALISED Yji' ) call CPU_TIME(t1) call prep (LL_iter(0)) !call preptest(LL) write (6,'(''Initial LogLike = '',t10,(g15.7))') LL_iter(0) call CPU_TIME(t0) write(6,'(''Time step0 is: '',t20,g12.5, '' seconds'')') t0-t1 !z call delta_sec ( 'PREPARE' ) ! timed in prep do iter = 1,itermax write(6,'(''Iteration ='',t10,i10)') iter call CPU_TIME(t01) call steptest(iter) ! step1(m) call CPU_TIME(t02) ttmp = t02-t01 write(6,*) '' write(6,'(''Time this loop= '',t20,f12.5, '' seconds'')') ttmp write(6,'(''Time cost= '',t20,f12.5,'' seconds'')') t02-t1 write(6,'(''Remain time = '',t20,f12.5,'' seconds'')') ttmp*(itermax-iter) call delta_sec ( 'cpu_time report' ) write(6,*) '' enddo call CPU_TIME (t2) call time_dhms (t2-t1,day,hour,minute,second) ! data analysis nest = 10 allocate(valmusigma(nest),errormusigma(nest),kappa(nest),valname(nest)) valname(1)='w1' valname(2)='w2' valname(3)='Mu1' valname(4)='Mu2' valname(5)='MuV' valname(6)='Sig1' valname(7)='Sig2' valname(8)='SigV' valname(9)='Sigma' valname(10)='LogLike' istart_begin = ceiling(tenth*real(itermax)) istart_end = ceiling(0.9_r8*real(itermax)) istart_step = ceiling(tenth*real(itermax)) !call get_musigma(1,valmusigma,errormusigma,kappa,nest) !error_old = errormusigma(10) !call get_musigma(istart_begin,valmusigma,errormusigma,kappa,nest) !error_old2 = errormusigma(10) !error_old = max(error_old1,error_old2) !write(6,*) 'errorold=', error_old !do istart = istart_begin, istart_end, istart_step ! call get_musigma(istart,valmusigma,errormusigma,kappa,nest) ! if (errormusigma(10) < error_old) then ! write (6,'(/, ''Start From Iteration '',i6,''/'',i6, '' ('',f6.2,''%) '' '' Final Result Is: '')') & ! istart,itermax,((real(itermax)-real(istart))/real(itermax))*100_r8 ! do i = 1,nest ! write (6,'(a8,'' = '',f12.5,1x,''+-'',f12.5, 10x, ''| kappa = '',f12.5)') & ! valname(i),valmusigma(i),errormusigma(i),kappa(i) ! enddo ! error_old=errormusigma(10) ! result_showed = .true. ! else ! if (result_showed) then ! write(6,*) 'Based on the current initial condition and iterations, the results might have been found.' ! write(6,*) 'However please make sure to try different intial conditions to find true global maxlkihood.' ! result_found = .true. ! exit ! else ! cycle ! endif ! endif !enddo if (.true.) then !if (.not.result_found) then ! just show all the results. call get_musigma_maxLL(LL_iter(1:itermax)) open(11,FILE='LL_iter.txt',FORM='FORMATTED') do iter = 0, itermax write(11,'(i10,1x,g15.7)') iter, LL_iter(iter) enddo close(11) do istart = istart_begin, istart_end, istart_step call get_musigma(istart,valmusigma,errormusigma,kappa,nest) write (6,'(/, "Start From Iteration ",i6,"/",i6, " (",f6.2,"%) ", " Final Result Is: ")') & istart,itermax,((real(itermax)-real(istart))/real(itermax))*100_r8 do i = 1,nest if (i /= 10) then write (6,'(a8,'' = '',f12.5,1x,''+-'',f12.5, 10x, ''| kappa = '',f12.5)') & valname(i),valmusigma(i),errormusigma(i),kappa(i) else write (6,'(a8,'' = '',g12.5,1x,''+-'',g12.5, 10x, ''| kappa = '',f12.5)') & valname(i),valmusigma(i),errormusigma(i),kappa(i) endif enddo enddo !if (.not.result_found) write(6,*) 'AS-IS results. Can try increasing samples, double itermax or change initial condition.' endif write (6,'(/,''Total time cost is: '',i10,'' days'',i10,'' hours'',i10,'' minutes'',f10.3,'' seconds'')') & day,hour,minute,second call delta_sec ( 'ANALYSED' ) call delta_sec ( '_FINISHED' ) write (6,*) ' calls to pYq_i_detail =', calls_pYq_i_detail call report_exp_chk write (*,*) 'Program end normally.' stop contains subroutine time_dhms (timeall, d,h,m,sec) ! timeall in seconds real(kind=r8) :: timeall, sec integer(kind=i4) :: d,h,m sec = timeall d = int( sec /60/60/24) ; sec = sec - d*60*60*24 h = int( sec /60/60 ) ; sec = sec - h*60*60 m = int( sec /60 ) ; sec = sec - m*60 return end subroutine time_dhms end program EM_mix subroutine delta_sec ( description ) character*(*) description ! integer*8 :: tick, rate integer*8 :: last_tick=-1 ! last tick delta_sec was called real*8 :: sec, all_sec = 0 ! logical :: do_summary = .false. integer*4, save :: nt=0, i character*30, save :: list_of_descriptions(50)=' ' integer*4, save :: num_calls(50)=0 real*8, save :: times(50)=0 ! Get ticks since last call call system_clock ( tick, rate ) if ( last_tick < 0 .or. description == '_START') then last_tick = tick all_sec = 0 do_summary = description == '_START' nt = 0 end if ! report this time interval : ignore if #.... sec = dble(tick-last_tick) / dble(rate) all_sec = all_sec + sec if ( description(1:1) /= '#' ) & write (6,11) description, sec, all_sec last_tick = tick ! if ( .not. do_summary ) return ! ! save all times for final summary report if selected do i = 1,nt if ( list_of_descriptions(i) /= description ) cycle times(i) = times(i) + sec num_calls(i) = num_calls(i) + 1 exit end do ! add if new description to list of descriptions if ( i > nt ) then if ( nt < size(times) ) nt = nt+1 list_of_descriptions(nt) = description times(nt) = times(nt) + sec num_calls(nt) = num_calls(nt) + 1 end if ! report summary times if finished if ( description == '_FINISHED') then write (6,10) nt times(nt) = sum(times(1:nt)) num_calls(nt) = sum(num_calls(1:nt)) do i = 1,nt write (6,12) i, list_of_descriptions(i), times(i), num_calls(i) times(i) = 0 num_calls(i) = 0 end do end if 10 format (/'#### Delta_Sec Summary #### ',i4// & ' Id Description Elapsed Calls') 11 format ('#### delta_sec #### ',a,t50,2f10.4) 12 format (i3,' ',a, f10.4, i9) end subroutine delta_sec