!###############################BEGIN GMCMC############################# ! General Markov Chain Monte Carlo routines ! author: Guthrie Miller guthriemiller@gmail.com ! reference: http://www.lulu.com/content/paperback-book/probabilistic-interpretation-of-data/13620353 module GMCMC_MOD character*10, save :: tapefile='tape.txt' ! real*8 esumu or esumv is written on tape file in hex ! real*8 thetu or thetv is written on tape file in hex ! one extra character to provide spacing and readability character*16, save :: cesum character*17, save, dimension(:), allocatable :: cthet character*39, save :: info ! date, time, size info about file character*39, save :: infoid ! file info about datafile character*80, save :: dataid ! name of datafile logical, save:: accept,ranwalk,slow0,slow=.true. logical, save, dimension(:), allocatable :: bad integer, save :: ialg,icalcvx,icycle,i 1 ,idorun,ifocus,ihypu,ihypv 1 ,incr(2),incracpt(2),incrhyp(2,2),iprbmax,irecmax(2),irpt 1 ,iseed,iv,ivmin,nburn,ndebug,ngrps,nhyp=1,nprm 1 ,nprmsgrp,nrpts,nsumlike,ntape,nv,nvgood integer, save ,dimension(:), allocatable::iv1good,nprmgrp,iwrap integer, save, dimension(:,:), allocatable::iiprm integer*8, save :: irun,irunacpt,irundbug 1 ,nrun,nrun0,nthin real*8, save :: burnfrac,c2(2,2),c2f,esumu,eminv 1 ,eminv0,c2ndref=1,esumref,esumrefp,fracrw,phyp(2) 1 ,proddthet,prodsum,pswitch=0.1d0,sumlike,xsec,xsec0 1 ,xrand,One=1 real*8 dranr250 external dranr250 real*8, save, dimension(:), allocatable :: esumv 1 ,pfocus,pfocusn,pfocusc,prbnd,prbcd,thetu,dthet real*8, save, dimension(:,:), allocatable :: thetv end module GMCMC_MOD !----------------------------------------------------------------------- ! the user calls subroutine GMCMC below and needs to supply following five routines: ! calcvx calculate extra v (candidate) variables xv, ! which are the natural variables of the problem, ! from the fundamental theta variables ! and calculate energy = chi**2/2 ! movex move xv variables to xu variables ! when candidate is accepted and chain moved ! mvbackvx(iv1) restore xv variables from xu for candidate iv1 ! increm increment tallies ! zero zero tallies ! theta variables written to tapefile (TAPE.TXT), ! which contains NTAPE records uniformly spaced throughout entire run ! Other output also written to file GMCMC.TXT as well as to terminal ! random numbers should be generated using e.g. xrand=dranr250(0) !----------------------------------------------------------------------- subroutine GMCMC(nrunIN,ntapeIN,nrptsIN,ndebugIN 1 ,mdataIN,iseedIN,icycleIN,nvIN,ialgIN 1 ,slowIN,burnfracIN,fracrwIN 1 ,nprmIN,ngrpsIN,nprmsgrpIN,nprmgrpIN,iiprmIN 1 ,pfocusIN,iwrapIN,thetIN,dthetIN 1 ,datafileIN) ! this program validates input parameters, initializes variables ! and initiates run by calling runchain ! arguments: ! NRUN number of chain iterations ! NTAPE number of chain iterations written on file tapefile ! NRPTS number of reports written out during the run ! NDEBUG for debugging, number of detailed print outs ! MDATA number of data points (for calculation of ChiSq/mDATA) ! ISEED random number seed ! ICYCLE = 1 for parameter movement by regular cycling through groups ! = 2 for randomly chosen next group ! = 3 for randomly chosen with probability of group PFOCUS(IGRP) ! NV number of candidates ! IALG = 1 for MRT algorithm ! = 2 for B algorithm ! = 3 for multiple candidate (MC) algorithm ! SLOW logical variable that if true uses slow but sure calculations ! without tricky speedups ! BURNFRAC the fraction of total iterations that are ignored in tallies ! because of possible initialization bias ! FRACRW fraction of time random walk is used rather than choosing ! candidate from entire space ! NPRM number of parameters ! NGRPS number of parameter groups ! NPRMSGRP maximum value of NPRMGRP(IGRP) for array allocation ! NPRMGRP(IGRP), IGRP=1,NGRPS number of parameters withing group IGRP ! IIPRM(IPRMGRP,IGRP), IPRMGRP=1,NPRMGRP(IGRP) parameters withing group IGRP ! PFOCUS(IGRP) if NGRPS>1, probability of focusing on group IGRP, when ICYCLE=3 ! IWRAP(IPRM) wraping of parameter IPRM, 1 no wrapping, 2 cyclical wrapping ! when parameter goes outside of theta interval (0,1) ! THET(IPRM) intial value of theta parameter IPRM ! DTHET(IPRM) random walk half width for parameter IPRM ! DATAFILE name of datafile, used with checksum to see if data have changed ! ! These parameters for NHYP=2 can be changed by referencing them in GMCMC_MOD ! NHYP = 1 for no alternate (reference) hypotheses (default) ! = 2 for reference hypothesis ! if NHYP = 2, these parameters used: ! C2NDREF equivalent Chisq/nData to give constant energy when in the reference hypothesis ! PSWITCH probability to have proposal switch hypothesis from normal to reference use GMCMC_MOD implicit none integer ialgIN,icycleIN,iseedIN,mdataIN,ndebugIN,ngrpsIN 1 ,nprmIN,nprmsgrpIN,nrptsIN,ntapeIN,nvIN,ireturn integer iiprmIN(nprmsgrpIN,ngrpsIN),iwrapIN(nprmIN) 1 ,nprmgrpIN(ngrpsIN) integer*8 nrunIN logical slowIN real*8 burnfracIN,fracrwIN real*8 dthetIN(nprmIN),pfocusIN(ngrpsIN),thetIN(nprmIN) character*(*) datafileIN character*80 ermess1,ermess2 logical ex integer ih,igrp,iprm,iprmgrp,iv1 open(1,file='gmcmc.out') write(*,*)'-------------------begin GMCMC---------------------' write(1,*)'-------------------begin GMCMC---------------------' if(iseedIN.eq.0) then write(ermess1,*)'iseed=',iseedIN ermess2='random number seed must be positive integer' go to 999 endif iseed=iseedIN if(ndebugIN.gt.0) then write(*,*)'GMCMC: ndebug=',ndebugIN 1 ,' number of debug printouts' write(1,*)'GMCMC: ndebug=',ndebugIN 1 ,' number of debug printouts' endif ndebug=ndebugIN if(slowIN) then write(*,*)'GMCMC: REMINDER slow=',slowIN 1 ,' means always use slow but sure calculations' write(1,*)'GMCMC: REMINDER slow=',slowIN 1 ,' means always use slow but sure calculations' endif slow=slowIN slow0=slow if(burnfracIN.eq.0.or.burnfracIN.gt.1.) then write(ermess1,*)'burnfrac=',burnfracIN ermess2='burn fraction should be small fraction, i.e. 0.05' go to 999 endif burnfrac=burnfracIN if(fracrwIN.lt.0.or.fracrwIN.gt.1.) then write(ermess1,*)'fracrw=',fracrwIN ermess2='fraction of time to random walk, between 0 and 1' go to 999 endif fracrw=fracrwIN if(nrunIN.eq.0) then write(ermess1,*)'nrun=',nrunIN ermess2='run length, number of mcmc iterations, not set' go to 999 endif nrun=nrunIN if(ntapeIN.eq.0.or.ntapeIN.gt.nrunIN) then write(ermess1,*)'ntape=',ntapeIN ermess2='number of iterations to write on tape bad' go to 999 endif ntape=ntapeIN if(nrptsIN.eq.0.or.nrptsIN.gt.ntapeIN) then write(*,*)'GMCMC: WARNING nrpts bad, nrpts=',nrptsIN write(1,*)'GMCMC: WARNING nrpts bad, nrpts=',nrptsIN endif nrpts=nrptsIN if(mdataIN.eq.0) then write(ermess1,*)'mdata=',mdataIN ermess2='number of data points' go to 999 endif ! c2f is coefficient that converts energy (sum of squared residuals) ! to chisq per data point c2f=2.d0/mdataIN if(nhyp.eq.2) then write(*,*)'GMCMC: nhyp=',nhyp 1 ,' zero-dimensional reference hypothesis' write(1,*)'GMCMC: nhyp=',nhyp 1 ,' zero-dimensional reference hypothesis' if(pswitch.lt.0.or.pswitch.gt.1) then write(ermess1,*)'ih,pswitch=',ih,pswitch ermess2='switch probabiltiy bad' 1 //' should be >= 0 and <= 1' go to 999 endif if(c2ndref.lt.0) then write(ermess1,*)'c2ndref=',c2ndref ermess2='need positive c2ndref=c2f*esumref' go to 999 endif esumref=c2ndref/c2f esumrefp=esumref+log(pswitch) ! esumrefp=-log(P0/pswitch) write(*,*)'GMCMC: c2ndref (Equivalent Chisq/nData) =' 1 ,c2ndref,' esumref=',esumref write(1,*)'GMCMC: c2ndref (Equivalent Chisq/nData) =' 1 ,c2ndref,' esumref=',esumref write(*,*)' esumrefp=esumref+log(pswitch)=',esumrefp 1 ,' Equivalent Chisq/nData=',c2f*esumrefp write(1,*)' esumrefp=esumref+log(pswitch)=',esumrefp 1 ,' Equivalent Chisq/nData=',c2f*esumrefp elseif(nhyp.eq.1) then else write(ermess1,*)'nhyp bad, nhyp=',nhyp ermess2='nhyp = 1 or 2' go to 999 endif c write(*,*)'datafile=',trim(datafile) if(trim(datafileIN).ne.'') then inquire(file=datafileIN,exist=ex) if(.not.ex) then write(ermess1,*)'datafile=',trim(datafileIN) ermess2="datafile doesn't exist" go to 999 endif else write(*,*)'GMCMC: WARNING datafile=' 1 ,trim(datafileIN),' blank' write(1,*)'GMCMC: WARNING datafile=' 1 ,trim(datafileIN),' blank' endif dataid=datafileIN c write(*,*)'ngrps=',ngrps if(nprmIN.le.0) then write(ermess1,*)'nprm=',nprmIN ermess2='nprm bad' go to 999 endif nprm=nprmIN if(nvIN.le.0) then write(ermess1,*)'nv=',nvIN ermess2='number of candidates bad' go to 999 endif nv=nvIN if(ngrpsIN.le.0) then write(ermess1,*)'ngrps=',ngrpsIN ermess2='ngrps bad' go to 999 endif ngrps=ngrpsIN c nogrps=.true. if(nprmsgrpIN.le.0) then write(ermess1,*)'nprmsgrp=',nprmsgrpIN ermess2='nprmsgrp bad' go to 999 endif nprmsgrp=nprmsgrpIN if(allocated(iiprm))deallocate(iiprm) if(allocated(nprmgrp))deallocate(nprmgrp) if(allocated(thetv))deallocate(thetv) if(allocated(thetu))deallocate(thetu) if(allocated(cthet))deallocate(cthet) if(allocated(bad))deallocate(bad) if(allocated(dthet))deallocate(dthet) if(allocated(iwrap))deallocate(iwrap) if(allocated(iv1good))deallocate(iv1good) if(allocated(esumv))deallocate(esumv) if(allocated(prbnd))deallocate(prbnd) if(allocated(prbcd))deallocate(prbcd) allocate(iiprm(nprmsgrp,ngrps),nprmgrp(ngrps)) allocate(thetv(nprm,nv),thetu(nprm),cthet(nprm),bad(nv)) allocate(dthet(nprm),iwrap(nprm)) allocate(iv1good(nv),esumv(nv)) allocate(prbnd(nv+1),prbcd(nv+1)) c write(*,*)'GMCMC: after allocate' if(ngrps.gt.1.and. 1 (icycleIN.le.0.or.icycleIN.gt.3)) then write(ermess1,*)'ngrps,icycle=' 1 ,ngrps,icycleIN ermess2='icycle=1 (cycle), 2 (random)'// 1 ', 3 (prob specified by pfocus)' go to 999 endif icycle=icycleIN igrpLoop: do igrp=1,ngrps c write(*,*)'igrp,nprmgrp=',igrp,nprmgrp(igrp) if(nprmgrpIN(igrp).eq.0) then write(ermess1,*)'nprmgrp(igrp), igrp=' 1 ,nprmgrpIN(igrp),igrp ermess2='at least 1 parameter in any group' go to 999 elseif(nprmgrpIN(igrp).gt.nprmsgrp) then write(ermess1,*)'nprmgrp(igrp), igrp, nprmsgrp=' 1 ,nprmgrpIN(igrp),igrp, nprmsgrp ermess2='too many parameters' go to 999 endif nprmgrp(igrp)=nprmgrpIN(igrp) do iprmgrp=1,nprmgrp(igrp) if(iiprmIN(iprmgrp,igrp).lt.1.or. 1 iiprmIN(iprmgrp,igrp).gt.nprm) then write(ermess1,*)'iiprm,iprmgrp,igrp,nprm=' 1 ,iiprmIN(iprmgrp,igrp),iprmgrp,igrp,nprm ermess2='iiprm needs to be between 1 and nprm' go to 999 endif iiprm(iprmgrp,igrp)=iiprmIN(iprmgrp,igrp) iprm=iiprm(iprmgrp,igrp) if(iwrapIN(iprm).lt.1.or.iwrapIN(iprm).gt.2) then write(ermess1,*)'iwrap(iprm), iprm=',iwrapIN(iprm),iprm ermess2='iwrap = 1 (no wrap), or 2 (wrap)' go to 999 endif enddo enddo igrpLoop iprmLoop: do iprm=1,nprm do igrp=1,ngrps do iprmgrp=1,nprmgrp(igrp) if(iprm.eq.iiprm(iprmgrp,igrp)) cycle iprmLoop enddo enddo write(ermess1,*)'iprm=',iprm ermess2='iprm does not appear in iiprm' go to 999 enddo iprmLoop if(icycle.eq.3) then allocate(pfocus(ngrps),pfocusn(ngrps),pfocusc(ngrps)) if(minval(pfocusIN(1:ngrps)).le.0) then write(ermess1,*)'minval(pfocus).le.0, pfocus=' 1 ,pfocusIN(1:ngrps) ermess2='for icycle=3, pfocus needs to be >0' go to 999 endif pfocus=pfocusIN(1:ngrps) pfocusn=pfocus/sum(pfocus) ! nomalized probability do igrp=1,ngrps pfocusc(igrp)=sum(pfocusn(:igrp)) ! cumulative probability enddo pfocusc(ngrps)=1 ! in case it isn't because of roundoff errors endif if(ialgIN.lt.1.or.ialgIN.gt.3) then write(ermess1,*)'ialg=',ialgIN ermess2='ialg = 1,2,3 for MRT, B, MC' go to 999 endif if(ialgIN.eq.3.and.nhyp.eq.2) then write(*,*)'GMCMC: WARNING ialg=',ialgIN 1 ,' should not be used with nhyp=2' write(1,*)'GMCMC: WARNING ialg=',ialgIN 1 ,' should not be used with nhyp=2' endif ialg=ialgIN ! initialize do iprm=1,nprm if(thetIN(iprm).lt.0.or.thetIN(iprm).gt.1) then write(ermess1,*)'thet(iprm), iprm=',thetIN(iprm),iprm ermess2='starting theta value required to be between 0 and 1' go to 999 endif if(dthetIN(iprm).le.0.or.dthetIN(iprm).gt.1) then write(ermess1,*)'iprm,dthet(iprm)=' 1 ,iprm,dthetIN(iprm) ermess2='dthet required to be > 0 and <= 1' go to 999 endif dthet(iprm)=dthetIN(iprm) iwrap(iprm)=iwrapIN(iprm) thetu(iprm)=thetIN(iprm) enddo do igrp=1,ngrps proddthet=1 do iprmgrp=1,nprmgrp(igrp) iprm=iiprm(iprmgrp,igrp) proddthet=proddthet*min(2*dthet(iprm),1d0) enddo enddo nsumlike=0 sumlike=0 write(*,*)' Parameterization:' write(1,*)' Parameterization:' if(icycle.eq.3) then write(*,*)' iprm igrp iprmgrp iwrap' 1 ,' pfocus thet0 dthet' write(1,*)' iprm igrp iprmgrp iwrap' 1 ,' pfocus thet0 dthet' else write(*,*)' iprm igrp iprmgrp iwrap' 1 ,' thet0 dthet' write(1,*)' iprm igrp iprmgrp iwrap' 1 ,' thet0 dthet' endif do igrp=1,ngrps do iprmgrp=1,nprmgrp(igrp) iprm=iiprm(iprmgrp,igrp) if(icycle.eq.3) then write(*,'(4i8,2f8.5,1p2e9.2)') 1 iprm,igrp,iprmgrp,iwrap(iprm) 1 ,pfocusn(igrp),thetu(iprm),dthet(iprm) write(1,'(4i8,2f8.5,1p2e9.2)') 1 iprm,igrp,iprmgrp,iwrap(iprm) 1 ,pfocusn(igrp),thetu(iprm),dthet(iprm) else write(*,'(4i8,2f8.5,1p2e9.2)') 1 iprm,igrp,iprmgrp,iwrap(iprm) 1 ,thetu(iprm),dthet(iprm) write(1,'(4i8,2f8.5,1p2e9.2)') 1 iprm,igrp,iprmgrp,iwrap(iprm) 1 ,thetu(iprm),dthet(iprm) endif enddo enddo write(*,*) write(1,*) ifocus=1 nvgood=1 ihypv=1 do iv1=1,nvgood bad(iv1)=.false. iv1good(iv1)=iv1 do igrp=1,ngrps do iprmgrp=1,nprmgrp(igrp) iprm=iiprm(iprmgrp,igrp) thetv(iprm,iv1)=thetu(iprm) c write(*,*)'GMCMC: iv1,iprm,thetv=',iv1,iprm,thetv(iprm,iv1) enddo enddo enddo slow=.true. write(*,*)'GMCMC: before first call to CALCVX' write(1,*)'GMCMC: before first call to CALCVX' call calcvx iv=1 write(*,*)'GMCMC: before first call to MOVE' write(1,*)'GMCMC: before first call to MOVE' call move ! move moves v variables for iv to u variables c write(*,*)'after move' ! restore slow slow=slow0 write(*,*)'GMCMC: first value, Chisq/nData=',c2f*esumu write(1,*)'GMCMC: first value, Chisq/nData=',c2f*esumu if(esumu.ne.esumu) then write(*,*)'GMCMC: Error, Chisq/nData bad, esumu=',esumu write(1,*)'GMCMC: Error, Chisq/nData bad, esumu=',esumu stop endif write(*,*)'GMCMC: initial parameter check complete!' write(1,*)'GMCMC: initial parameter check complete!' c write(*,*)'slow=',slow0 call chktape(ireturn) if(ireturn.eq.1) return call runchain call finalze return 999 continue write(*,*)'GMCMC: ERROR ',trim(ermess1),' ',trim(ermess2) write(1,*)'GMCMC: ERROR ',trim(ermess1),' ',trim(ermess2) stop end !----------------------------------------------------------------------- subroutine chktape(ireturn) ! checks for the existence of a tape file ! if one exists and the data or requested number of chain ! iterations haven't changed, just read from tape file ! without recalculating chain use GMCMC_MOD, only : burnfrac,dataid,idorun,iseed,xrand,dranr250 1 ,info,infoid,irun,irundbug,ndebug,nrun,ntape,tapefile implicit none integer ichk,ireturn,if1,ib1,ibuf1 1 ,nrunfac,ntapetape,nthintape,nprmtape character*80 dataidtape character*120 cmd character*50 mess logical chkfile integer*8 nruncalc,nruntape common/r250com/if1,ib1,ibuf1(250) c write(*,*) c 1 'CHKTAPE: ntape,nrun,nrun0,nrpts,irpt,nthin,nburn=' c 1 ,ntape,nrun,nrun0,nrpts,irpt,nthin,nburn ireturn=0 ! set up debug printout if(ndebug.gt.0) then irundbug=nrun/ndebug endif c write(*,*)'dataid=',dataid if(trim(dataid).ne.'') then ! get identification of data used cmd='chksum$ '//trim(dataid)// ' > dataid.txt' c write(*,*)'cmd = ',cmd call shellcmd(cmd) if(chkfile('dataid.txt',infoid)) then open(167,file='dataid.txt') read(167,'(a80)')dataid write(*,*)'CHKTAPE: dataid=',dataid write(1,*)'CHKTAPE: dataid=',dataid close(167) endif endif idorun=0 write(*,*)'CHKTAPE: Check for existing ',trim(tapefile) write(1,*)'CHKTAPE: Check for existing ',trim(tapefile) if(chkfile(trim(tapefile),info)) then open(170,file=trim(tapefile),action='read') write(*,*)'CHKTAPE: reading file: ',trim(tapefile) 1 ,' info=',info write(1,*)'CHKTAPE: reading file: ',trim(tapefile) 1 ,' info=',info read(170,*)nthintape,ntapetape,nruntape,nprmtape read(170,'(a80)')dataidtape write(*,*)' dataidtape=',trim(dataidtape) write(1,*)' dataidtape=',trim(dataidtape) close(170) nrunfac=nrun/nruntape nruncalc=nrunfac nruncalc=nruncalc*nruntape c write(*,*)'nrun,nruntape,nrunfac=',nrun,nruntape,nrunfac c write(*,*)'nrun,nruncalc=',nrun,nruncalc if(dataid.eq.dataidtape.and.nrun.eq.nruntape) then ichk=1 mess=' with same nrun, dataid' else ichk=0 mess=', use it' endif if(ichk.eq.1) then ! use existing tape, skip runchain write(*,*)'CHKTAPE: existing ',trim(tapefile),trim(mess) write(1,*)'CHKTAPE: existing ',trim(tapefile),trim(mess) ! initialize random number generator r250 to correspond to end of chain open(169,file=tapefile) read(169,*) ! skip two lines read(169,*) 100 read(169,*,err=101)irun ! error will be caused by presence of text " rand250 initialization:" go to 100 101 continue ! the is the rand250 info needed for continuation of sequence of random numbers read(169,*)if1,ib1,ibuf1 close(169) call finalze ireturn=1 return elseif(nrun.eq.nruncalc.and. 1 dataid.eq.dataidtape) then ! extend existing tape write(*,'('' CHKTAPE: existing '',a9,'', nrun ='',i2 1 ,'' x nrun tape same dataid'')')trim(tapefile),nrunfac write(1,'('' CHKTAPE: existing '',a9,'', nrun ='',i2 1 ,'' x nrun tape same dataid'')')trim(tapefile),nrunfac idorun=2 return else ! rerun, making new tape write(*,*)'CHKTAPE: ',trim(tapefile) 1 ,' exists but nrun or dataid not right for continuation' write(*,*)' on tape nrun = ',nruntape 1 ,' data id = ',dataidtape write(1,*)'CHKTAPE: ',trim(tapefile) 1 ,' exists but nrun or dataid not right for continuation' write(1,*)' on tape nrun = ',nruntape 1 ,' data id = ',dataidtape write(*,*)' deleting ',trim(tapefile),', new run' write(1,*)' deleting ',trim(tapefile),', new run' cmd='del '//trim(tapefile) call shellcmd(cmd) idorun=1 return endif endif idorun=1 return end !----------------------------------------------------------------------- subroutine runchain ! start or continue mcmc run ! this routine handles the "tape file" ! deciding whether to start new run or continue old one use GMCMC_MOD, only : burnfrac,dataid,idorun 1 ,ihypu,ihypv,irpt,irun,irundbug,iseed,iv,iv1good 1 ,nburn,ndebug,nprm,nrpts,nrun,ntape,nthin 1 ,nv,nvgood,slow,slow0,tapefile,xsec0 implicit none integer time character*80 dataidtape real*8 xrand,dranr250 character*120 cmd integer i,ibuf1,ib1,iend,if1,iwrt0 1 ,nrunfac,nrun0,ncopy,nnew,nvsav,nwrt0 integer*8 irunrpt common/r250com/if1,ib1,ibuf1(250) if(idorun.eq.1) then go to 1 elseif(idorun.eq.2) then go to 2 else write(*,*)'RUNCHAIN: should not come here, idorun=',idorun write(1,*)'RUNCHAIN: should not come here, idorun=',idorun stop endif 1 continue ! come here to start new run nthin=nrun/ntape nrun=nthin nrun=nrun*ntape nburn=nint(burnfrac*nrun) c write(*,*)'nthin,nrun,nburn=',nthin,nrun,nburn write(*,*)'RUNCHAIN: new MCMC run of ',nrun,' iterations' write(*,'('' write on tape'',i5 1 ,'' records total (one every '' 1 ,i10,'' iterations)'')')ntape,nthin write(1,*)'RUNCHAIN: new MCMC run of ',nrun,' iterations' write(1,'('' write on tape'',i5 1 ,'' records total (one every '' 1 ,i10,'' iterations)'')')ntape,nthin xsec0=time() c write(*,*)'nrpts=',nrpts irunrpt=nrun/nrpts c write(*,*)'irunrpt,nrpts=',irunrpt,nrpts irpt=0 irun=0 ! open new tape file open(168,file='new'//trim(tapefile)) write(168,*)nthin,ntape,nrun,nprm write(168,'(a80)')dataid ! initialize random number generator r250 xrand=dranr250(iseed) write(*,*)'RUNCHAIN: first random number=',xrand 1 ,' iseed=',iseed write(1,*)'RUNCHAIN: first random number=',xrand 1 ,' iseed=',iseed slow=.true. ! this do loop is the actual mcmc run do i=1,nrun call run1 ! first run with slow = true if(i.eq.1) slow=slow0 if(mod(irun,nthin).eq.0) then ! write record on tape call wrt1u(168) call increm0 endif if(mod(irun,irunrpt).eq.0) then ! write one of nrpts reports call wrtrpt endif enddo ! write rand250 initialization info to allow continuation write(168,*)'rand250 initialization:' write(168,*)if1,ib1,ibuf1 close(168) write(*,*)'RUNCHAIN: copy preliminary new'//trim(tapefile) 1 //' -> '//trim(tapefile) write(1,*)'RUNCHAIN: copy preliminary new'//trim(tapefile) 1 //' -> '//trim(tapefile) cmd='copy new'//trim(tapefile)//' '//trim(tapefile)//' /y' call shellcmd(cmd) cmd='del new'//trim(tapefile) call shellcmd(cmd) return 2 continue ! come here if previous run with good tapefile open(169,file=tapefile) read(169,*)iwrt0,nwrt0,nrun0 read(169,'(a80)')dataidtape write(*,*)'RUNCHAIN: reading file: ',trim(tapefile) 1 ,' dataid=',trim(dataidtape) write(1,*)'RUNCHAIN: reading file: ',trim(tapefile) 1 ,' dataid=',trim(dataidtape) ! new run will be some integer nrunfac times as long nrunfac=nrun/nrun0 c write(*,*)'nrunfac,nrunfac*nrun0=',nrunfac,nrunfac*nrun0 if(nrunfac.gt.1.and.nrun.eq.nrunfac*nrun0.and. 1 dataid.eq.dataidtape) then write(*,*)'RUNCHAIN: ',trim(tapefile),' OK for continuation' write(1,*)'RUNCHAIN: ',trim(tapefile),' OK for continuation' else write(*,*)'RUNCHAIN: ',trim(tapefile) 1 ,': nrun or dataid not right to allow continuation,' 1 ,' start new run' write(*,*)' nrun needs to be ',nrunfac*nrun0 write(*,*)' data id from tape = ',dataidtape write(1,*)'RUNCHAIN: ',trim(tapefile) 1 ,': nrun or dataid not right to allow continuation,' 1 ,' start new run' write(1,*)' nrun needs to be ',nrunfac*nrun0 write(1,*)' data id from tape = ',dataidtape close(169) write(*,*)' deleting ',trim(tapefile) write(1,*)' deleting ',trim(tapefile) cmd='del '//trim(tapefile) call shellcmd(cmd) go to 1 endif nthin=nrun/ntape nrun=nthin nrun=nrun*ntape nburn=nint(burnfrac*nrun) c write(*,*)'nthin,nrun,nburn=',nthin,nrun,nburn open(168,file='new'//trim(tapefile)) write(168,*)nthin,ntape,nrun,nprm write(168,'(a80)')dataidtape ncopy=ntape/nrunfac nnew=ntape-ncopy write(*,*)'RUNCHAIN: extending mcmc run to ',nrun,' iterations' write(*,'('' write '',i5,'' records total (one every '' 1 ,i10,'' iterations)'')')ntape,nthin write(*,'('' copy '',i4,'' records from old tape'')')ncopy write(*,'('' write '',i4,'' new records'')')nnew write(1,*)'RUNCHAIN: extending mcmc run to ',nrun,' iterations' write(1,'('' write '',i5,'' records total (one every '' 1 ,i10,'' iterations)'')')ntape,nthin write(1,'('' copy '',i4,'' records from old tape'')')ncopy write(1,'('' write '',i4,'' new records'')')nnew xsec0=time() irunrpt=nrun/nrpts c write(*,*)'nrun,nrpts,nrpts,irunrpt=',nrun,nrpts,nrpts,irunrpt irpt=0 nvsav=nv nv=1 iv=1 nvgood=1 iv1good(1)=1 ! read record from tape 100 read(169,*,err=101)irun ! error will be caused by presence of text " rand250 initialization:" c write(*,*)'reading record from old tape' c write(*,*)'irun=',irun ! get full info from tape every nthin records if(mod(irun,nthin).eq.0) then backspace(169) if(irun.eq.nrun0) then call rd1tov(169,iend) slow=.true. call calcvx call move else call rd1tou(169,iend) endif call wrt1u(168) call increm0 endif go to 100 101 continue ! initialize rand250 ! the is the rand250 info needed for continuation of sequence ! of random numbers read(169,*)if1,ib1,ibuf1 close(169) if(irun.ne.nrun0) then write(*,*)'RUNCHAIN: irun.ne.nrun0, irun,nrun0=',irun,nrun0 write(1,*)'RUNCHAIN: irun.ne.nrun0, irun,nrun0=',irun,nrun0 stop endif nv=nvsav ! continue mcmc run do i=nrun0+1,nrun ! this loop continues old run call run1 ! first run slow, because v variables not set for all iv1 ! can't restore for igrp=ifocus ! after first run restore slow if(irun.eq.nrun0+1) slow=slow0 if(mod(irun,nthin).eq.0) then ! write record on tape call wrt1u(168) call increm0 endif if(mod(irun,irunrpt).eq.0) call wrtrpt enddo ! write rand250 initialization info to allow run continuation write(168,*)'rand250 initialization:' write(168,*)if1,ib1,ibuf1 close(168) close(169) write(*,*)'RUNCHAIN: copy preliminary new'//trim(tapefile) 1 //' -> ',trim(tapefile) write(1,*)'RUNCHAIN: copy preliminary new'//trim(tapefile) 1 //' -> ',trim(tapefile) cmd='copy new'//trim(tapefile)//' '//trim(tapefile) call shellcmd(cmd) cmd='del new'//trim(tapefile) call shellcmd(cmd) write(*,*)'file ',trim(tapefile),' written' write(1,*)'file ',trim(tapefile),' written' end !----------------------------------------------------------------------- subroutine increm0 ! increment talleys for wrtrpt use GMCMC_MOD, only : accept,c2,c2f,esumu,ihypu,incr 1 ,incracpt,incrhyp,irun,nburn,nhyp implicit none integer ih c write(*,*)'INCREM: irun,nburn=',irun,nburn incr(1)=incr(1)+1 if(accept) incracpt(1)=incracpt(1)+1 ih=ihypu+nhyp-1 c write(*,*)'INCREM0: ih=',ih incrhyp(ih,1)=incrhyp(ih,1)+1 c2(ih,1)=c2(ih,1)+c2f*esumu if(irun.le.nburn) return incr(2)=incr(2)+1 if(accept) incracpt(2)=incracpt(2)+1 incrhyp(ih,2)=incrhyp(ih,2)+1 c2(ih,2)=c2(ih,2)+c2f*esumu c write(*,*)'INCREM: incr,ihypu,c2/incrhyp=',incr,ihypu,c2/incrhyp end !----------------------------------------------------------------------- subroutine zero0 ! zero built-in tallies use GMCMC_MOD, only : c2,incr,incracpt,incrhyp 1 ,irun,irunacpt,nhyp,irpt implicit none irun=0 irunacpt=0 incr=0 incracpt=0 c2=0 incrhyp=0 end !----------------------------------------------------------------------- subroutine wrtrpt use GMCMC_MOD, only : c2,incr,incracpt,incrhyp 1 ,irpt,irun,nhyp,nburn,nrpts,nrun,nv,xsec,xsec0 implicit none ! write report. This is done nrpts times throughout the entire run integer time integer ih integer,save ::incracpt0,incrhyp0(2),incr10 real*8 ,save :: c2sum0(2),ttaken real*4 tend if(irpt.eq.0) then write(*,*)'WRTRPT: Number of chain iterations = ',nrun write(*,*)' Number of candidates per chain iteration = ',nv write(*,*)' Number reports = ',nrpts write(1,*)'WRTRPT: Number of chain iterations = ',nrun write(1,*)' Number of candidates per chain iteration = ',nv write(1,*)' Number reports = ',nrpts if(nhyp.eq.2) then write(*,'('' rpt# irun Chisq/ND(hyp1,0) '' 1 ,''time(min) #tape'' 1 ,'' acptd NormalHyp RefHyp'')') write(1,'('' rpt# irun Chisq/ND(hyp1,0) '' 1 ,''time(min) #tape'' 1 ,'' acptd NormalHyp RefHyp'')') else write(*,'('' rpt# irun Chisq/ND time(min) #tape'' 1 ,'' acptd'')') write(1,'('' rpt# irun Chisq/ND time(min) #tape'' 1 ,'' acptd'')') endif tend=0 incr10=0 incrhyp0=0 incracpt0=0 c2sum0=0 endif xsec=time() ! system call giving clock time in seconds irpt=irpt+1 tend=(nrun-irun)*(xsec-xsec0)/(60*irun) c write(*,*)'irpt,xsec,incr(1),incr10=',irpt,xsec,incr(1),incr10 if(nhyp.eq.2) then write(*,'(i10,i11,1p2g10.3,0pf10.2,2i9,2i10)') 1 irpt,irun 1 ,((c2(ih,1)-c2sum0(ih))/max((incrhyp(ih,1)-incrhyp0(ih)),1) 1 ,ih=nhyp,1,-1) 1 ,tend 1 ,incr(1)-incr10,incracpt(1)-incracpt0 1 ,(incrhyp(ih,1)-incrhyp0(ih),ih=nhyp,1,-1) write(1,'(i10,i11,1p2g10.3,0pf10.2,2i9,2i10)') 1 irpt,irun 1 ,((c2(ih,1)-c2sum0(ih))/max((incrhyp(ih,1)-incrhyp0(ih)),1) 1 ,ih=nhyp,1,-1) 1 ,tend 1 ,incr(1)-incr10,incracpt(1)-incracpt0 1 ,(incrhyp(ih,1)-incrhyp0(ih),ih=nhyp,1,-1) else write(*,'(i10,i11,1p1g10.3,0pf10.2,2i9,2i10)') 1 irpt,irun 1 ,((c2(ih,1)-c2sum0(ih))/max((incrhyp(ih,1)-incrhyp0(ih)),1) 1 ,ih=nhyp,1,-1) 1 ,tend 1 ,incr(1)-incr10,incracpt(1)-incracpt0 write(1,'(i10,i11,1p1g10.3,0pf10.2,2i9,2i10)') 1 irpt,irun 1 ,((c2(ih,1)-c2sum0(ih))/max((incrhyp(ih,1)-incrhyp0(ih)),1) 1 ,ih=nhyp,1,-1) 1 ,tend 1 ,incr(1)-incr10,incracpt(1)-incracpt0 endif incracpt0=incracpt(1) incrhyp0=incrhyp(:,1) c2sum0=c2(:,1) incr10=incr(1) if(irun.ne.nrun) return write(*,*)' Totals (excluding burnfrac):' write(1,*)' Totals (excluding burnfrac):' if(nhyp.eq.2) then write(*,'(i10,i11,1p2g10.3,0pf10.2,2i9,2i10)') 1 irpt,irun-nburn 1 ,(c2(ih,2)/max(incrhyp(ih,2),1),ih=nhyp,1,-1) 1 ,tend 1 ,incr(2),incracpt(2) 1 ,(incrhyp(ih,2),ih=nhyp,1,-1) write(1,'(i10,i11,1p2g10.3,0pf10.2,2i9,2i10)') 1 irpt,irun-nburn 1 ,(c2(ih,2)/max(incrhyp(ih,2),1),ih=nhyp,1,-1) 1 ,tend 1 ,incr(2),incracpt(2) 1 ,(incrhyp(ih,2),ih=nhyp,1,-1) else write(*,'(i10,i11,1pg10.3,0pf10.2,2i9,2i10)') 1 irpt,irun-nburn 1 ,(c2(ih,2)/max(incrhyp(ih,2),1),ih=nhyp,1,-1) 1 ,tend 1 ,incr(2),incracpt(2) write(1,'(i10,i11,1pg10.3,0pf10.2,2i9,2i10)') 1 irpt,irun-nburn 1 ,(c2(ih,2)/max(incrhyp(ih,2),1),ih=nhyp,1,-1) 1 ,tend 1 ,incr(2),incracpt(2) endif xsec=time() ttaken=(xsec-xsec0)/60 write(*,*)' Elapsed time to run chain(min) = ',ttaken write(1,*)' Elapsed time to run chain(min) = ',ttaken end !----------------------------------------------------------------------- subroutine finalze ! at end of run the "tape file" is read and tallies incremented use GMCMC_MOD, only : burnfrac,dataid,idorun,iv,iv1good 1 ,irecmax,irun,nburn,irunacpt,ngrps,nhyp,nrun,nsumlike,ntape 1 ,nthin,nv,nvgood,proddthet,prodsum,slow,sumlike 1 ,tapefile,thetu,thetv,xsec0,One,esumref 1 ,incrhyp,esumrefp,esumu,ihypu,ihypv,c2f implicit none integer time,iend real*8 xlike if(idorun.ge.1)then write(*,*)'FINALZE: Number of moves=',irunacpt 1 ,' Average acceptance fraction (excluding burnfrac)=' 1 ,One*irunacpt/(nrun-nburn) write(1,*)'FINALZE: Number of moves=',irunacpt 1 ,' Average acceptance fraction (excluding burnfrac)=' 1 ,One*irunacpt/(nrun-nburn) if(irunacpt.lt.ntape) then write(*,*)'WARNING: NMOVES < NTAPE, increase NRUN' write(1,*)'WARNING: NMOVES < NTAPE, increase NRUN' endif if(ngrps.eq.1) then prodsum=sumlike*proddthet/nsumlike write(*,*)'FINALZE: Direct estimate of integral of' 1 ,' likelihood over prior=',prodsum 1 ,' for',nsumlike,' trials' write(*,*)' Equivalent Chisq/nData=',-c2f*log(prodsum) write(1,*)'FINALZE: Direct estimate of integral of' 1 ,' likelihood over prior=',prodsum 1 ,' for',nsumlike,' trials' write(1,*)' Equivalent Chisq/nData=',-c2f*log(prodsum) endif endif nv=1 iv=1 nvgood=1 iv1good(1)=1 nburn=nint(burnfrac*nrun) c write(*,*)'FINALZE: nrun,ntape,burnfrac=',nrun,ntape,burnfrac xsec0=time() ! read and process events from tape write(*,*)'FINALZE: reading from ',trim(tapefile) write(1,*)'FINALZE: reading from ',trim(tapefile) open(168,file=trim(tapefile),action='read') call zero0 call zero slow=.true. ! first record contains nthin,ntape,nrun read(168,*)nthin,ntape,nrun if(ntape.eq.0) then write(*,*)'no records on tape' write(1,*)'no records on tape' stop endif ! second record contains dataid read(168,*) ! read thru to find irecmax irecmax(1)=0 irecmax(2)=0 100 read(168,*,end=101,err=101)irun c write(*,*)'irun,nburn=',irun,nburn irecmax(1)=irecmax(1)+1 if(irun.gt.nburn) irecmax(2)=irecmax(2)+1 go to 100 101 continue c write(*,*)'irecmax=',irecmax rewind(168) read(168,*) read(168,*) 1 call rd1tov(168,iend) c write(*,*)'irun=',irun if(iend.eq.1) go to 2 c write(*,*)'before calcvx' if(nhyp.eq.2) then if(ihypv.eq.0) then esumu=esumrefp ihypu=ihypv else call calcvx call move ! move information from v to u variables endif else call calcvx call move endif call increm0 ! increm operates on u variables if(irun.gt.nburn) call increm go to 1 2 continue if(nhyp.eq.2) then xlike=(exp(-esumref)*incrhyp(2,2))/incrhyp(1,2) write(*,*)'FINALZE: Estimate of integral of likelihood' 1 ,' over prior from exp(-esumref)*nNormalHyp/nRefHyp' write(*,*)' =',xlike,' Equivalent Chisq/nData=',-c2f*log(xlike) write(1,*)'FINALZE: Estimate of integral of likelihood' 1 ,' over prior from exp(-esumref)*nNormalHyp/nRefHyp' write(1,*)' =',xlike,' Equivalent Chisq/nData=',-c2f*log(xlike) endif write(*,*)'--------------------end GMCMC----------------------' write(1,*)'--------------------end GMCMC----------------------' close(168) close(1) end !----------------------------------------------------------------------- subroutine run1 ! generate candidate points in parameter space ! test for acceptance, if accepted then move to new point use GMCMC_MOD, only : accept,icalcvx,irun,irunacpt,irundbug 1 ,ndebug,nburn implicit none irun=irun+1 if(ndebug.gt.0.and.mod(irun,irundbug).eq.0) 1 write(*,*)'RUN1:irun=',irun ! generate candidates call gencand c write(*,*)'RUN1: before calcacpt' call calcacpt if(accept) then if(ndebug.gt.0.and.mod(irun,irundbug).eq.0) 1 write(*,*)'RUN1: moving' call move0 if(icalcvx.eq.1) then call move1 call movex endif if(irun.gt.nburn) irunacpt=irunacpt+1 c write(*,*)'RUN1:accept=.true.,nburn=',nburn endif return end !----------------------------------------------------------------------- subroutine move0 ! move basic variable from v to u use GMCMC_MOD, only : esumu,esumv,ihypu,ihypv,iv implicit none ihypu=ihypv esumu=esumv(iv) c write(*,*)'MOVE0: iv,ihypu,esumu=',iv,ihypu,esumu return end !----------------------------------------------------------------------- subroutine move1 ! move theta variables from v to u use GMCMC_MOD, only : ifocus,iiprm,iv,nprm,nprmgrp 1 ,slow,thetu,thetv implicit none integer igrp,iprm,iprmgrp if(slow) then thetu=thetv(:,iv) else igrp=ifocus do iprmgrp=1,nprmgrp(igrp) iprm=iiprm(iprmgrp,igrp) c write(*,*)'igrp,iprmgrp,iprm=',igrp,iprmgrp,iprm thetu(iprm)=thetv(iprm,iv) c write(*,*)'MOVE: iprm,thetu=',iprm,thetu(iprm) enddo endif end !----------------------------------------------------------------------- subroutine move call move0 call move1 call movex return end !----------------------------------------------------------------------- subroutine gencand ! generate candidates for next point in parameter space denoted by v ! current point is denoted by u ! theta variables are generated ! If there are groups of parameters, candidates generated ! only for group ifocus use GMCMC_MOD, only : bad,dthet,esumrefp,esumu,esumv 1 ,fracrw,ialg,icalcvx,icycle,ifocus 1 ,ihypu,ihypv,iiprm,irun,irundbug,iv1good,iwrap 1 ,ndebug,ngrps,nhyp,nprm,nprmgrp,nv,nvgood 1 ,pfocusc,proddthet,pswitch,ranwalk,slow,thetu,thetv implicit none logical nogrps data nogrps /.true./ integer icalcd,igrp,igrp1,igrp2,ih,ihyp,iprm,iprmgrp,ivgood 1 ,iv1 real*8 xrand,dranr250 if(ndebug.gt.0.and.mod(irun,irundbug).eq.0) then if(nhyp.eq.2) then write(*,*)'GENCAND: nv=',nv,' ihypu=',ihypu else write(*,*)'GENCAND: nv=',nv endif endif ! restore starting v values to u values for igrp = ifocus ! which was moved last iteration ! if no v values changed last iteration this is not neccessary ! and icalcvx=0 c write(*,*)'GENCAND: slow,nogrps,nv,nprm=' c 1 ,slow,nogrps,nv,nprm if(slow.or.nogrps) then ! slow but sure, restore everthing ! nogrps=.true. everything changed last time do ivgood=1,nv thetv(:,ivgood)=thetu ! set to u value enddo elseif(icalcvx.eq.1) then do ivgood=1,nv ! only parms for ifocus changed do iprmgrp=1,nprmgrp(ifocus) iprm=iiprm(iprmgrp,ifocus) ! set to u value thetv(iprm,ivgood)=thetu(iprm) c write(*,*)'GENCAND: ivgood,iprm,thetv=',ivgood,iprm,thetu(iprm) enddo ! restore dependent vx variables (extra v variables) call mvbackvx(ivgood) enddo endif c write(*,*)'thetv=',thetv(:,1) c write(*,*)'thetu=',thetu if(irun.lt.1.or.ndebug.eq.0) go to 100 ! check restore c write(*,*)'check restore, irun=',irun do 90 ivgood=1,nv do iprm=1,nprm ! compare to u value if(thetv(iprm,ivgood).ne.thetu(iprm)) then write(*,*)'GENCAND: restore problem irun=',irun write(*,*)'thetv=',thetv(:,ivgood) write(*,*)'thetu=',thetu go to 90 endif enddo 90 continue 100 continue ! end check ! Is random walk to be used? ! use random walk for fraction fracrw of time ! otherwise move variables over entire range ! ranwalk also used in routine CALCPRBD if(fracrw.eq.0.) then ranwalk=.false. elseif(fracrw.eq.1.) then ranwalk=.true. else xrand=dranr250(0) c write(*,*)'xrand=',xrand if(xrand.lt.fracrw) then ranwalk=.true. else ranwalk=.false. endif endif ! change ifocus if grouping used if(ngrps.gt.1) then ! parameter grouping used, calculate ifocus if(icycle.eq.1) then ! cycle ifocus: 1, 2, 3 etc ifocus=mod(ifocus,ngrps)+1 elseif(icycle.eq.2) then ! choose ifocus randomly with uniform probability xrand=dranr250(0) ifocus=ngrps*xrand+1 else ! choose ifocus randomly with some specified probability distribution ! requires pfocusc, the cumulative probability of igrp xrand=dranr250(0) ifocus=icalcd(ngrps,pfocusc(1),xrand) endif igrp1=ifocus igrp2=ifocus else igrp1=1 igrp2=ngrps endif nogrps=.false. ! nogrps used in restore operation at beginning of this subroutine ! normally=.false., except for nhyp=2 when can switch ! from grouping with some number ngrps to no grouping icalcvx=1 ! icalcvx = 1 means that new energy calculation will be done ! icalcvx sometimes set to 0 ! icalcvx=1 indicates that thetv's will be moved and need to be restored ! at beginning of next iteration ih=ihypu+nhyp-1 ! ih = 1, Ref hyp, ihypu=0, ih = 2, normal hyp ! nhyp = 1, no alternate hypothesis ! nhyp = 2 alt hyp with no parameters if(nhyp.eq.1) goto 1 ! come here for nhyp=2 ! Ref hypothesis, ihypu = 0, has no parameters, esum=esumref ! Normal hypothesis, ihypu = 1 if(ihypu.eq.1) then ! normal hypothesis if(dranr250(0).lt.pswitch) then ! probability of switching from normal hyp to ref hyp is pswitch ihypv=0 ! reference hypothesis nvgood=1 iv1good(1)=1 esumv(1)=esumrefp ! corresponds to -log(P0/pswitch) c write(*,*)'GENCAND: ihypv=0,esumv,ihypu,esumu=' c 1 ,esumv(1),ihypu,esumu icalcvx=0 ! no thetv's moved go to 2 else ! ihyp1 -> ihyp1 normal chain movement ihypv=1 ! normal hypothesis go to 1 endif else ! ihypu=0 ref hyp immediately switches ihyp0 -> ihyp1 ihypv=1 ! normal hypothesis ranwalk=.false. ! transitions to entire space ! ranwalk also used in routine CALCPRBD nogrps=.true. ! no grouping used, important for next restore igrp1=1 igrp2=ngrps go to 1 endif 1 continue ! calculate theta's for candidates nvgood=0 do 1000 iv1=1,nv bad(iv1)=.false. ! start with clean badness slate do igrp=igrp1,igrp2 do iprmgrp=1,nprmgrp(igrp) iprm=iiprm(iprmgrp,igrp) ! move theta possibly with random walk ! or possibly chosing from entire space, depending on RANWALK if(ranwalk) then if(dthet(iprm).ge.1.) then thetv(iprm,iv1)=dranr250(0) else thetv(iprm,iv1)=thetu(iprm) 1 +(2*dranr250(0)-1)*dthet(iprm) if(iwrap(iprm).eq.2) then c write(*,*)'before wrap=',thetv(iprm,iv1) thetv(iprm,iv1)=thetv(iprm,iv1)+10 1 -int(thetv(iprm,iv1)+10) c write(*,*)'after wrap=',thetv(iprm,iv1) elseif(thetv(iprm,iv1).lt.0.or. 1 thetv(iprm,iv1).gt.1) then bad(iv1)=.true. c write(*,*)'iprm,thetv=',iprm,thetv(iprm,iv1) if(nv.eq.1) go to 999 go to 1000 endif endif else thetv(iprm,iv1)=dranr250(0) endif enddo enddo c write(*,*)'iv1,iprm,thetv,thetu,dthet,deviation=',iv1,iprm c 1 ,thetv(iprm,iv1),thetu(iprm),dthet(iprm) c 1 ,(thetv(iprm,iv1)-thetu(iprm))/dthet(iprm) nvgood=nvgood+1 iv1good(nvgood)=iv1 1000 continue 2 continue ! skip normal candidate calculation above 999 continue if(ndebug.gt.0.and.mod(irun,irundbug).eq.0) then write(*,*)'iv1good iv1 ihyp igrp iprm thetv' 1 ,' dthet deviation ranwalk bad' ihyp=ihypv do ivgood=1,nvgood iv1=iv1good(ivgood) c write(*,*)'iv1,ngrps,igrp,nprmgrp=' c 1 ,iv1,ngrps,ifocus,nprmgrp(igrp) igrp=ifocus do igrp=igrp1,igrp2 do iprmgrp=1,nprmgrp(igrp) iprm=iiprm(iprmgrp,igrp) write(*,'(5i8,2f12.8,g12.4,2l10)') 1 iv1,iv1,ihyp,igrp,iprm 1 ,thetv(iprm,iv1),dthet(iprm) 1 ,(thetv(iprm,iv1)-thetu(iprm))/dthet(iprm) 1 ,ranwalk,bad(iv1) enddo enddo enddo c write(*,*)'GENCAND: before call calcvx' endif if(icalcvx.eq.1) call calcvx end !----------------------------------------------------------------------- subroutine calcacpt ! call calcprbd to calculate the unnormalized probability prbnd of ! iv1good=1,nvgood new configurations ! including the current configuration as iv1good=0 ! These are stored in prbcd(i), i=1,nvgood+1 ! Generate ivgoodnew from prbcd, which determines whether ! to stay (ivgoodnew=0) ! or move to a new configuration use GMCMC_MOD, only : accept,ihypu,ihypv,irun,irundbug 1 ,iv,iv1good,ndebug,nprm,nvgood,prbcd,thetv implicit none integer i,igoodnew,ihyp,ip,iprm,nprb real*8 xrand,dranr250 ! ivgood=1,nvgood are new configurations with bad.ne..true. if(nvgood.eq.0) then c write(*,*)'CALCACPT: nvgood=',nvgood igoodnew=0 else call calcprbd nprb=nvgood+1 ! igoodnew=igenpd(nprb,prbcd)-1 xrand=dranr250(0) c write(*,*)'IGENPD:xrand=',xrand do ip=1,nprb if(xrand.le.prbcd(ip)) go to 10 enddo ip=nprb write(*,*)'CALCACPT ERROR: should not come here' 1 ,', irun=',irun,', xrand,prbcd=',xrand,(prbcd(i),i=1,nprb) 10 continue igoodnew=ip-1 endif if(igoodnew.gt.0) then accept=.true. iv=iv1good(igoodnew) else accept=.false. iv=0 endif if(ndebug.gt.0.and.mod(irun,irundbug).eq.0) then if(iv.eq.0) then ihyp=ihypu else ihyp=ihypv endif write(*,*)'CALCACPT: accept=',accept,' iv=',iv 1 ,' ihyp=',ihyp if(accept) then write(*,*)' iprm thetv' do 12 iprm=1,nprm write(*,'(i8,f12.8)')iprm,thetv(iprm,iv) 12 continue endif endif end !----------------------------------------------------------------------- subroutine calcprbd ! calculates relative probability of new points in parameter space use GMCMC_MOD, only : bad,eminv,esumu,esumv 1 ,ialg,ifocus,ihypu,ihypv,iprbmax,irun,irundbug,iv1good,nburn 1 ,ndebug,ngrps,nhyp,nsumlike ,nvgood,prbcd,prbnd,pswitch 1 ,sumlike,ranwalk implicit none integer igood,iiv,iprb,ivgood,ivgoodmin,ivmax,ivmin,iv1 real*8 eminv0,esum,sum,sumc logical dosum ivmin=iv1good(1) ivgoodmin=1 ! find minimum energy of new candidates do ivgood=1,nvgood iv1=iv1good(ivgood) if(esumv(iv1).lt.esumv(ivmin)) then c write(*,*)'iv1,esumv(iv1),esumv(ivmin)=' c 1 ,iv1,esumv(iv1),esumv(ivmin) ivgoodmin=ivgood ivmin=iv1 endif enddo eminv=esumv(ivmin) if(esumu.lt.eminv) then iprbmax=1 eminv0=esumu else iprbmax=ivgoodmin+1 eminv0=eminv endif sum=0 do iprb=2,nvgood+1 ivgood=iprb-1 iv1=iv1good(ivgood) if(ialg.eq.1) then if(esumv(iv1).lt.esumu) then prbnd(iprb)=1.d0/nvgood else prbnd(iprb)=dexp(esumu-esumv(iv1))/nvgood endif elseif(ialg.eq.2) then prbnd(iprb)=1/((1+dexp((esumv(iv1)-esumu)))*nvgood) else prbnd(iprb)=dexp(eminv0-esumv(iv1)) endif sum=sum+prbnd(iprb) enddo if(ialg.eq.3) then ! multiple candidate algorithm prbnd(1)=dexp(eminv0-esumu) sum=sum+prbnd(1) do iprb=1,nvgood+1 prbnd(iprb)=prbnd(iprb)/sum enddo else ! calculate probability for iprb=1, chain stays put prbnd(1)=1-sum endif sumc=0 dosum=ranwalk.and.ngrps.eq.1.and.ihypv.eq.1 1 .and.irun.gt.nburn do iprb=1,nvgood+1 sumc=sumc+prbnd(iprb) ! calculate the cummulative probability over the candidates prbcd(iprb)=sumc ! This code is used to calculate sumlike, which can be used ! to directly estimate the normalization integral ! only works for ranwalk=.true. and ngrps=1 if(dosum.and.iprb.gt.1) then ivgood=iprb-1 iv1=iv1good(ivgood) esum=esumv(iv1) sumlike=sumlike+dexp(-esum) nsumlike=nsumlike+1 endif enddo if(ndebug.gt.0.and.mod(irun,irundbug).eq.0) then if(iprbmax.eq.1) then ivmax=0 else ivmax=iv1good(iprbmax-1) endif write(*,*)'CALCPRBD: esumu,eminv,ivMaxProb=' 1 ,esumu,eminv,ivmax write(*,*)' iv igood esum prb' do iprb=1,nvgood+1 if(iprb.eq.1) then esum=esumu iiv=0 igood=0 else igood=iprb-1 iiv=iv1good(igood) esum=esumv(iiv) endif write(*,'(2i8,g12.3,f12.8)')iiv,igood,esum,prbnd(iprb) enddo endif end !----------------------------------------------------------------------- function igenpd(n,pc) ! generate a random integer from 1 to n ! with probability given by the cumulative probability pc implicit none integer i,igenpd,n real*8 xrand,dranr250 real*8 pc(n) xrand=dranr250(0) c write(*,*)'IGENPD:xrand=',xrand do igenpd=1,n if(xrand.le.pc(igenpd)) return enddo igenpd=n write(*,*)'IGENPD ERROR: should not come here' 1 ,', xrand,pc=',xrand,(pc(i),i=1,n) end !----------------------------------------------------------------------- function icalcd(n,pc,theta) ! given cumulative probability distribution pc(i) over i=1,n ! calculate the integer icalcd that has cumulative probability theta implicit none integer icalcd,n real*8 theta,pc(n) do icalcd=1,n if(theta.le.pc(icalcd)) return enddo icalcd=n write(*,'('' ICALCD ERROR: should not come here'' 1 ,'' theta,pc(n)='',2z17.16)') theta,pc(n) write(1,'('' ICALCD ERROR: should not come here'' 1 ,'' theta,pc(n)='',2z17.16)') theta,pc(n) end !----------------------------------------------------------------------- subroutine wrt1u(iunit) ! write one u record on "tape" use GMCMC_MOD, only:cesum,cthet,esumu,ifocus 1 ,ihypu,irun,nprm,thetu implicit none integer iunit write(cesum,'(z16)') esumu write(cthet,'(z16)') thetu write(iunit,*)irun 1 ,ihypu,cesum,ifocus,cthet end !----------------------------------------------------------------------- subroutine rd1tou(iunit,iend) ! read one record from "tape" to u arrays use GMCMC_MOD, only:cesum,cthet,esumu,ifocus,ihypu,irun 1 ,ngrps,nprm,thetu implicit none integer iend,iunit iend=0 read(iunit,*,end=99)irun 1 ,ihypu,cesum,ifocus,cthet read(cesum,'(z16)') esumu read(cthet,'(z16)') thetu if(ngrps.eq.1) ifocus=1 return 99 continue iend=1 end !----------------------------------------------------------------------- subroutine rd1tov(iunit,iend) ! read one record from "tape" to v arrays use GMCMC_MOD, only:cesum,cthet,esumv,ifocus,ihypv,irun 1 ,iv,ngrps,nprm,thetv implicit none integer iend,iunit iend=0 read(iunit,*,end=99,err=99)irun 1 ,ihypv,cesum,ifocus,cthet read(cesum,'(z16)') esumv(iv) read(cthet,'(z16)') thetv(:,iv) if(ngrps.eq.1) ifocus=1 c write(*,*)'RD1TOV: iv,ihypv,thetv(1,iv)=' c 1 ,iv,ihypv,thetv(1,iv) return 99 continue iend=1 end !----------------------------------------------------------------------- logical function chkfile(filename,info) ! check if a file with path filename of length >0 exists. implicit none integer i character*(*) filename,info character*120 cmd ! use dir command to get file info, put in file FILEINFO.TXT cmd='dir '//trim(filename)//' > fileinfo.txt' c write(*,*)'cmd=',cmd call shellcmd(cmd) ! if file exists, file info is in 6th record c write(*,*)'before open' open(221,file='fileinfo.txt') c write(*,*)'before do i=1,6' do i=1,6 read(221,'(a39)',end=99,err=99) info c write(*,*)'i,info=',i,info enddo close(221) chkfile=.true. ! check for 0 file size if(info(37:38).eq.' 0') chkfile=.false. c write(*,*)'chkfile,info=',chkfile,info return 99 continue chkfile=.false. c write(*,*)'ISTHERE: chkfile,info=',chkfile,info close(221) end !----------------------------------------------------------------------- subroutine shellcmd(cmd) implicit none character*(*) cmd integer system,iexit cmd='cmd /c '//trim(cmd) iexit=system(trim(cmd)) c write(*,*)'iexit=',iexit if(iexit.ne.0.and.iexit.ne.1) then write(*,*)'SHELLCMD: abnormal ending, iexit=' 1 ,iexit,', cmd=',cmd write(1,*)'SHELLCMD: abnormal ending, iexit=' 1 ,iexit,', cmd=',cmd endif return end !------------------------------Begin r250------------------------------- real*8 function dranr250(ir) if(ir.gt.0) call r250init(ir) dranr250=irand250()/2147483648.0 ! largest 31 bit integer + 1 = 2**31 return end !----------------------------------------------------------------------- subroutine rtest ! Test the random number generators Print *,'The EXPECTED values are in parenthesis' ! First test the LCM generator ix = lcmrand( 1 ) Do 10 I = 1, 9999 ix = lcmrand( -1 ) 10 EndDo Print *,'Final LCM value: (1043618065)', ix ! Now initialize and test R250 Call R250Init( 1 ) Do 20 I = 1, 10000 ix = irand250() 20 EndDo Print *,'Final R250 value: (267771767)', ix END !----------------------------------------------------------------------- ! R250.F77 The R250 Pseudo-random number generator ! algorithm from: ! Kirkpatrick, S., and E. Stoll, 1981; A Very Fast Shift-Register ! Sequence Random Number Generator, Journal of Computational Physics, ! V. 40. p. 517 ! see also: ! Maier, W.L., 1991; A Fast Pseudo Random Number Generator, ! Dr. Dobb's Journal, May, pp. 152 - 157 ! Uses the Linear Congruential Method, ! the "minimal standard generator" ! Park & Miller, 1988, Comm of the ACM, 31(10), pp. 1192-1201 ! for initialization ! For a review of BOTH of these generators, see: ! Carter, E.F, 1994; Generation and Application of Random Numbers, ! Forth Dimensions, Vol. XVI, Numbers 1,2 May/June, July/August ! $Author: skip $ ! $Workfile: r250.f $ ! $Revision: 1.1 $ ! $Date: 07 Nov 1996 01:23:06 $ ! =================================================================== Function lcmrand(ix) ! The minimal standard PRNG for 31 bit unsigned integers ! designed with automati! overflow protection ! uses ix as the seed value if it is greater than zero ! otherwise it is ignored Integer*4 ix Integer*4 a, b, m, q, r Integer*4 hi, lo, test Integer*4 x SAVE x Parameter (a = 16807, b = 0, m = 2147483647) Parameter (q = 127773, r = 2836) If ( ix .gt. 0 ) x = ix hi = x / q lo = mod( x, q ) test = a * lo - r * hi if ( test .gt. 0 ) then x = test else x = test + m endif lcmrand = x return End ! =================================================================== C ! R250, call R250Init with the desired initial seed BEFORE ! the first invocation of IRAND250() C ! =================================================================== Subroutine R250Init(iseed) Integer k, mask, msb Integer indexf, indexb, buffer(250) Common/R250COM/indexf,indexb,buffer Integer ms_bit, all_bits, half_range, step Parameter ( ms_bit = Z'40000000') Parameter ( half_range = Z'20000000' ) Parameter ( all_bits = Z'7FFFFFFF' ) Parameter ( step = 7 ) indexf = 1 indexb = 104 k = iseed Do 10 i = 1, 250 buffer(i) = lcmrand( k ) k = -1 10 EndDo Do 20 i = 1, 250 if ( lcmrand( -1 ) .gt. half_range ) then buffer(i) = ior( buffer(i), ms_bit ) endif 20 EndDo msb = ms_bit mask = all_bits Do 30 i = 0,30 k = step * i + 4 buffer(k) = iand( buffer(k), mask ) buffer(k) = ior( buffer(k), msb ) msb = msb / 2 mask = mask / 2 30 EndDo Return END !----------------------------------------------------------------------- Function irand250() ! R250 PRNG, run after R250_Init Integer newrand Integer indexf, indexb, buffer(250) Common/R250COM/indexf, indexb,buffer newrand = ieor( buffer(indexf), buffer(indexb) ) buffer(indexf) = newrand indexf = indexf + 1 if ( indexf .gt. 250 ) indexf = 1 indexb = indexb + 1 if ( indexb .gt. 250 ) indexb = 1 irand250= newrand return End !----------------------------end r250----------------------------------- !###############################END GMCMC############################### !----------------------------------------------------------------------- !########################test program follows########################### !----------------------------------------------------------------------- module testMCMC_mod integer, save :: ndim,nsum real*8, save :: x0,sig0 real*8, save, dimension(:), allocatable :: xu real*8, save, dimension(:,:), allocatable :: xv,xsum end module testMCMC_mod !----------------------------------------------------------------------- program testMCMC ! test of MCMC basic programs ! Markov Chain Monte Carlo, writing tape file of events ! In the core MCMC prog, all parameters are expressed ! in terms of theta variables, which run from 0 to 1 ! All priors on theta variables are independent and uniform ! The actual variables ("x" variables) are expressed as ! functions of the theta variables, where the function is such that ! actual variable has the desired prior probability distribution when ! theta variable is unifomly distributed from 0 to 1 ! set up call to gmcmc use GMCMC_MOD, only: xsec,xsec0 use testMCMC_mod implicit none integer ialg,icycle,igrp,iseed,j,iprmgrp,iprm,iwrap 1 ,mdata,nburn,ndebug,ngrps,nprm 1 ,nprmsgrp,nrpts,ntape,nv,time real*8 dthetrw,fracrw,thet0,xliketrue,avg integer, dimension(:), allocatable :: iwrap1,nprmgrp integer, dimension(:,:), allocatable :: iiprm real*8, dimension(:), allocatable:: pfocus,thet real*8, dimension(:,:), allocatable:: dthet real*8 burnfrac character*70 datafile data datafile/' '/ character*10 tapefile logical slow,exist data slow/.false./ integer*8 nrun namelist/l/iseed,thet0,nrun,ntape,nrpts,ndebug,burnfrac 1 ,fracrw,dthetrw,iwrap,icycle,ngrps,nv,ialg,ndim,x0,sig0 tapefile='tape.txt' ! name of tape file ntape=100 ! number of iterations written on tape file ialg=1 ! MCMC algorithm: 1 for MRT, 2 for B, 3 for MC nv=1 ! number of multiple candidates nrun=1000000 ! number of chain iterations in MCMC run iwrap=1 ! wrap variable if random walk goes outside of (0,1) range icycle=1 ! cycling of parameter groups that are varied: 1 for cyclic, 2 random, ! 3 specified probability supplied in pfocus ngrps=6 ! number of parameter groups thet0=0 ! starting theta value (for all dimensions) burnfrac=.05d0 ! initial discarded fraction of chain iterations dthetrw=.001d0 ! random walk Delta theta fracrw=0.5d0 ! fraction of time random walk is used rather than generating candidates over entire space nrpts=10 ! number of reports during course of run iseed=1 ! random number generator seed ndebug=0 ! >0 for debug printout ! Specification of desired steady state distribution: ndim=6 ! dimensionality x0=0.5 ! location of mode(s) sig0=1.d-6 ! width of mode(s) inquire(file='gmcmc.in',exist=exist) if(.not.exist)goto 998 open(10,file='gmcmc.in') read(10,l,err=999) close(10) 998 continue write(*,l) nprm=ndim nprmsgrp=ndim allocate (iiprm(nprmsgrp,ngrps),nprmgrp(ngrps),iwrap1(nprm)) allocate (pfocus(nprm),thet(nprm),dthet(nprm,2)) allocate (xv(ndim,nv),xu(ndim),xsum(3,ndim)) if(ngrps.eq.1) then nprmgrp(1)=ndim igrp=1 pfocus(igrp)=1 do iprmgrp=1,nprmgrp(igrp) iiprm(iprmgrp,igrp)=iprmgrp enddo elseif(ngrps.eq.ndim) then do igrp=1,ngrps pfocus(igrp)=1 nprmgrp(igrp)=1 iiprm(1,igrp)=igrp enddo else write(*,*)'bad ngrps=',ngrps stop endif nburn=nint(burnfrac*nrun) do igrp=1,ngrps do iprmgrp=1,nprmgrp(igrp) iprm=iiprm(iprmgrp,igrp) iwrap1(iprm)=iwrap dthet(iprm,1)=dthetrw dthet(iprm,2)=dthet(iprm,1) thet(iprm)=thet0 enddo enddo xliketrue=1 do j=1,ndim xliketrue=xliketrue*sqrt(2*3.14159256)*sig0 enddo write(*,*)'integral of likelihood function over prior=' 1 ,xliketrue mdata=ndim call GMCMC(nrun,ntape,nrpts,ndebug 1 ,mdata,iseed,icycle,nv,ialg 1 ,slow,burnfrac,fracrw 1 ,nprm,ngrps,nprmsgrp,nprmgrp,iiprm 1 ,pfocus,iwrap1,thet,dthet 1 ,datafile) write(*,*)'After GMCMC: nsum=',nsum write(*,*)'j xavg xsd' do j=1,ndim avg=xsum(1,j)/nsum write(*,*)j,avg,sqrt(xsum(2,j)/nsum-avg**2) enddo write(*,*) xsec=time() write(*,*)'Elapsed time for post processing (sec) =' c ,(xsec-xsec0) stop 999 write(*,*)'ERROR involving namelist read of GMCMC.IN' go to 998 end !----------------------------------------------------------------------- subroutine zero use testMCMC_mod, only:xsum,nsum implicit none ! zero tallies xsum=0 nsum=0 end !----------------------------------------------------------------------- subroutine increm ! increment talleys use testMCMC_mod, only:xsum,xu,ndim,nsum implicit none integer j do j=1,ndim xsum(1,j)=xsum(1,j)+xu(j) xsum(2,j)=xsum(2,j)+xu(j)**2 enddo nsum=nsum+1 end !----------------------------------------------------------------------- subroutine calcvx ! calculate dependent v variables from theta variables use GMCMC_MOD, only:esumv,iv1good,nvgood,thetv use testMCMC_mod, only:ndim,xv,x0,sig0 implicit none integer iv1,ivgood,iprm,j real*8 chi ! here is the calculation of the likelihood function ! normally done in two possible ways, slow but sure and as fast as possible ! but hopefully also sure! ! This simulates an energy calculation in terms of the actual parameters ! rather than the theta parameters do ivgood=1,nvgood iv1=iv1good(ivgood) esumv(iv1)=0 do j=1,ndim iprm=j xv(j,iv1)=thetv(iprm,iv1) ! normally there is some general functional relationship between the problem ! variables xv and the MCMC theta variables (theta from 0 to 1 with uniform prior probability) chi=(xv(j,iv1)-x0)/sig0 ! normally some general functional relationship between problem variables and chi esumv(iv1)=esumv(iv1)+chi**2/2 enddo enddo end !----------------------------------------------------------------------- subroutine mvbackvx(iv1) ! restore extra dependent v variables use testMCMC_mod, only:ndim,xu,xv implicit none integer iv1,j do j=1,ndim xv(j,iv1)=xu(j) enddo end !----------------------------------------------------------------------- subroutine movex ! If a candidate Xv is accepted, chain is moved and Xv->Xv ! This routine moves v variables to u variables ! if ngrps > 0 for the group ifocus only, the minimum required ! after gencand, calcacpt, because other variables have not changed use GMCMC_MOD, only: iv use testMCMC_mod, only:ndim,xu,xv implicit none integer j ! move other dependent variables do j=1,ndim xu(j)=xv(j,iv) enddo end !-----------------------------------------------------------------------