! Program to process .txt chain files produced by CosmoMC (cosmologist.info/cosmomc) ! Calculates statistics, thins, produces 1D and 2D marginalized plots and scatter plots ! and writes out useful information. ! Option to adjust priors and map parameters - write code below ! Joe Zuntz: Added pylab output and parameter scaling. !Output files are (depending on options) ! * file_root.margestats - the mean, limits and stddev from the 1D marginalized distributions ! * file_root.likestats - the best fit, and limits from the N-D confidence region ! * file_root.m - MatLab .m file to produce 1D marginalized plots ! * file_root.sm - SuperMongo .sm file to produce 1D marginalized plots ! * file_root_2D.m - MatLab .m file to produce 2D marginalized plots ! * file_root_3D.m - MatLab .m file to produce 2D sample plots colored by a third parameter ! * file_root_tri.m - MatLab .m file to produce triangle plots (1D on diagonal, 2D off-diagonal) ! * file_root_thin.txt - thinned merged versions of input files ! * file_root_single.txt - File of weight-1 single samples generated by probabilistic importance sampling ! * file_root.covmat - covariance matrix for parameters ! * file_root.corr - covariance matrix for parameters normalized to have unit variance ! * file_root.PCA - human-readable file giving details of PCA, constrained parameters, etc ! * file_root.converge - Various statistics to help assess chain convergence/sampling error ! You can easily edit the .m and .sm files produced to produce custom layouts of plots, etc. ! Data for plots are exported to the plot_data folder !This version July 2005 !March 03: fixed bug computing the limits in the .margestats file !May 03: Added support for triangle plots !Dec 03: Support for non-chain samples, auto-correlation convergence, auto_label, ! split tests to quanify sampling error on upper and lower quantiles !Jan 05: MatLab 7 options, various minor changes !Jul 05: Added limits info to .margestats output module MCSamples use settings implicit none integer, parameter :: gp = KIND(1.d0) real(gp), dimension(:,:), target, allocatable :: coldata !(col_index, row_index) integer, dimension(:), allocatable :: thin_ix integer thin_rows integer, parameter :: max_rows = 500000 integer, parameter :: max_cols = 500 integer, parameter :: max_chains = 100 integer, parameter :: max_split_tests = 4 integer chain_indices(max_chains), num_chains_used integer nrows, ncols real numsamp integer ND_cont1, ND_cont2 integer covmat_dimension integer colix(max_cols), num_vars !Parameters with non-blank labels real mean(max_cols), sddev(max_cols) real, dimension(:,:), allocatable :: corrmatrix character(LEN=120) rootname character(LEN=120) labels(max_cols) logical has_limits(max_cols), has_limits_top(max_cols),has_limits_bot(max_cols) integer num_contours real matlab_version real contours(2), contour_levels(2) real mean_mult, max_mult integer :: indep_thin = 0 integer chain_numbers(max_chains) logical isused(max_cols) logical force_twotail logical smoothing, plot_meanlikes, plot_NDcontours logical shade_meanlikes, make_single_samples integer single_thin,cust2DPlots(max_cols**2) integer ix_min(max_cols),ix_max(max_cols) real center(max_cols),width(max_cols) real, dimension(:,:), allocatable :: TheBins, bins2D, bin2Dlikes, bin2Dmax real meanlike, maxlike, maxmult logical BW,do_shading character(LEN=80) ComparePlots(20) integer Num_ComparePlots logical :: prob_label = .false. logical :: jz_have_factors real :: jz_factors(max_cols) contains subroutine AdjustPriors !Can adjust the multiplicity of each sample in coldata(1, rownum) for new priors !Be careful as this code is parameterisation dependent integer i real H0,chisq stop 'You need to write the AdjustPriors subroutine in GetDist.F90 first!' write (*,*) 'Adjusting priors' do i=0, nrows-1 !E.g. bbn prior H0 = coldata(22,i) chisq = (H0 - 72)**2/8**2 coldata(1,i) = coldata(1,i)*exp(-chisq/2) coldata(2,i) = coldata(2,i) + chisq/2 end do end subroutine AdjustPriors subroutine MapParameters(invars) real invars(1:ncols) integer i ! map parameters in invars: eg. invars(3)=invars(3)*invars(4) do i=1,ncols-2 invars(i+2)=invars(i+2)*jz_factors(i) end do end subroutine MapParameters subroutine CoolChain(cool) real, intent(in) :: cool integer i real maxL, newL write (*,*) 'Cooling chains by ', cool MaxL = minval(coldata(2,0:nrows-1)) do i=0, nrows-1 newL = coldata(2,i)*cool coldata(1,i) = coldata(1,i)*exp(-(newL - coldata(2,i)) - MaxL*(1-cool) ) coldata(2,i) = newL end do end subroutine CoolChain subroutine DeleteZeros integer i,ii ii=0 do i=0, nrows-1 if (coldata(1,i)/=0) then coldata(:,ii) = coldata(:,i) ii=ii+1 end if end do if (ii==0) stop 'Prior has removed all models!' if (ii /= nrows) write (*,*) 'Prior removed ',nrows-ii,'models' nrows = ii end subroutine DeleteZeros subroutine SortColData(bycol) !Sort coldata in order of likelihood Type(double_pointer), dimension(:), allocatable :: rowptrs integer, intent(in) :: bycol real(gp), dimension(:,:), pointer :: tmp integer i allocate(tmp(ncols,0:nrows-1)) allocate(rowptrs(0:nrows-1)) do i = 0, nrows -1 rowptrs(i)%p => coldata(1:,i) end do call QuickSortArr(rowptrs(0), 1, nrows, bycol) !Make new table using sorted pointers do i = 0, nrows-1 tmp(1:ncols,i) = rowptrs(i)%p(1:ncols) end do deallocate(rowptrs) do i=0, nrows-1 !Write this out explicitly to avoid stack overflow problems coldata(1:ncols,i) = tmp(1:ncols,i) end do deallocate(tmp) end subroutine SortColData subroutine MakeSingleSamples(single_thin) !Make file of weight-1 samples by choosing samples with probability given by their weight use Random integer i, single_thin character(LEN=20) :: fmt real maxmult Feedback = 0 call initRandom() fmt = trim(numcat('(',ncols))//'E16.7)' open(unit=50,file=trim(rootname)//'_single.txt',form='formatted',status='replace') maxmult = maxval(coldata(1,0:nrows-1)) do i= 0, nrows -1 if (ranmar() <= coldata(1,i)/maxmult/single_thin) write (50,fmt) 1.0, coldata(2:ncols,i) end do close(50) end subroutine MakeSingleSamples subroutine WriteThinData(fname,cool) character(LEN=*), intent(in) :: fname real,intent(in) :: cool character(LEN=20) :: fmt integer i real MaxL, NewL if (cool /= 1) write (*,*) 'Cooled thinned output with temp: ', cool MaxL = minval(coldata(2,0:nrows-1)) fmt = trim(numcat('(',ncols))//'E16.7)' open(unit=50,file=fname,form='formatted',status='replace') do i=0, thin_rows -1 if (cool/=1) then newL = coldata(2,thin_ix(i))*cool write (50,fmt) exp(-(newL - coldata(2,thin_ix(i))) - MaxL*(1-cool) ), newL, & coldata(3:ncols,thin_ix(i)) else write (50,fmt) 1., coldata(2:ncols,thin_ix(i)) end if end do write (*,*) 'Wrote ',thin_rows, 'thinned samples' close(50) end subroutine WriteThinData subroutine ThinData(fac, ix1,ix2) !Make thinned samples integer, intent(in) :: fac integer, intent(in), optional :: ix1,ix2 integer i, tot, nout, nend, mult if (allocated(thin_ix)) deallocate(thin_ix) allocate(thin_ix(0:nint(numsamp)/fac)) tot= 0 nout=0 i = 0 if (present(ix1)) i=ix1 nend = nrows if (present(ix2)) nend = ix2+1 mult = coldata(1,i) do while (i< nend) if (abs(nint(coldata(1,i)) - coldata(1,i)) > 1e-4) & stop 'non-integer weights in ThinData' if (mult + tot < fac) then tot = tot + mult i=i+1 if (i< nend) mult = nint(coldata(1,i)) else thin_ix(nout) = i nout= nout+1 if (mult == fac - tot) then i=i+1 if (i< nend) mult = nint(coldata(1,i)) else mult = mult - (fac -tot) end if tot = 0 end if end do thin_rows = nout close(50) end subroutine ThinData subroutine GetCovMatrix integer i, j real mean(max_cols) real scale real, dimension(:,:), allocatable :: covmatrix allocate(corrmatrix(ncols-2,ncols-2)) corrmatrix = 0 do i=1, ncols-2 if (isused(i+2)) then mean(i) = sum(coldata(1,0:nrows-1)*coldata(i+2,0:nrows-1))/numsamp end if end do do i = 1, ncols-2 if (isused(i+2)) then do j = i, ncols-2 if (isused(j+2)) then corrmatrix(i,j) = sum(coldata(1,0:nrows-1)*((coldata(2+i,0:nrows-1)-mean(i))* & (coldata(2+j,0:nrows-1)-mean(j))))/numsamp corrmatrix(j,i) = corrmatrix(i,j) end if end do end if end do if (covmat_dimension /= 0) then if (covmat_dimension > ncols -2) stop 'covmat_dimension larger than number of parameters' write (*,*) 'Writing covariance matrix for ',covmat_dimension,' parameters' allocate(covmatrix(covmat_dimension,covmat_dimension)) covmatrix=corrmatrix(1:covmat_dimension,1:covmat_dimension) call WriteSqMatrix(trim(rootname) //'.covmat',covmatrix, covmat_dimension) deallocate(covmatrix) end if do i=1, ncols-2 if (corrmatrix(i,i) > 0) then scale = sqrt(corrmatrix(i,i)) corrmatrix(i,:) = corrmatrix(i,:)/scale corrmatrix(:,i) = corrmatrix(:,i)/scale end if end do call WriteSqMatrix(trim(rootname) //'.corr',corrmatrix, ncols -2) end subroutine GetCovMatrix function MostCorrelated2D(i1,i2, direc) integer, intent(in) :: i1,i2, direc integer MostCorrelated2D !Find which parameter is most correllated with the degeneracy in ix1, ix2 integer pars(2) real mat2D(2,2), evals(2), u(2,2) real corrs(2,ncols-2) if (direc /= 0 .and. direc /= -1) stop 'Invalid 3D color parameter' pars(1)=i1 pars(2)=i2 mat2D = corrmatrix(pars,pars) u = mat2D call Diagonalize(u, evals, 2) corrs = matmul(transpose(u), corrmatrix(pars,:)) corrs(:,pars) = 0 MostCorrelated2D = MaxIndex(abs(corrs(2+direc,colix(1:num_vars)-2)), num_vars) MostCorrelated2D = colix(MostCorrelated2D) -2 end function MostCorrelated2D subroutine GetFractionIndices(fraction_indices,n) integer, intent(in) :: n integer num, fraction_indices(*), i real tot, aim tot = 0 aim=numsamp/n num = 1 fraction_indices(1) = 0 fraction_indices(n+1) = nrows do i=0, nrows-1 tot = tot + coldata(1,i) if (tot > aim) then num=num+1 fraction_indices(num) = i if (num==n) exit aim=aim+numsamp/n end if end do end subroutine GetFractionIndices function ConfidVal(ix,limfrac,upper,ix1,ix2) integer, intent(IN) :: ix real, intent(IN) :: limfrac logical, intent(IN) :: upper integer, intent(IN), optional :: ix1,ix2 real ConfidVal, samps real try, lasttry, try_t, try_b integer l,t !find upper and lower bounds l=0 t=nrows-1 if (present(ix1)) l=ix1 if (present(ix2)) t = ix2 try_b = minval(coldata(ix,l:t)) try_t = maxval(coldata(ix,l:t)) samps = sum(coldata(1,l:t)) lasttry = -1 if (upper) then do try = sum(coldata(1,l:t),mask = coldata(ix,l:t) > (try_b + try_t)/2) if (try > samps*limfrac) then try_b = (try_b+try_t)/2 else try_t = (try_b+try_t)/2 end if if (try == lasttry) exit lasttry = try end do else do try = sum(coldata(1,l:t),mask = coldata(ix,l:t) < (try_b + try_t)/2) if (try > samps*limfrac) then try_t = (try_b+try_t)/2 else try_b = (try_b+try_t)/2 end if if (try == lasttry) exit lasttry = try end do end if ConfidVal = try_t end function ConfidVal subroutine PCA(pars,n,param_map,normparam_num) !Perform principle component analysis !In other words, get eigenvectors and eigenvalues for normalized variables !with optional (log) mapping integer, intent(in) :: n, pars(n), normparam_num character(LEN=*), intent(IN) :: param_map integer i, j, normparam, locarr(1:1) character(LEN =60) fmt, tmpS real PCmean(n), sd(n), newmean(n), newsd(n) real evals(n) real corrmatrix(n,n), u(n,n) real, dimension(:,:), allocatable :: PCdata character(LEN=100) PClabs(n), div, expo logical doexp write (*,*) 'Doing PCA for ',n,' parameters' call CreateTxtFile(trim(rootname) //'.PCA', 40) write (40,*) 'PCA for parameters: ' if (normparam_num /=0) then normparam = IndexOf(normparam_num,pars,n) if (normparam ==0) stop 'Invalid PCA normalization parameter' else normparam = 0 end if allocate(PCdata(n,0:nrows-1)) PCdata(:,:) = coldata(2+pars,0:nrows-1) fmt = trim(numcat('',n))//'f8.4)' doexp =.false. do i = 1, n if (param_map(i:i)=='L') then doexp = .true. PCdata(i,:) = log(PCdata(i,:)) PClabs(i) = 'ln(' //trim(labels(pars(i)+2))//')' elseif (param_map(i:i)=='M') then doexp = .true. PCdata(i,:) = log(-PCdata(i,:)) PClabs(i) = 'ln(-' //trim(labels(pars(i)+2))//')' else PClabs(i) = trim(labels(pars(i)+2)) end if write (40,*) pars(i),':',trim(PClabs(i)) PCmean(i) = sum(coldata(1,0:nrows-1)*PCdata(i,:))/numsamp PCdata(i,:) = PCdata(i,:) - PCmean(i) sd(i) = sqrt(sum(coldata(1,0:nrows-1)*PCdata(i,:)**2)/numsamp) if (sd(i)/=0) PCdata(i,:) = PCdata(i,:)/sd(i) corrmatrix(i,i) = 1 end do write (40,*) '' write (40,*) 'Correlation matrix for reduced parameters' do i = 1, n do j = i, n corrmatrix(i,j) = sum(coldata(1,0:nrows-1)*PCdata(i,:)*PCdata(j,:))/numsamp corrmatrix(j,i) = corrmatrix(i,j) end do write (40,'(1I4,'': '','//trim(fmt)) pars(i), corrmatrix(i,:) end do u = corrmatrix call Diagonalize(u, evals, n) write (40,*) '' write (40,*) 'e-values of correlation matrix' do i = 1, n write (40,'(''PC'',1I2,'': '','//trim(fmt)) i, evals(i) end do write (40,*) '' write (40,*) 'e-vectors' do i = 1, n write (40,'(1I3,'': '','//trim(fmt)) pars(i), u(i,:) end do if (normparam /= 0) then !Set so parameter normparam has exponent 1 do i=1, n u(:,i) = u(:,i)/u(normparam,i)*sd(normparam) end do else !Normalize so main component has exponent 1 do i=1, n locarr(1:1) = maxloc(abs(u(:,i))) u(:,i) = u(:,i)/u(locarr(1),i)*sd(locarr(1)) end do end if do i = 0, nrows -1 PCdata(:,i) = matmul(transpose(u), PCdata(:,i)) if (doexp) PCdata(:,i) = exp(PCdata(:,i)) end do write (40,*) '' write (40,*) 'Principle components' do i = 1, n tmpS = trim(numcat('PC',i))//' (e-value: '//trim(RealToStr(evals(i)))//')' !ifc gives recursive IO error if you use a write within a write, !even if separate or internal files write (40,*) trim(tmpS) do j=1,n if (param_map(j:j)=='L' .or. param_map(j:j)=='M' ) then expo = RealToStr(1/sd(j)*u(j,i) ) if (param_map(j:j)=='M') then div = RealToStr( -exp(PCmean(j))) else div = RealToStr( exp(PCmean(j))) end if write(40,*) '['//trim(RealToStr(u(j,i)))//'] ('//trim(labels(pars(j)+2)) & //'/'//trim(div)//')^{'//trim(expo)//'}' else expo = RealToStr(sd(j)/u(j,i)) if (doexp) then write(40,*) '['//trim(RealToStr(u(j,i)))//'] exp(('// & trim(labels(pars(j)+2))//'-'//trim(RealToStr(PCmean(j)))// ')/'// trim(expo)//')' else write(40,*) '['//trim(RealToStr(u(j,i)))//'] ('// & trim(labels(pars(j)+2))//'-'//trim(RealToStr(PCmean(j)))// ')/'// trim(expo) end if end if end do newmean(i) = sum(coldata(1,0:nrows-1)*PCdata(i,:))/numsamp newsd(i) = sqrt(sum(coldata(1,0:nrows-1)*(PCdata(i,:)-newmean(i))**2)/numsamp) write (40,*) ' = '//trim(RealToStr(newmean(i)))// ' +- '//trim(RealToStr(newsd(i))) write (40,'('' ND limits: '',4f9.3)') minval(PCdata(i,0:ND_cont1)),maxval(PCdata(i,0:ND_cont1)), & minval(PCdata(i,0:ND_cont2)),maxval(PCdata(i,0:ND_cont2)) write (40,*) '' end do !Find out how correlated these components are with other parameters write (40,*) 'Correlations of principle components' write (40,trim(numcat('('' '',',n))//'I8)') (I, I=1, n) do i=1, n PCdata(i,:) = (PCdata(i,:) - newmean(i)) /newsd(i) end do do j=1,n write (40,'(''PC'',1I2)', advance = 'no') j do i=1,n write (40,'(1f8.3)', advance ='no') sum(coldata(1,0:nrows-1)*PCdata(i,:)* & PCdata(j,:))/numsamp end do write (40,*) '' end do do j=1, num_vars write (40,'(1I4)', advance = 'no') colix(j)-2 do i=1,n write (40,'(1f8.3)', advance ='no') sum(coldata(1,0:nrows-1)*PCdata(i,:)* & (coldata(colix(j),0:nrows-1)-mean(j))/sddev(j))/numsamp end do write (40,*) ' ('//trim(labels(colix(j)))//')' end do close(40) deallocate(PCdata) end subroutine PCA subroutine GetUsedCols integer j do j = 3, ncols isused(j) = any(coldata(j,0:nrows-1)/=coldata(j,0)) end do end subroutine GetUsedCols subroutine DoConvergeTests real chain_means(max_chains,max_cols),chain_samp(max_chains) real between_chain_var(max_cols), in_chain_var(max_cols) integer frac(max_split_tests+1), split_n, chain_start(max_chains) integer i,j,k,jj,kk,ix, maxoff real split_tests(max_split_tests) real mean(max_cols), fullmean(max_cols),fullvar(max_cols) real, parameter :: cutfrac = 0.5 real usedsamps,evals(max_cols),sc,R, maxsamp real, dimension(:,:), allocatable :: cov, meanscov integer usedvars(max_cols), num, thin_fac(max_chains), markov_thin(max_chains) real(gp) u, g2, fitted, focus, alpha, beta, probsum, tmp1 real confid integer i1,i2,i3, nburn(max_chains), hardest, endb,hardestend integer tran(2,2,2), tran2(2,2), off integer, dimension(:), allocatable :: binchain real, parameter :: epsilon = 0.001 double precision, dimension(:,:), allocatable :: corrs character(LEN=10) :: typestr integer autocorr_thin ! Get statistics for individual chains, and do split tests on the samples call CreateTxtFile(trim(rootname) //'.converge', 40) if (num_chains_used > 1) write (*,*) 'Number of chains used = ',num_chains_used chain_indices(num_chains_used+1) = nrows do i=1, num_chains_used chain_start(i) = chain_indices(i) + nint((chain_indices(i+1)-chain_indices(i))*cutfrac) chain_samp(i) = sum(coldata(1,chain_start(i):chain_indices(i+1)-1)) end do usedsamps = sum(chain_samp(1:num_chains_used)) maxsamp = maxval(chain_samp(1:num_chains_used)) num = 0 do j = 3, ncols if (isused(j)) then mean(j)=0 fullmean(j) = sum(coldata(1,0:nrows-1)*coldata(j,0:nrows-1))/numsamp if (num_chains_used> 1) then num=num+1 usedvars(num)=j do i=1, num_chains_used mean(j) = mean(j) + sum(coldata(1,chain_start(i):chain_indices(i+1)-1)*& coldata(j,chain_start(i):chain_indices(i+1)-1)) end do mean(j) = mean(j)/usedsamps end if end if end do do j = 3, ncols if (isused(j)) & fullvar(j)= sum(coldata(1,0:nrows-1)*(coldata(j,0:nrows-1)-fullmean(j))**2)/numsamp end do if (num_chains_used > 1) then write (40,*) '' write(40,*) 'Variance test convergence stats using last half chains' write (40,*) 'param var(chain mean)/mean(chain var)' write (40,*) '' do j = 3, ncols between_chain_var(j) = 0 in_chain_var(j) = 0 if (isused(j)) then if (num_chains_used > 1) then !Get stats for individual chains - the variance of the means over the mean of the variances do i=1, num_chains_used chain_means(i,j) = sum(coldata(1,chain_start(i):chain_indices(i+1)-1)* & coldata(j,chain_start(i):chain_indices(i+1)-1))/chain_samp(i) between_chain_var(j) = between_chain_var(j) + & ! chain_samp(i)/maxsamp* & !Weight for different length chains (chain_means(i,j) - mean(j))**2 in_chain_var(j) = in_chain_var(j) + & !chain_samp(i)/maxsamp *& sum(coldata(1,chain_start(i):chain_indices(i+1)-1)* & (coldata(j,chain_start(i):chain_indices(i+1)-1)-chain_means(i,j))**2) end do between_chain_var(j) = between_chain_var(j)/(num_chains_used-1) !(usedsamps/maxsamp -1) in_chain_var(j) = in_chain_var(j)/usedsamps write (40,'(1I3,f9.4," '//trim(labels(j))//'")') j-2, & between_chain_var(j) /in_chain_var(j) end if end if end do end if if (num_chains_used > 1 .and. covmat_dimension>0) then !Assess convergence in the var(mean)/mean(var) in the worst eigenvalue !c.f. Brooks and Gelman 1997 do while (usedvars(num) > covmat_dimension+2) num= num - 1 end do allocate(meanscov(num,num)) allocate(cov(num,num)) do jj=1,num j=usedvars(jj) do kk=jj, num k = usedvars(kk) meanscov(jj,kk) =0 cov(jj,kk)= 0 do i= 1, num_chains_used cov(jj,kk) = cov(jj,kk) + & !sqrt(chain_samp(jj)*chain_samp(kk))/maxsamp * & sum(coldata(1,chain_start(i):chain_indices(i+1)-1)* & (coldata(j,chain_start(i):chain_indices(i+1)-1)-chain_means(i,j))* & (coldata(k,chain_start(i):chain_indices(i+1)-1)-chain_means(i,k))) meanscov(jj,kk) = meanscov(jj,kk)+ & !sqrt(chain_samp(jj)*chain_samp(kk))/maxsamp * & (chain_means(i,j)-mean(j))*(chain_means(i,k)-mean(k)) end do meanscov(kk,jj) = meanscov(jj,kk) cov(kk,jj) = cov(jj,kk) end do end do meanscov = meanscov/(num_chains_used-1) !(usedsamps/maxsamp -1) cov = cov / usedsamps do jj=1,num sc = sqrt(cov(jj,jj)) cov(jj,:) = cov(jj,:) / sc cov(:,jj) = cov(:,jj) / sc meanscov(jj,:) = meanscov(jj,:) /sc meanscov(:,jj) = meanscov(:,jj) /sc end do call Diagonalize(meanscov, evals, num) cov = matmul(matmul(transpose(meanscov),cov),meanscov) write (40,*) '' write (40,*) 'var(mean)/mean(var) for eigenvalues of covariance of means' R = 0 do jj=1,num write (40,'(1I3,f9.4)') jj,evals(jj)/cov(jj,jj) R = max(R,evals(jj)/cov(jj,jj)) end do !R is essentially the Gelman and Rubin statistic write (*,'(" var(mean)/mean(var), 1/2 chains, worst e-value: R-1 = ",f9.4)') R deallocate(cov,meanscov) end if !Do tests for robustness under using splits of the samples !Return the rms ([change in upper/lower quantile]/[standard deviation]) !when data split into 2, 3,.. sets write (40,*) '' write (40,*) 'Split tests: rms_n([delta(upper/lower quantile)]/sd) n={2,3,4}:' write(40,*) 'i.e. mean sample splitting change in the quantiles in units of the st. dev.' write (40,*) '' do j = 3, ncols if (isused(j)) then do endb =0,1 do split_n = 2,max_split_tests call GetFractionIndices(frac,split_n) split_tests(split_n) = 0 confid = ConfidVal(j,(1-contours(num_contours))/2,endb==0,0,nrows-1) do i=1,split_n split_tests(split_n) = split_tests(split_n) + & (ConfidVal(j,(1-contours(num_contours))/2,endb==0,frac(i),frac(i+1)-1)-confid)**2 end do !i split_tests(split_n) = sqrt(split_tests(split_n)/split_n/fullvar(j)) end do if (endb==0) then typestr = 'upper' else typestr = 'lower' end if write (40,'(1I3,'//trim(IntToStr(max_split_tests-1)) // 'f9.4," ' & //trim(labels(j))//' '//trim(typestr)//'")') j-2, split_tests(2:max_split_tests) end do !endb end if end do ! Now do Raftery and Lewis method ! See http://www.stat.washington.edu/tech.reports/raftery-lewis2.ps ! Raw non-importance sampled chains only if (all(abs(coldata(1,0:nrows-1) - nint(coldata(1,0:nrows-1)))<1e-4)) then nburn = 0 do ix=1, num_chains_used thin_fac(ix) = nint(maxval(coldata(1,chain_indices(ix):chain_indices(ix+1)-1))) do j = 3, covmat_dimension+2 if (isused(j) .and. (force_twotail .or. .not. has_limits(j))) then do endb =0,1 !Get binary chain depending on whether above or below confidence value u = ConfidVal(j,(1-contours(num_contours))/2,endb==0,& chain_indices(ix),chain_indices(ix+1)-1) do !thin_fac call ThinData(thin_fac(ix),chain_indices(ix),chain_indices(ix+1)-1) if (thin_rows < 2) exit allocate(binchain(0:thin_rows-1)) where (coldata(j,thin_ix(0:thin_rows-1)) >= u) binchain = 1 elsewhere binchain = 2 endwhere tran = 0 !Estimate transitions probabilities for 2nd order process do i = 2, thin_rows-1 tran(binchain(i-2),binchain(i-1),binchain(i)) = & tran(binchain(i-2),binchain(i-1),binchain(i)) +1 end do deallocate(binchain) !Test whether 2nd order is better than Markov using BIC statistic g2 = 0 do i1=1,2 do i2=1,2 do i3=1,2 if (tran(i1,i2,i3)/=0) then fitted = dble(tran(i1,i2,1) + tran(i1,i2,2)) * & (tran(1,i2,i3) + tran(2,i2,i3)) / dble( tran(1,i2,1) + & tran(1,i2,2) + tran(2,i2,1) + tran(2,i2,2) ) focus = dble( tran(i1,i2,i3) ) g2 = g2 + log( focus / fitted ) * focus end if end do !i1 end do !i2 end do !i3 g2 = g2 * 2 if (g2 - log( dble(thin_rows-2) ) * 2 < 0) exit thin_fac(ix) = thin_fac(ix) + 1 end do !thin_fac !Get Markov transition probabilities for binary processes if (sum(tran(:,1,2))==0 .or. sum(tran(:,2,1))==0) then thin_fac(ix) = 0 goto 203 end if alpha = sum(tran(:,1,2))/dble(sum(tran(:,1,1))+sum(tran(:,1,2))) beta = sum(tran(:,2,1))/dble(sum(tran(:,2,1))+sum(tran(:,2,2))) probsum = alpha + beta tmp1 = log(probsum * epsilon / max(alpha,beta))/ log( dabs(1.0d0 - probsum) ) if (int( tmp1 + 1 ) * thin_fac(ix) > nburn(ix)) then nburn(ix) = int( tmp1 + 1 ) * thin_fac(ix) hardest = j hardestend = endb end if end do end if end do !j markov_thin(ix) = thin_fac(ix) !Get thin factor to have independent samples rather than Markov hardest = max(hardest,1) u = ConfidVal(hardest,(1-contours(num_contours))/2,hardestend==0) thin_fac(ix) = thin_fac(ix) + 1 do !thin_fac call ThinData(thin_fac(ix),chain_indices(ix),chain_indices(ix+1)-1) if (thin_rows < 2) exit allocate(binchain(0:thin_rows-1)) where (coldata(hardest,thin_ix(0:thin_rows-1)) > u) binchain = 1 elsewhere binchain = 2 endwhere tran2 = 0 !Estimate transitions probabilities for 2nd order process do i = 1, thin_rows-1 tran2(binchain(i-1),binchain(i)) = & tran2(binchain(i-1),binchain(i)) +1 end do deallocate(binchain) !Test whether independence is better than Markov using BIC statistic g2 = 0 do i1=1,2 do i2=1,2 if (tran2(i1,i2)/=0) then fitted = dble( (tran2(i1,1) + tran2(i1,2)) * (tran2(1,i2) + & tran2(2,i2)) ) / dble(thin_rows -1) focus = dble( tran2(i1,i2) ) if (fitted <= 0 .or. focus <= 0) then write (*,*) 'Raftery and Lewis estimator had problems' return end if g2 = g2 + dlog( focus / fitted ) * focus end if end do !i1 end do !i2 g2 = g2 * 2 if (g2 - log( dble(thin_rows-1) ) < 0) exit thin_fac(ix) = thin_fac(ix) + 1 end do !thin_fac 203 if (thin_rows < 2) thin_fac(ix) = 0 end do !chains write (40,*) '' write (40,*) 'Raftery&Lewis statistics' write (40,*) '' write (40,*) 'chain markov_thin indep_thin nburn' do ix = 1, num_chains_used if (thin_fac(ix)==0) then write (40,'(1I4," Not enough samples")') chain_numbers(ix) else write (40,'(1I4,3I12)') chain_numbers(ix), markov_thin(ix), & thin_fac(ix), nburn(ix) end if end do if (any(thin_fac(1:num_chains_used)==0)) then write (*,*) 'RL: Not enough samples to estimate convergence stats' else call writeS ('RL: Thin for Markov: '//& Trim(IntToStr(maxval(markov_thin(1:num_chains_used))))) indep_thin = maxval(thin_fac(1:num_chains_used)) call writeS('RL: Thin for indep samples: '// & trim(IntToStr(indep_thin))) call WriteS('RL: Estimated burn in steps: '//& trim(IntToStr(maxval(nburn(1:num_chains_used))))//' ('//& trim(IntToStr(nint(maxval(nburn(1:num_chains_used))/mean_mult)))//' rows)') end if !!Get correlation lengths write (40,*) '' write (40,*) 'Parameter auto-correlations as function of step separation' write (40,*) '' if (indep_thin ==0) then autocorr_thin = 20 elseif (indep_thin <= 30) then autocorr_thin = 5 else autocorr_thin = 5* (indep_thin/30) end if call ThinData(autocorr_thin,0,nrows-1) maxoff = min(15,thin_rows/(autocorr_thin*num_chains_used)) allocate(Corrs(ncols,maxoff)) corrs = 0 do off =1,maxoff do i=off, thin_rows-1 do j = 3, ncols if (isused(j)) & corrs(j,off) = corrs(j,off) + (coldata(j,thin_ix(i))-fullmean(j))* & (coldata(j,thin_ix(i-off)) - fullmean(j)) end do end do do j = 3, ncols if (isused(j)) & corrs(j,off) = corrs(j,off)/(thin_rows-off)/fullvar(j) end do end do write (40,'(" ",'//trim(IntToStr(maxoff)) // 'I8)') & (/(I, I=autocorr_thin, maxoff*autocorr_thin,autocorr_thin)/) do j = 3, ncols if (isused(j)) then write (40,'(1I3,'//trim(IntToStr(maxoff)) // 'f8.3," ' & //trim(labels(j))//'")') j-2, corrs(j,:) end if end do deallocate(Corrs) end if close(40) end subroutine DoConvergeTests subroutine Get2DPlotData(j,j2) integer, intent(in) :: j,j2 integer i,ix1,ix2,wx,wy real binweight, dist, norm, maxbin real try_b, try_t,try_sum, try_last character(LEN=120) :: plotfile, filename,numstr bins2D = 0 bin2Dlikes = 0 bin2Dmax = 1e30 do i = 0, nrows-1 ix1=nint((coldata(colix(j),i)-center(j))/width(j)) ix2=nint((coldata(colix(j2),i)-center(j2))/width(j2)) if (smoothing) then do wx = max(ix_min(j),ix1-2),min(ix_max(j),ix1+2) dist = (coldata(colix(j),i) - (center(j) + wx*width(j)))**2/width(j)**2 do wy = max(ix_min(j2),ix2-2),min(ix_max(j2),ix2+2) binweight = coldata(1,i)*exp( - (& dist +(coldata(colix(j2),i) - (center(j2) + wy*width(j2)))**2/width(j2)**2 )/2) bins2D(wx,wy) = bins2D(wx,wy) + binweight bin2Dlikes(wx,wy) = bin2Dlikes(wx,wy) + binweight*exp(meanlike - coldata(2,i)) end do end do else bins2D(ix1,ix2) = bins2D(ix1,ix2) + coldata(1,i) bin2Dlikes(ix1,ix2) = bin2Dlikes(ix1,ix2) + coldata(1,i)*exp(meanlike-coldata(2,i)) end if bin2Dmax(ix1,ix2) = min(bin2Dmax(ix1,ix2),real(coldata(2,i))) end do do ix1=ix_min(j),ix_max(j) do ix2 =ix_min(j2),ix_max(j2) if (bins2D(ix1,ix2) >0) bin2Dlikes(ix1,ix2) = bin2Dlikes(ix1,ix2)/bins2D(ix1,ix2) if (plot_NDcontours) then !Map maximum likelihood in each bin into a significance value from the full N-D distribution if (bin2DMax(ix1,ix2) == 1e30) then bin2DMax(ix1,ix2) = 0 else bin2DMax(ix1,ix2) = sum(coldata(1,0:nrows-1), mask = coldata(2,0:nrows-1) > bin2DMax(ix1,ix2))/numsamp end if end if end do end do if (has_limits(colix(j)) .or. has_limits(colix(j2))) then !Fix up underweighting near edges. Note this makes edge pixels noisier. do ix1 = ix_min(j), ix_max(j) do ix2 = ix_min(j2),ix_max(j2) if (ix1 ==ix_min(j) .and. has_limits_bot(colix(j))) bins2D(ix1,ix2) = bins2D(ix1,ix2)*2 if (ix1 ==ix_max(j) .and. has_limits_top(colix(j))) bins2D(ix1,ix2) = bins2D(ix1,ix2)*2 if (ix2 ==ix_min(j2).and. has_limits_bot(colix(j2))) bins2D(ix1,ix2) = bins2D(ix1,ix2)*2 if (ix2 ==ix_max(j2).and. has_limits_top(colix(j2))) bins2D(ix1,ix2) = bins2D(ix1,ix2)*2 if (ix1 ==ix_min(j)+1.and. has_limits_bot(colix(j))) bins2D(ix1,ix2) = bins2D(ix1,ix2)/0.84 if (ix1 ==ix_max(j)-1.and. has_limits_top(colix(j))) bins2D(ix1,ix2) = bins2D(ix1,ix2)/0.84 if (ix2 ==ix_min(j2)+1.and. has_limits_bot(colix(j2))) bins2D(ix1,ix2) = bins2D(ix1,ix2)/0.84 if (ix2 ==ix_max(j2)-1.and. has_limits_top(colix(j2))) bins2D(ix1,ix2) = bins2D(ix1,ix2)/0.84 end do end do end if ! Get contour containing contours(:) of the probability allocate(TheBins(ix_max(j)-ix_min(j)+1,ix_max(j2)-ix_min(j2)+1)) TheBins = bins2D(ix_min(j):ix_max(j),ix_min(j2):ix_max(j2)) norm = sum(TheBins) do ix1 = 1, num_contours try_t = maxval(TheBins) try_b = 0 try_last = -1 do try_sum = sum(TheBins,mask = TheBins < (try_b + try_t)/2) if (try_sum > (1-contours(ix1))*norm) then try_t = (try_b+try_t)/2 else try_b = (try_b+try_t)/2 end if if (try_sum == try_last) exit try_last = try_sum end do contour_levels(ix1) = (try_b+try_t)/2 end do deallocate(TheBins) plotfile = numcat(trim(numcat(trim(rootname)//'_2D_',colix(j)-2))//'_',colix(j2)-2) filename = 'plot_data/'//trim(plotfile) open(unit=49,file=filename,form='formatted',status='replace') do ix1 = ix_min(j), ix_max(j) write (49,trim(numcat('(',ix_max(j2)-ix_min(j2)+1))//'E16.7)') bins2D(ix1,ix_min(j2):ix_max(j2)) end do close(49) open(unit=49,file=trim(filename)//'_cont',form='formatted',status='replace') write(numstr,*) contour_levels(1:num_contours) if (num_contours==1) numstr = trim(numstr)//' '//trim(numstr) write (49,*) trim(numstr) close(49) if (shade_meanlikes) then open(unit=49,file=trim(filename)//'_likes',form='formatted',status='replace') maxbin = maxval(bin2Dlikes(ix_min(j):ix_max(j),ix_min(j2):ix_max(j2))) do ix1 = ix_min(j), ix_max(j) write (49,trim(numcat('(',ix_max(j2)-ix_min(j2)+1))//'E16.7)') & bin2Dlikes(ix1,ix_min(j2):ix_max(j2))/maxbin end do close(49) end if if (plot_NDcontours) then open(unit=49,file=trim(filename)//'_confid',form='formatted',status='replace') do ix1 = ix_min(j), ix_max(j) write (49,trim(numcat('(',ix_max(j2)-ix_min(j2)+1))//'E16.7)') & bin2Dmax(ix1,ix_min(j2):ix_max(j2)) end do close(49) end if end subroutine Get2DPlotData function PlotContMATLAB(aunit,aroot,j,j2, DoShade) character(LEN=*), intent(in) :: aroot integer, intent(in) :: j, j2,aunit logical, intent(in) :: DoShade character(LEN=120) fname,numstr,plotfile logical PlotContMATLAB PlotContMATLAB= .false. plotfile = numcat(trim(numcat(trim(aroot)//'_2D_',colix(j)-2))//'_',colix(j2)-2) if (.not. FileExists('plot_data/'//trim(plotfile))) return PlotContMATLAB= .true. write (aunit,'(a)') 'pts=load (fullfile(''plot_data'',''' // trim(plotfile) //'''));' fname =trim(trim(aroot)//trim(numcat('_p',colix(j2)-2))) write (aunit,'(a)') 'tmp = load (fullfile(''plot_data'',''' // trim(fname)// '.dat''));' write (aunit,*) 'x1 = tmp(:,1);' fname =trim(trim(aroot)//trim(numcat('_p',colix(j)-2))) write (aunit,'(a)') 'tmp = load (fullfile(''plot_data'',''' // trim(fname)// '.dat''));' write (aunit,*) 'x2 = tmp(:,1);' ! fmt = trim(numcat('(''x1 = ['',',ix_max(j2)-ix_min(j2)+1))//'E16.7,''];'')' ! write (aunit,fmt) (/(I, I=ix_min(j2), ix_max(j2), 1)/)*width(j2)+center(j2) ! fmt = trim(numcat('(''x2 = ['',',ix_max(j)-ix_min(j)+1))//'E16.7,''];'')' ! write (aunit,fmt) (/(I, I=ix_min(j), ix_max(j), 1)/)*width(j)+center(j) if (DoShade) then if (shade_meanlikes) then write (aunit,'(a)') 'load (fullfile(''plot_data'',''' // trim(plotfile) //'_likes''));' write (aunit,*) 'contourf(x1,x2,'//trim(plotfile)//'_likes,64);' else write (aunit,*) 'contourf(x1,x2,pts,64);' end if write (aunit,*) 'set(gca,''climmode'',''manual''); shading flat; hold on;' end if if (num_contours /= 0) then write (aunit,'(a)') 'cnt = load (fullfile(''plot_data'',''' // trim(plotfile) //'_cont''));' if (plot_NDcontours) then write (aunit,'(a)') 'load (fullfile(''plot_data'',''' // trim(plotfile) //'_confid''));' write(numstr,*) 1-contours(1:num_contours) if (num_contours==1) numstr = trim(numstr)//' '//trim(numstr) write (aunit,'(a)') 'contour(x1,x2,'//trim(plotfile)//'_confid,['//trim(numstr)//'],''r'');' end if end if end function PlotContMATLAB function PlotContPyLAB(aunit,aroot,j,j2, DoShade) character(LEN=*), intent(in) :: aroot integer, intent(in) :: j, j2,aunit logical, intent(in) :: DoShade character(LEN=120) fname,numstr,plotfile logical PlotContPyLAB integer i character(len=20) nchar PlotContPyLAB= .false. plotfile = numcat(trim(numcat(trim(aroot)//'_2D_',colix(j)-2))//'_',colix(j2)-2) if (.not. FileExists('plot_data/'//trim(plotfile))) return PlotContPyLAB= .true. write (aunit,'(a)') 'pts=load(''plot_data/'// trim(plotfile) //''')' fname =trim(trim(aroot)//trim(numcat('_p',colix(j2)-2))) write (aunit,'(a)') 'tmp = load(''plot_data/' // trim(fname)// '.dat'')' write (aunit,'(a)') 'x1 = tmp[:,0]' fname =trim(trim(aroot)//trim(numcat('_p',colix(j)-2))) write (aunit,'(a)') 'tmp = load(''plot_data/' // trim(fname)// '.dat'')' write (aunit,'(a)') 'x2 = tmp[:,0]' ! fmt = trim(numcat('(''x1 = ['',',ix_max(j2)-ix_min(j2)+1))//'E16.7,''];'')' ! write (aunit,fmt) (/(I, I=ix_min(j2), ix_max(j2), 1)/)*width(j2)+center(j2) ! fmt = trim(numcat('(''x2 = ['',',ix_max(j)-ix_min(j)+1))//'E16.7,''];'')' ! write (aunit,fmt) (/(I, I=ix_min(j), ix_max(j), 1)/)*width(j)+center(j) if (DoShade) then write (aunit,'(a)') 'mycnt=((double(range(102)))*0.01-0.01)*pts.max()' if (shade_meanlikes) then write (aunit,'(a)') 'cdat=load(''plot_data/' // trim(plotfile) //'_likes'')' write (aunit,'(a)') 'contourf(x1,x2,cdat,mycnt)' else write (aunit,'(a)') 'contourf(x1,x2,pts,mycnt)' end if end if if (num_contours /= 0) then write (aunit,'(a)') 'cnt = load(''plot_data/' // trim(plotfile) //'_cont'')' if (plot_NDcontours) then write (aunit,'(a)') 'confid=load(''plot_data/' // trim(plotfile) //'_confid'')' numstr="[" do i=1,num_contours-1 write(nchar,*) i numstr=numstr//"confid["//trim(nchar)//"]," end do write(nchar,*) num_contours numstr=numstr//'confid['//trim(nchar)//']]' write (aunit,'(a)') 'contour(x1,x2,'//trim(numstr)//'],''r'')' end if end if end function PlotContPyLAB subroutine Write2DPlotMATLAB(aunit,j,j2, DoLabelx,DoLabely) integer, intent(in) :: aunit,j,j2 logical, intent(in) :: DoLabelx,DoLabely character(LEN=120) fmt integer i if (PlotContMATLAB(aunit,rootname,j,j2,do_shading) .and. & num_contours /= 0) then write (aunit,'(a)') '[C h] = contour(x1,x2,pts,cnt,lineM{1});' write (aunit,'(a)') 'set(h,''LineWidth'',lw1);' write (aunit,*) 'hold on; axis manual; ' do i = 1, Num_ComparePlots fmt = trim(numcat('lineM{',i+1)) // '}' if (PlotContMATLAB(aunit,ComparePlots(i),j,j2,.false.)) & write (aunit,'(a)') '[C h] = contour(x1,x2,pts,cnt,'//trim(fmt)//');' write (aunit,'(a)') 'set(h,''LineWidth'',lw2);' end do end if write (aunit,*) 'hold off; set(gca,''Layer'',''top'',''FontSize'',axes_fontsize);' fmt = ''',''FontSize'',lab_fontsize);' if (DoLabelx) write (aunit,*) 'xlabel('''//trim(labels(colix(j2)))//trim(fmt) if (DoLabely) write (aunit,*) 'ylabel(''',trim(labels(colix(j)))//trim(fmt) end subroutine Write2DPlotMATLAB subroutine Write2DPlotPyLAB(aunit,j,j2, DoLabelx,DoLabely) integer, intent(in) :: aunit,j,j2 logical, intent(in) :: DoLabelx,DoLabely character(LEN=120) fmt integer i if (PlotContPyLAB(aunit,rootname,j,j2,do_shading) .and. & num_contours /= 0) then write (aunit,'(a)') 'contour(x1,x2,pts,cnt,colors=lineM[0][1],linewidth=lw1)' do i = 1, Num_ComparePlots fmt = trim(numcat('lineM[',i)) // ']' if (PlotContPyLAB(aunit,ComparePlots(i),j,j2,.false.)) & write (aunit,'(a)') 'contour(x1,x2,pts,cnt,colors='//trim(fmt)//',linewidth=lw2)' write (aunit,'(a)') 'setp(h,''LineWidth'',lw2)' end do end if ! write (aunit,*) 'hold off; set(gca,''Layer'',''top'',fontsize=axes_fontsize)' fmt = ''',fontsize=lab_fontsize)' if (DoLabelx) write (aunit,'(a)') 'xlabel('''//trim(labels(colix(j2)))//trim(fmt) if (DoLabely) write (aunit,'(a)') 'ylabel('''//trim(labels(colix(j)))//trim(fmt) write(aunit,'(a)') 'xticks(fontsize=axes_fontsize)' write(aunit,'(a)') 'yticks(fontsize=axes_fontsize)' end subroutine Write2DPlotPyLAB subroutine WriteMatLabInit(unit,sm) integer, intent(in) :: unit logical, intent(in) :: sm if (sm) then write(unit,*) 'lab_fontsize = 9; axes_fontsize = 9;' else write(unit,*) 'lab_fontsize = 12; axes_fontsize = 12;' end if if (BW) then if (plot_meanlikes) then write(unit,*) 'lineM = {''-k'',''-r'',''-b'',''-m'',''-g'',''-y''};'; write(unit,*) 'lineL = {'':k'','':r'','':b'','':m'','':g'','':y''};'; write(unit,*) 'lw1=3;lw2=1;' !Line Widths else write(unit,*) 'lineM = {''-k'',''--r'',''-.b'','':m'',''--g'',''-.y''};'; write(unit,*) 'lw1=1;lw2=1;' !Line Widths end if else write(unit,*) 'lineM = {''-k'',''-r'',''-b'',''-m'',''-g'',''-y''};'; write(unit,*) 'lineL = {'':k'','':r'','':b'','':m'','':g'','':y''};'; write(unit,*) 'lw1=1;lw2=1;' !Line Widths end if end subroutine WriteMatLabInit subroutine WritePyLabInit(unit,sm) integer, intent(in) :: unit logical, intent(in) :: sm write(unit,'(a)') 'import matplotlib' write(unit,'(a)') 'import numpy' write(unit,'(a)') 'import matplotlib.pylab' write(unit,'(a)') 'from pylab import *' write(unit,'(a)') 'from matplotlib import rc' write(unit,'(a)') 'matplotlib.use(''Agg'')' write(unit,'(a)') 'rc(''font'', family=''serif'', style=''normal'','// & 'variant=''normal'',weight=''bold'', stretch=''normal'', size=''small'')' if (sm) then write(unit,'(a)') 'lab_fontsize = 6' write(unit,'(a)') 'axes_fontsize = 6' else write(unit,'(a)') 'lab_fontsize = 8' write(unit,'(a)') 'axes_fontsize = 8' end if if (BW) then if (plot_meanlikes) then write(unit,'(a)') 'lineM = [''-k'',''-r'',''-b'',''-m'',''-g'',''-y'']' write(unit,'(a)') 'lineL = ['':k'','':r'','':b'','':m'','':g'','':y'']' write(unit,'(a)') 'lw1=3' write(unit,'(a)') 'lw2=1' !Line Widths else write(unit,'(a)') 'lineM = [''-k'',''--r'',''-.b'','':m'',''--g'',''-.y'']' write(unit,*) 'lw1=1' write(unit,'(a)')'lw2=1' !Line Widths end if else write(unit,'(a)') 'lineM = [''-k'',''-r'',''-b'',''-m'',''-g'',''-y'']' write(unit,'(a)') 'lineL = ['':k'','':r'','':b'','':m'','':g'','':y'']' write(unit,'(a)') 'lw1=1' write(unit,'(a)') 'lw2=1' !Line Widths end if end subroutine WritePyLabInit subroutine Write1DplotMatLab(aunit,j) integer, intent(in) :: aunit, j character(LEN=120) fname character(LEN=60) fmt integer ix1 fname =trim(trim(rootname)//trim(numcat('_p',colix(j)-2))) write (aunit,*) 'load (fullfile(''plot_data'',''' // trim(fname)// '.dat''));' write (aunit,*) 'pts='//trim(fname)//';' write (aunit,*) 'plot(pts[:,1),pts(:,2),lineM{1},''LineWidth'',lw1);' write (aunit,*) 'axis([-Inf,Inf,0,1.1]);axis manual;' write (aunit,*) 'set(gca,''FontSize'',axes_fontsize); hold on;' if (plot_meanlikes) then write (aunit,*) 'load (fullfile(''plot_data'',''' // trim(fname)// '.likes''));' write (aunit,*) 'pts='//trim(fname)//';' write (aunit,*) 'plot(pts(:,1),pts(:,2),lineL{1},''LineWidth'',lw1);' end if do ix1 = 1, Num_ComparePlots fname = trim(trim(ComparePlots(ix1))//trim(numcat('_p',colix(j)-2))) if (FileExists('plot_data/' // trim(fname)//'.dat')) then write (aunit,*) 'pts = load (fullfile(''plot_data'',''' // trim(fname)// '.dat''));' fmt = trim(numcat('lineM{',ix1+1)) // '}' write (aunit,*) 'plot(pts(:,1),pts(:,2),'//trim(fmt)//',''LineWidth'',lw2);' if (plot_meanlikes) then write (aunit,*) 'pts=load (fullfile(''plot_data'',''' //trim(fname)//'.likes''));' fmt = trim(numcat('lineL{',ix1+1)) // '}' write (aunit,*) 'plot(pts(:,1),pts(:,2),'//trim(fmt)//',''LineWidth'',lw2);' end if end if end do end subroutine Write1DplotMatLab subroutine Write1DplotPyLab(aunit,j) integer, intent(in) :: aunit, j character(LEN=120) fname character(LEN=60) fmt integer ix1 fname =trim(trim(rootname)//trim(numcat('_p',colix(j)-2))) write (aunit,'(a)') 'pts=load(''plot_data/' // trim(fname)// '.dat'')' write (aunit,'(a)') 'plot(pts[:,0],pts[:,1],lineM[0],linewidth=lw1)' write (aunit,'(a)') 'axis([min(pts[:,0]),max(pts[:,0]),0,1.1])' write (aunit,'(a)') 'locs,labels=xticks()' write (aunit,'(a)') 'setp(labels, fontsize=lab_fontsize)' if (plot_meanlikes) then write (aunit,'(a)') 'pts=load(''plot_data/' // trim(fname)//'.likes'')' write (aunit,'(a)') 'plot(pts[:,0],pts[:,1],lineM[0],linewidth=lw1)' end if do ix1 = 1, Num_ComparePlots fname = trim(trim(ComparePlots(ix1))//trim(numcat('_p',colix(j)-2))) if (FileExists('plot_data/' // trim(fname)//'.dat')) then write (aunit,'(a)') 'pts = load (''plot_data/' // trim(fname)// '.dat'')' fmt = trim(numcat('lineM[',ix1)) // ']' write (aunit,'(a)') 'plot(pts[:,0],pts[:,1],'//trim(fmt)//',''linewidth=lw2)' if (plot_meanlikes) then write (aunit,'(a)') 'pts=load (''plot_data/'//trim(fname)//'.likes'')' fmt = trim(numcat('lineL[',ix1)) // ']' write (aunit,'(a)') 'plot(pts[:,0],pts[:,1],'//trim(fmt)//',''linewidth=lw2)' end if end if end do end subroutine Write1DplotPyLab end module MCSamples subroutine CheckMatlabAxes(afile) integer, intent(in) :: afile write (afile,*) 'ls =get(gca,''XTick'');sz=size(ls,2);' write (afile,*) 'if(sz > 4)' write (afile,*) ' set(gca,''XTick'',ls(:,1:2:sz));' write (afile,*) 'end;' end subroutine CheckMatlabAxes subroutine CheckPylabAxes(afile) integer, intent(in) :: afile write (afile,'(a)') 'ls,loc=xticks()' write (afile,'(a)') 'sz=len(ls)' write (afile,'(a)') 'if sz > 4 :' write (afile,'(a)') ' zindex=ls.searchsorted(0)' write (afile,'(a)') ' odd=(zindex%2==1)' write (afile,'(a)') ' if odd:' write (afile,'(a)') ' newls=numpy.arange(ls[1],ls.max(),(ls[2]-ls[0]))' write (afile,'(a)') ' else:' write (afile,'(a)') ' newls=numpy.arange(ls[0],ls.max(),(ls[2]-ls[0]))' write (afile,'(a)') ' xticks(newls)' write (afile,'(a)') '' end subroutine CheckPylabAxes program GetDist use IniFile use MCSamples implicit none real bincounts(-1000:1000),binlikes(-1000:1000),binmaxlikes(-1000:1000) character(LEN=80) InputFile, numstr character(LEN=300) fname character(LEN=10000) InLine ! ,LastL,OutLine real invars(max_cols) integer ix, ix1,i,ix2,ix3, nbins, wx real ignorerows logical bad, adjust_priors real amax, amin character(LEN=120) filename, infile, in_root,out_root,sm_file character(LEN=120) matlab_col, InS1,InS2,fmt, tmpstr integer plot_row, plot_col, chain_num, first_chain,chain_ix, plot_num integer x,y,j,j2, num_2D_plots, num_cust2D_plots integer num_3D_plots character(LEN=80) contours_str, plot_3D(20), bin_limits integer plot_x, plot_y, thin_factor, outliers real limmin(max_cols),limmax(max_cols), maxbin integer plot_2D_param, j2min real try_b, try_t, binweight real cont_lines(max_cols,2,2), dist, distweight, limfrac integer bestfit_ix integer chain_exclude(max_chains), num_exclude logical map_params logical :: triangle_plot = .false. real counts real cool integer PCA_num, PCA_normparam integer PCA_params(max_cols) integer plotparams_num, plotparams(max_cols) character(LEN=max_cols) PCA_func logical Done2D(max_cols,max_cols) logical samples_are_chains, no_plots real thin_cool logical no_tests, auto_label ! Convert WMAP format file ! LastL = '' ! i=0 ! open(unit=50,file='D:\Work\cosmomc\WMAP_chains\LCDMWMAP.dat',form='formatted',status='old') ! open(unit=51,file='D:\Work\cosmomc\WMAP_chains\LCDMWMAP.txt',form='formatted',status='replace') ! do ! read (50,'(a)',end=111) InLine ! read(50,*) amax,amin ! i=i+1 ! write (outLine,'(1e15.5,a,1e15.5)') -amin, trim(InLine), amax ! if (OutLine /= LastL) then ! write (51,'(1I5,a)') i, trim(LastL) ! LastL = OutLine ! i=0 ! end if ! end do !111 close(50) ! close(51) ! stop InputFile = GetParam(1) if (InputFile == '') stop 'No parameter input file' call Ini_Open(InputFile, 1, bad, .false.) if (bad) stop 'Error opening parameter file' Ini_fail_on_not_found = .true. in_root = Ini_Read_String('file_root') rootname = ExtractFileName(in_root) chain_num = Ini_Read_Int('chain_num') ncols = Ini_Read_Int('columnnum',0) allocate(coldata(ncols,0:max_rows)) nbins = Ini_Read_Int('num_bins') ignorerows = Ini_Read_Double('ignore_rows',0.d0) adjust_priors = Ini_Read_Logical('adjust_priors',.false.) Ini_fail_on_not_found = .false. matlab_version = Ini_Read_Real('matlab_version',6.) labels(1) = 'mult' labels(2) = 'likelihood' auto_label = Ini_Read_Logical('auto_label',.false.) do ix=3, ncols write (numstr,*) ix-2 if (auto_label) then labels(ix) = trim(adjustl(numstr)) else labels(ix) = Ini_Read_String('lab'//trim(adjustl(numstr))) end if end do !JZ jz_have_factors=Ini_Read_Logical('jz_have_factors',.false.) if (jz_have_factors) then write(*,*) "Using multiplicative factors on parameters:" do ix=1,ncols-2 write (numstr,*) ix jz_factors(ix)=Ini_Read_Real('jzfac'//trim(adjustl(numstr)),1.0) write(*,*) jz_factors(ix) end do end if prob_label = Ini_Read_Logical('prob_label',.false.) triangle_plot = Ini_Read_Logical('triangle_plot',.false.) samples_are_chains = Ini_Read_Logical('samples_are_chains',.true.) no_tests = Ini_Read_Logical('no_tests',.false.) no_plots = Ini_Read_Logical('no_plots',.false.) thin_factor = Ini_Read_Int('thin_factor',0) thin_cool = Ini_read_Real('thin_cool',1.) first_chain = Ini_Read_Int('first_chain',1) make_single_samples = Ini_Read_logical('make_single_samples',.false.) single_thin = Ini_Read_Int('single_thin',1) cool = Ini_Read_Real('cool',1.) Num_ComparePlots = Ini_Read_Int('compare_num',0) do ix = 1, Num_ComparePlots ComparePlots(ix) = ExtractFileName(Ini_Read_String(numcat('compare',ix))) end do has_limits_top = .false. has_limits_bot = .false. bin_limits = Ini_Read_String('all_limits') do ix=3, ncols write (numstr,*) ix-2 if (bin_limits /= '') then InLine = bin_limits else InLine = Ini_Read_String('limits'//trim(adjustl(numstr))) end if if (InLine /= '') then read (InLine,*) InS1, InS2 if (trim(adjustl(InS1)) /= 'N') then has_limits_bot(ix) = .true. read(InS1,*) limmin(ix) end if if (trim(adjustl(InS2)) /= 'N') then has_limits_top(ix) = .true. read(InS2,*) limmax(ix) end if end if end do has_limits = has_limits_top .or. has_limits_bot plot_2D_param = Ini_Read_Int('plot_2D_param',0) if (plot_2D_param /= 0) then plot_2D_param = plot_2D_param + 2 num_cust2D_plots = 0 else !Use custom array of specific plots num_cust2D_plots = Ini_Read_Int('plot_2D_num',0) do ix = 1, num_cust2D_plots InLine = Ini_Read_String(numcat('plot',ix)) read (InLine, *) ix1,ix2 cust2DPLots(ix) = ix1+2 + (ix2+2)*100 end do end if InS1 =Ini_Read_String('exclude_chain') num_exclude = 0 chain_exclude = 0 read (InS1, *, end =20) chain_exclude 20 do while (chain_exclude(num_exclude+1)/=0) num_exclude = num_exclude + 1 end do map_params = Ini_Read_Logical('map_params',.false.) if (map_params) & write (*,*) 'WARNING: Mapping params - .covmat file is new params.' shade_meanlikes = Ini_Read_Logical('shade_meanlikes',.false.) matlab_col = Ini_Read_String('matlab_colscheme') smoothing = Ini_Read_Logical('smoothing',.true.) out_root = Ini_Read_String('out_root') if (out_root /= '') then rootname = out_root write (*,*) 'producing files with root '//trim(out_root) end if sm_file = trim(rootname) // '.sm' contours(1) = Ini_Read_Real('contour1',0.) contours(2) = Ini_Read_Real('contour2',0.) contours_str = trim(Ini_Read_String('contour1')) // ' ; '// & trim(Ini_Read_String('contour2')) num_contours = count(contours/= 0) if (contours(1)==0) contours(1) = contours(2) force_twotail = Ini_Read_Logical('force_twotail',.false.) if (force_twotail) write (*,*) 'Computing two tail limits' covmat_dimension = Ini_Read_Int('cov_matrix_dimension',0) if (covmat_dimension == -1) covmat_dimension = ncols-2 plot_meanlikes = Ini_Read_Logical('plot_meanlikes',.false.) plot_NDcontours = Ini_Read_Logical('plot_NDcontours',.false.) PCA_num = Ini_Read_Int('PCA_num',0) if (PCA_num /= 0) then if (PCA_num <2) stop 'Can only do PCA for 2 or more parameters' InLine = Ini_Read_String('PCA_params') PCA_func =Ini_Read_String('PCA_func') !Characters representing functional mapping if (PCA_func == '') PCA_func(1:PCA_num) = 'N' !no mapping if (InLine == 'ALL' .or. InLine == 'all') then PCA_params(1:PCA_num) = (/ (i, i=1,PCA_num)/) else read (InLine,*) PCA_params(1:PCA_num) end if PCA_NormParam = Ini_Read_Int('PCA_normparam',0) end if plotparams_num = Ini_Read_Int('plotparams_num',0) if (plotparams_num /= 0) then InLine = Ini_Read_String('plotparams') read(InLine,*) plotparams(1:plotparams_num) end if num_3D_plots = Ini_Read_Int('num_3D_plots',0) do ix =1, num_3D_plots plot_3D(ix) = Ini_Read_String(numcat('3D_plot',ix)) end do BW = Ini_Read_Logical('B&W',.false.) do_shading = Ini_Read_Logical('do_shading',.true.) call Ini_Close allocate(bins2D(-1000:1000,-1000:1000),bin2Dlikes(-1000:1000,-1000:1000),bin2Dmax(-1000:1000,-1000:1000)) !Read in the chains nrows = 0 num_chains_used =0 do chain_ix = first_chain, first_chain-1 + max(1,chain_num) do ix = 1, num_exclude if (chain_exclude(ix) == chain_ix) goto 30 end do if (chain_num == 0) then infile = trim(in_root) // '.txt' else write (numstr,*) chain_ix infile = trim(in_root) //'_'//trim(adjustl(numstr))// '.txt' end if write (*,*) 'reading ' // trim(infile) open(unit=50,file=infile,form='formatted',status='old', err=28) if (ignorerows >=1) then do ix = 1, nint(ignorerows) read (50,'(a)',end=1) InLine end do end if num_chains_used = num_chains_used + 1 if (num_chains_used > max_chains) stop 'Increase max_chains in GetDist' chain_indices(num_chains_used) = nrows chain_numbers(num_chains_used) = chain_ix do read (50,'(a)',end=1) InLine if (samples_are_chains) then read(InLine,*,end=1, err=2) invars(1:ncols) else read(InLine,*,end=1, err=2) invars(3:ncols) invars(1) = 1 invars(2) = 1 end if if (map_params) call MapParameters(invars) coldata(1:ncols, nrows) = invars(1:ncols) nrows = nrows + 1 if (nrows > max_rows) stop 'need to increase max_rows' end do 2 write (*,*) 'error reading line ', nrows + int(ignorerows), ' - skipping rest of file' 1 close(50) if (ignorerows<1 .and. ignorerows/=0) then i = chain_indices(num_chains_used) j = nint((nrows-i-1)*ignorerows) do ix = i,nrows-j-1 coldata(:,ix) = coldata(:,ix+j) end do nrows = nrows - j end if goto 30 28 write (*,'(" chain ",1I4," missing")') chain_ix 30 end do if (nrows == 0) stop 'No un-ignored rows! (check number of chains/burn in)' if (cool /= 1) call CoolChain(cool) !Adjust weights if requested if (adjust_priors) then call AdjustPriors end if call GetUsedCols !See which parameters are fixed mean_mult = sum(coldata(1,0:nrows-1))/nrows max_mult = (mean_mult*nrows)/min(nrows/2,500) outliers = count(coldata(1,0:nrows-1) > max_mult) if (outliers /=0) write (*,*) 'outlier fraction ', real(outliers)/nrows maxmult = maxval(coldata(1,0:nrows-1)) numsamp = sum(coldata(1,0:nrows-1)) if (.not. no_tests) call DoConvergeTests if (adjust_priors) call DeleteZeros write (*,*) 'mean input multiplicity = ',mean_mult !Output thinned data if requested !Must do this with unsorted output if (thin_factor /= 0) then call ThinData(thin_factor) call WriteThinData(trim(rootname)//'_thin.txt',thin_cool) end if !Produce file of weight-1 samples if requested if (num_3D_plots/=0 .and. .not. make_single_samples ) then make_single_samples = .true. single_thin = max(1,nint(numsamp/maxmult)/2000) end if if (make_single_samples) call MakeSingleSamples(single_thin) !Sort data in order of likelihood of points call SortColData(2) numsamp = sum(coldata(1,0:nrows-1)) !Get ND confidence region (index into sorted coldata) counts = 0 ND_cont1=-1; ND_cont2=-1 do j=0, nrows -1 counts = counts + coldata(1,j) if (counts > numsamp*contours(1) .and. ND_cont1==-1) then ND_cont1 = j if (ND_cont2 /= -1) exit end if if (counts > numsamp*contours(2) .and. ND_cont2==-1) then ND_cont2 = j end if end do !Only use variables whose label's are not empty (and in list of plotparams if plotparams_num /= 0) num_vars = 0 do ix = 3,ncols if (labels(ix) /= '' .and. isused(ix)) then if (plotparams_num == 0 .or. any(plotparams(1:plotparams_num)==ix-2)) then num_vars = num_vars + 1 colix(num_vars) = ix mean(num_vars) = sum(coldata(1,0:nrows-1)*coldata(ix,0:nrows-1))/numsamp sddev(num_vars) = sqrt(sum(coldata(1,0:nrows-1)*(coldata(ix,0:nrows-1) & -mean(num_vars))**2)/numsamp) end if end if end do triangle_plot = triangle_plot .and. (num_vars > 1) write (*,*) 'using ',nrows,' rows, processing ',num_vars,' parameters' if (indep_thin/=0) then write (*,*) 'Approx indep samples: ', nint(numsamp/indep_thin) else write (*,*) 'effective number of samples (assuming indep): ', nint(numsamp/maxmult) end if !Get covariance matrix and correlation matrix call GetCovMatrix if (PCA_num>0) call PCA(PCA_params,PCA_num,PCA_func, PCA_NormParam) !Find best fit, and mean likelihood bestfit_ix = 0 !Since we have sorted the lines maxlike = coldata(2,bestfit_ix) write (*,*) 'Best fit -Ln(like) = ', maxlike if (coldata(2,nrows-1) - maxlike < 30) then meanlike = log(sum(exp((coldata(2,0:nrows-1) -maxlike))*coldata(1,0:nrows-1)) & / numsamp) + maxlike write (*,*) 'Ln(mean 1/like) = ', meanlike end if meanlike = sum(coldata(2,0:nrows-1)*coldata(1,0:nrows-1)) / numsamp write (*,*) 'mean(-Ln(like)) = ', meanlike meanlike = -log(sum(exp(-(coldata(2,0:nrows-1) -maxlike))*coldata(1,0:nrows-1)) & / numsamp) + maxlike write (*,*) '-Ln(mean like) = ', meanlike !Get marge stats do j = 1,num_vars ix = colix(j) do ix1 = 1, num_contours limfrac = (1-contours(ix1)) if (.not. has_limits(ix) .or. force_twotail) then limfrac = limfrac/2 !two tail end if cont_lines(j,2,ix1) = ConfidVal(ix,limfrac,.true.) cont_lines(j,1,ix1) = ConfidVal(ix,limfrac,.false.) end do !contour lines end do !write out stats open(unit=50,file=trim(rootname)//'.margestats',form='formatted',status='replace') write(50,'(a)') 'param mean sddev lower1 upper1 lower2 upper2' open(unit=51,file='plot_data/'//trim(rootname)//'_params',form='formatted',status='replace') do j=1, num_vars write(50,'(1I5,6E14.6," '//trim(labels(colix(j)))//'")') colix(j)-2, mean(j), sddev(j), & cont_lines(j,1:2,1),cont_lines(j,1:2,2) write (51,*) colix(j)-2, trim(labels(colix(j))) end do close(51) write (50,*) '' write (50,'(a)') 'Limits are: ' // trim(contours_str) if (.not. force_twotail) then do j=1, num_vars if (has_limits(colix(j))) then write(50,'(1I5," one tail; '//trim(labels(colix(j)))//'")') colix(j)-2 else write(50,'(1I5," two tail; '//trim(labels(colix(j)))//'")') colix(j)-2 end if end do else write (50,*) 'All limits are two tail' end if close(50) if (no_plots) stop !Output files for 1D plots plot_col = nint(sqrt(num_vars/1.4)) plot_row = (num_vars +plot_col-1)/plot_col !So we don't forget to run sm, delete current .ps file call DeleteFile(trim(rootname)//'.ps') open(unit=51,file=trim(rootname) // '.m',form='formatted',status='replace') !MatLab file for 1D plots call WriteMatLabInit(51,num_vars>3) open(unit=82,file=trim(rootname) // '.py', form='formatted',status='replace') call WritePyLabInit(82,num_vars>3) open(unit=50,file=sm_file,form='formatted',status='replace') !SuperMongo .sm file for 1D plots write (50,*) 'DEFINE TeX_strings 1' write (50,*) 'device postfile '//trim(rootname)//'.ps' write (50,*) 'DEFINE default_font rm' if (num_vars > 6) then write (50,*) 'expand 0.8' write (50, *)'define x_gutter 0.6' write (50,*) 'define y_gutter 0.8' else write (50,*) 'expand 1.1' endif write (50,*) 'location 2500 31000 3500 31000' if (BW) then write (50,*) 'lweight 4' else write (50,*) 'lweight 3' end if if (triangle_plot) then open(unit=52,file=trim(rootname) // '_tri.m',form='formatted',status='replace') open(unit=83,file=trim(rootname) // '_tri.py',form='formatted',status='replace') call WriteMatLabInit(52,num_vars>4) call WritePyLabInit(83,num_vars>4) end if cont_lines = 0 !Do 1D bins do j = 1,num_vars ix = colix(j) bincounts = 0 binlikes = 0 binmaxlikes = 0 if (has_limits_bot(ix)) then amin = limmin(ix) else amin = ConfidVal(ix,0.001,.false.) end if if (has_limits_top(ix)) then amax = limmax(ix) else amax = ConfidVal(ix,0.001,.true.) end if width(j) = (amax-amin)/(nbins+1) if (width(j)==0) cycle if (amin <= 0) then center(j) = amax else center(j) = amin end if ix_min(j) = nint((amin - center(j))/width(j)) ix_max(j) = nint((amax - center(j))/width(j)) if (smoothing .and. .not. has_limits_bot(ix) .and. amin > 0) ix_min(j) = ix_min(j)-1 if (smoothing .and. .not. has_limits_top(ix)) ix_max(j) = ix_max(j)+1 do i = 0, nrows-1 ix2=nint((coldata(ix,i)-center(j))/width(j)) if (smoothing .and. ix_min(j) /= ix_max(j)) then do wx = max(ix_min(j),ix2-2),min(ix_max(j),ix2+2) dist = coldata(colix(j),i) - (center(j) + wx*width(j)) distweight = exp(- dist**2/width(j)**2/2) binweight = coldata(1,i)*distweight bincounts(wx) = bincounts(wx) + binweight binlikes(wx) = binlikes(wx) + binweight*exp(meanlike-coldata(2,i)) binmaxlikes(wx) = max(binmaxlikes(wx), & real(exp(maxlike-coldata(2,i))*distweight)) end do else bincounts(ix2) = bincounts(ix2) + coldata(1,i) binlikes(ix2) = binlikes(ix2) + coldata(1,i)*exp(meanlike-coldata(2,i)) binmaxlikes(ix2) = max(binmaxlikes(ix2), real(exp(maxlike-coldata(2,i)))) end if end do do i=ix_min(j),ix_max(j) if (bincounts(i) >0) binlikes(i) = binlikes(i)/bincounts(i) end do if (smoothing .and. ix_min(j) /= ix_max(j)) then !account for underweighting near edges if (has_limits_bot(ix)) then bincounts(ix_min(j)) = bincounts(ix_min(j))*2 bincounts(ix_min(j)+1) = bincounts(ix_min(j)+1)/0.84 end if if (has_limits_top(ix)) then bincounts(ix_max(j)) = bincounts(ix_max(j))*2 bincounts(ix_max(j)-1) = bincounts(ix_max(j)-1)/0.84 end if end if fname =trim(trim(rootname)//trim(numcat('_p',colix(j)-2))) filename = 'plot_data/'// trim(fname) open(unit=49,file=trim(filename)//'.dat',form='formatted',status='replace') maxbin = maxval(bincounts(ix_min(j):ix_max(j))) if (maxbin==0) then fname = IntToStr(colix(j)-2) write (*,*) 'no samples in bin, param:'//trim(fname) stop end if do i = ix_min(j), ix_max(j) write (49,'(3E16.7)') center(j) + i*width(j), bincounts(i)/maxbin, binmaxlikes(i) end do if (ix_min(j) == ix_max(j)) write (49,'(1E16.7,'' 0'')') center(j) + ix_min*width(j) close(49) !Detect ends which don't go to zero if (.not. force_twotail) then if (bincounts(ix_max(j)) > maxbin*.15) cont_lines(j,2,:) = amax if (bincounts(ix_min(j)) > maxbin*.15) cont_lines(j,1,:) = amin end if if (plot_meanlikes) then open(unit=49,file=trim(filename)//'.likes',form='formatted',status='replace') maxbin = maxval(binlikes(ix_min(j):ix_max(j))) do i = ix_min(j), ix_max(j) write (49,'(2E16.7)') center(j) + i*width(j), binlikes(i)/maxbin end do close(49) end if plot_x = mod(j+plot_col-1,plot_col)+1 plot_y = plot_row - (j+plot_col-1)/plot_col + 1 !jz write (51,*) 'subplot(',plot_row,',',plot_col,',',j,')' call Write1DplotMatLab(51,j); write(82,'(''subplot('',1I5,'','',1I5,'','',1I5,'')'')') plot_row,plot_col,j call Write1DplotPyLab(82,j) call CheckPylabAxes(82) !JZ if (matlab_version < 7) then if (plot_col*plot_row> 4) call CheckMatlabAxes(51) write (51,*) ' p = get(gca,''Position'');' if (plot_col*plot_row > 8) then write (51,*) 'set(gca,''Position'',[p(1) (p(2)+p(4)/4) p(3) p(4)*4/5]);' else write (51,*) 'set(gca,''Position'',[p(1) (p(2)+p(4)/8) p(3) p(4)*8/9]);' end if end if if (prob_label) then write (51,*) 'ylabel(''Probability'')' write (82,'(a)') 'ylabel(''Probability'')' end if write(51,*) 'xlabel('''// trim(labels(ix))//''',''FontSize'',lab_fontsize);' write(82,'(a)') 'xlabel('''// trim(labels(ix))//''',fontsize=lab_fontsize)' write(82,'(a)') 'xticks(fontsize=axes_fontsize)' write (51,*) 'set(gca,''ytick'',[]);hold off;' write(82,'(a)') 'yticks([])' if (triangle_plot) then write (52,*) 'subplot(',num_vars,',',num_vars,',',(j-1)*num_vars+j,')' write(83,'(''subplot('',1I5,'','',1I5,'','',1I5,'')'')') num_vars,num_vars,(j-1)*num_vars+j call Write1DplotMatLab(52,j) call Write1DPlotPylab(83,j) if (num_vars > 2) call CheckMatlabAxes(52) if (num_vars > 2) call CheckPylabAxes(52) if (prob_label) then write (52,*) 'ylabel(''Probability'')' write (83,*) 'ylabel(''Probability'')' end if if (j==num_vars) then write(52,*) 'xlabel('''// trim(labels(ix))//''',''FontSize'',lab_fontsize);' write(83,'(a)') 'xlabel('''// trim(labels(ix))//''',fontsize=lab_fontsize)' end if write (52,*) 'set(gca,''ytick'',[]);hold off;' write(83,'(a)') 'yticks([])' end if write (50,'(''window'',4I4)') plot_col,plot_row, plot_x, plot_y write (50,*) 'data "' // trim(filename) //'.dat"' write (50,*) 'read { x 1 y 2 mx 3}' if (ix_min(j) == ix_max(j)) then write (50,'(''limits '',2E16.7,'' y'')') center(j) - abs(center(j))/2-0.1, & center(j) + abs(center(j))/2+0.1 else if (has_limits(ix)) then write (50,'(''limits '',2E16.7,'' 0 1.1'')') amin, amax else write (50,*) 'limits x 0 1.1' end if end if write (50,*) 'box 1 3 3 3' write (50,*) 'connect x y' ! write (50,*) 'connect x mx ' !uncomment to plot maximum likelihood values for each bin if (plot_meanlikes) then write (50,*) 'data "' // trim(filename)//'.likes"' write (50,*) 'read { x 1 y 2 }' write (50,*) 'ltype 1' write (50,*) 'connect x y' end if do ix1 = 1, Num_ComparePlots fname = trim(trim(ComparePlots(ix1))//trim(numcat('_p',colix(j)-2))) if (FileExists('plot_data/' // trim(fname)//'.dat')) then write (50,*) 'data "plot_data/' // trim(fname) //'.dat"' write (50,*) 'read { x 1 y 2 mx 3}' write (50,*) 'ltype 0' if (BW) then write (50,*) 'lweight 2' if (ix1>1) write (50,*) 'ltype 2' else write (50,*) 'ctype ',ix1+2 end if write (50,*) 'connect x y' if (plot_meanlikes) then write (50,*) 'data "plot_data/'//trim(fname)//'.likes"' write (50,*) 'read { x 1 y 2}' write (50,*) 'ltype 1' write (50,*) 'connect x y' end if end if end do if (Num_ComparePlots/=0) write (50,*) 'ctype black' if (Num_ComparePlots/=0 .and. BW) write (50,*) 'lweight 4' ! write (50,*) 'ltype 1' ! do ix1 = 1, num_contours ! do i2 =1,2 ! write (numstr,*) cont_lines(j,i2,ix1) ! write (50,*) 'set x = {'//trim(numstr)//' '//trim(numstr)//'}' ! write (50,*) 'set y = { 0 1.1}' !! write (50,*) 'connect x y' !uncomment this if you want dooted lines on 1D plots for condifence regions !Note 1D confidence regions are computed from full samples using different definition to !2D confidence regions ! end do ! end do write (50,*) 'ltype 0' write (50,*) 'xlabel '// trim(labels(ix)) write (50,*) write (50,*) '#Next...' end do write (50,*) 'quit' close(50) write (51,*) 'set(gcf, ''PaperUnits'',''inches'');' ! write(82,*) 'setp(gcf(),''PaperUnits'',''inches'')' !JZ think this is unnecessary if (plot_row < plot_col .and. plot_row*plot_col>9) then write (51,*) 'set(gcf, ''PaperPosition'',[ 0 0 8 10]);'; else write (51,*) 'set(gcf, ''PaperPosition'',[ 0 0 8 8]);'; end if write (51,*) 'print -dpsc2 '//trim(rootname)//'.ps;' write(82,'(a)') 'savefig('''//trim(rootname)//'.ps'')' close(51) close(82) !do 2D bins and produce matlab file if (plot_2D_param == 0 .and. num_cust2D_plots==0) then !In this case output the most correlated variable combinations write (*,*) 'doing 2D plots for most correlated variables' try_t = 1e5 x=0;y=0; num_cust2D_plots = 12 do j = 1, num_cust2D_plots try_b = -1e5 do ix1 =1 ,num_vars do ix2 = ix1+1,num_vars if (abs(corrmatrix(colix(ix1)-2,colix(ix2)-2)) < try_t .and. & abs(corrmatrix(colix(ix1)-2,colix(ix2)-2)) > try_b) then try_b = abs(corrmatrix(colix(ix1)-2,colix(ix2)-2)) x = ix1; y = ix2; end if end do end do if (try_b == -1e5) then num_cust2D_plots = j-1 exit end if try_t = try_b cust2Dplots(j) = colix(x) + colix(y)*100 end do end if if (num_cust2D_plots == 0) then num_2D_plots = 0 do j= 1, num_vars if (ix_min(j) /= ix_max(j)) then do j2 = j+1, num_vars if (ix_min(j2) /= ix_max(j2)) then if (plot_2D_param==0 .or. plot_2D_param == colix(j) .or. plot_2D_param == colix(j2)) & num_2D_plots = num_2D_plots + 1 end if end do end if end do else num_2D_plots = num_cust2D_plots end if done2D= .false. if (num_2D_plots > 0) then write (*,*) 'Producing ',num_2D_plots,' 2D plots' filename = trim(rootname)//'_2D.m' open(unit=50,file=filename,form='formatted',status='replace') filename = trim(rootname)//'_2D.py' open(unit=81,file=filename,form='formatted',status='replace') call WriteMatLabInit(50,num_2D_plots >=7) call WritePyLabInit(81,num_2D_plots >=7) plot_col = nint(sqrt(num_2D_plots/1.4)) plot_row = (num_2D_plots +plot_col-1)/plot_col plot_num = 0 do j= 1, num_vars if (ix_min(j) /= ix_max(j)) then if (plot_2D_param/=0 .or. num_cust2D_plots /= 0) then if (colix(j) == plot_2D_param) cycle j2min = 1 else j2min = j+1 end if do j2 = j2min, num_vars if (ix_min(j2) /= ix_max(j2)) then if (plot_2D_param/=0 .and. colix(j2) /= plot_2D_param) cycle if (num_cust2D_plots /= 0 .and. & count(cust2Dplots(1:num_cust2D_plots) == colix(j)*100 + colix(j2))==0) cycle plot_num = plot_num + 1 done2D(j,j2) = .true. call Get2DPlotData(j,j2) write (50,'(''subplot('',1I5,'','',1I5,'','',1I5,'');'')') plot_row,plot_col,plot_num write (81,'(''subplot('',1I5,'','',1I5,'','',1I5,'');'')') plot_row,plot_col,plot_num call Write2DPlotMATLAB(50,j,j2,.true.,.true.) call Write2DPlotPyLAB(81,j,j2,.true.,.true.) if (plot_row*plot_col > 4) call CheckMatlabAxes(50) if (plot_row*plot_col > 4) call CheckPylabAxes(81) end if end do end if end do write (50,*) 'set(gcf, ''PaperUnits'',''inches'');' if (matlab_col/='') write (50,*) trim(matlab_col) if (num_2D_plots < 5) then write (50,*) 'set(gcf, ''PaperPosition'',[ 0 0 6 6]);'; else write (50,*) 'set(gcf, ''PaperPosition'',[ 0 0 8 10]);'; end if write (50,*) 'print -dpsc2 '//trim(rootname)//'_2D.ps;' write (81,'(a)') 'savefig('''//trim(rootname)//'_2d.ps'')' close(50) end if if (triangle_plot) then !Add the off-diagonal 2D plots do j = 1, num_vars do j2 = j+1, num_vars if (.not. Done2D(j2,j)) call Get2DPlotData(j2,j) write (52,'(''subplot('',1I5,'','',1I5,'','',1I5,'');'')') num_vars, num_vars, & (j2-1)*num_vars + j write (83,'(''subplot('',1I5,'','',1I5,'','',1I5,'');'')') num_vars, num_vars, & (j2-1)*num_vars + j call Write2DPlotMATLAB(52,j2,j,j2==num_vars, j==1) call Write2DPlotPylab(83,j2,j,j2==num_vars, j==1) if (num_vars > 2) call CheckMatlabAxes(52) if (num_vars > 2) call CheckPylabAxes(52) end do end do write (52,*) 'set(gcf, ''PaperUnits'',''inches'');' write (52,*) 'set(gcf, ''PaperPosition'',[ 0 0 8 8]);'; write (52,*) 'print -dpsc2 '//trim(rootname)//'_tri.ps;' write(83,'(a)') 'savefig('''//trim(rootname)//'_tri.ps'')' close(52) close(83) end if !Do 3D plots (i.e. 2D scatter plots with coloured points) if (num_3D_plots /=0) then write (*,*) 'producing ',num_3D_plots, '2D colored scatter plots' filename = trim(rootname)//'_3D.m' open(unit=50,file=filename,form='formatted',status='replace') filename = trim(rootname)//'_3D.py' !JZ open(unit=80,file=filename,form='formatted',status='replace') !JZ write(80,'(a)') 'import matplotlib' write(80,'(a)') 'import numpy' write(80,'(a)') 'import matplotlib.pylab' write(80,'(a)') 'from pylab import *' write(80,'(a)') 'from matplotlib import rc' write(80,'(a)') 'matplotlib.use(''Agg'')' write(80,'(a)') 'rc(''font'', family=''serif'', style=''normal'','// & 'variant=''normal'',weight=''bold'', stretch=''normal'', size=''small'')' write (50,*) 'colormap(''jet'');' if (num_3D_plots ==1 ) then write(50,*) 'lab_fontsize = 12; axes_fontsize = 12;' write (80,'(a)') 'lab_fontsize = 12; axes_fontsize = 12;'!JZ else write(50,*) 'lab_fontsize = 9; axes_fontsize = 9;' write (80,'(a)') 'lab_fontsize = 9; axes_fontsize = 9;' !JZ end if if (mod(num_3D_plots,2)==0 .and. num_3D_plots < 11) then plot_col = num_3D_plots/2 else plot_col = nint(sqrt(1.*num_3D_plots)) end if plot_row = (num_3D_plots +plot_col-1)/plot_col write(50,*) 'pts = load('''//trim(rootname)//'_single.txt'');' write (80,'(a)') 'pts = load(''' // trim(rootname)// '_single.txt'')' !JZ do j=1, num_3D_plots read(plot_3D(j),*) ix1, ix2, ix3 !x, y, color if (ix3<1) ix3 = MostCorrelated2D(ix1,ix2,ix3) write (50,'(''subplot('',1I5,'','',1I5,'','',1I5,'');'')') plot_row,plot_col,j write (80,'(''subplot('',1I5,'','',1I5,'','',1I5,'');'')') plot_row,plot_col,j !JZ write (50,*) '%Do params ',ix1,ix2,ix3 write (50,*) 'scatter(pts(:,',2+ix1,'),pts(:,',2+ix2,'),3,pts(:,',2+ix3,'));' write(80,'(''scatter(pts[:,'',1I5,''],pts[:,'',1I5,''],c=pts[:,'',1I5,''],cmap=cm.jet)'')'), 2+ix1,2+ix2,2+ix3 ! write(80,'(a)') 'scatter(pts[:,'//trim(2+ix1)//'],pts[:,'//trim(2+ix2)//')],3,c=pts[:,'//trim(2+ix3)//'],cmap=cm.jet)' fmt = ''',''FontSize'',lab_fontsize);' write (50,*) 'xlabel('''//trim(labels(ix1+2))//trim(fmt) write (50,*) 'ylabel('''//trim(labels(ix2+2))//trim(fmt) write (50,*) 'set(gca,''FontSize'',axes_fontsize); pa = gca;' write (50,*) 'hbar = colorbar(''horiz'');axes(hbar);' write (50,*) 'xlabel('''//trim(labels(ix3+2))//trim(fmt) write (50,*) 'set(gca,''FontSize'',axes_fontsize);' fmt = ''',fontsize=lab_fontsize)' write (80,'(a)') 'xlabel('''//trim(labels(ix1+2))//trim(fmt) !JZ write (80,'(a)') 'ylabel('''//trim(labels(ix2+2))//trim(fmt) !JZ write (80,'(a)') 'xticks(fontsize=axes_fontsize)' !JZ write (80,'(a)') 'hbar = colorbar(orientation=''horizontal'')' !JZ write (80,'(a)') 'hbar.set_xlabel('''//trim(labels(ix3+2))//trim(fmt) !JZ !JZ Not done this bit if (num_3D_plots > 2 .and. matlab_version < 7) then write (50,*) ' p = get(pa,''Position'');' write (50,*) 'set(pa,''Position'',[p(1) (p(2)+p(4)/8) p(3) p(4)]);' end if end do !JZ Don't need this write (50,*) 'set(gcf, ''PaperUnits'',''inches'');' !JZ Agg handles this too, I think. tmpstr = RealToStr((plot_col*9.)/plot_row) write (50,*) 'set(gcf, ''PaperPosition'',[ 0 0 '// & trim(tmpstr) //' 8]);' write (50,*) 'print -dpsc2 '//trim(rootname)//'_3D.ps;' write(80,'(a)') 'savefig('''//trim(rootname)//'_3D.ps'')' close(50) close(80) end if !Limits from global likelihood open(unit=50,file=trim(rootname)//'.likestats',form='formatted',status='replace') write (50,*) 'Best fit sample -log(Like) = ',coldata(2,bestfit_ix) write (50,*) '' write(50,'(a)') 'param bestfit lower1 upper1 lower2 upper2' do j=1, num_vars write(50,'(1I5,5E14.6," '//trim(labels(colix(j)))//'")') colix(j)-2, coldata(colix(j),bestfit_ix),& minval(coldata(colix(j),0:ND_cont1)), & maxval(coldata(colix(j),0:ND_cont1)), & minval(coldata(colix(j),0:ND_cont2)), & maxval(coldata(colix(j),0:ND_cont2)) end do close(50) end program GetDist