diff --git a/README b/README new file mode 100644 index 0000000..a220f52 --- /dev/null +++ b/README @@ -0,0 +1,76 @@ +********************************************************************** +Adds trans-Planckian like modulation to primordial power spectra +expanded at first order in Hubble and sound flow parameters (valid for +k-inflation). A switch allows the inclusion of second order +corrections but in the Hubble flow functions only (valid for +canonically normalised fields). +********************************************************************** + +-These patches modify the cosmomc or camb codes and therefore may have +bugs that are not present in the originals. The authors of the +original versions should not be bothered with problems resulting in +the use of these patches. Feel free to insult me, or at least contact +me, in case of bugs! + +-Apply the patch in the cosmomc/ directory: +patch -p1 < purpose_MMYY.patch + +*Modifications traced by the comment "wig" in the source files. + +*Please note that the genuine Makefiles have been heavily changed and +may require edition to work on your system. They assume the existence +of a directory "cosmomc/WMAP" which points to the WMAP team likelihood +code directory. The WMAP data are also assumed to reside in the +"cosmomc/data" directory. + + +*Modified files and settings + +- camb/cmbmain.f90 +New switches "kSrcSampleBoost" and "kIntSampleBoost" to check accuracy. + +- camb/modules.f90 +New switches "DoClsInterpol" and "ComputeEveryl". AccuracyBoost and +lAccuracyBoost default values modified. points=10000 to accurately +sample oscillatory matter power spectra. + +- camb/power_sr.f90 +Implementation of Slow-Roll/k-inflation power spectrum and its +trans-Planckian extension. + +- camb/inidriver.F90 +New keys for the above material. + +- camb/params.ini +New keys for the above material. + +- cosmomc/CMB_Cls_simple.f90 +Added screen output to check camb accuracy default values. New +switches kSrcSampleLevel, kIntSampleLevel, CambComputeEveryl to +control the camb accuracy from cosmomc. + +- cosmomc/cmbtypes.f90 +Default value for "write_all_cls" to true. + +- cosmomc/params_CMB.f90 +Take care of the new power spectra parameters + +- cosmomc/settings.f90 +num_init_power modified + +- cosmomc/MCMC.f90 +Add MpiRank output in the logs + +- cosmomc/GetDist.f90 +3d transparent plots of the mean likelihoods. Try colormap('copper'). + +- cosmomc/driver.F90 +2nd instance in mpirun shifts the file names. + +- cosmomc/params.ini +New keys for the above material. + +- cosmomc/distparams.ini +New keys for the above material. + + diff --git a/camb/Makefile b/camb/Makefile index 5c09a67..65d2274 100644 --- a/camb/Makefile +++ b/camb/Makefile @@ -1,83 +1,115 @@ -#CAMB Makefile +# >>> DESIGNED FOR GMAKE <<< +#CAMB System unified Makefile + +ext=$(shell uname | cut -c1-3) +ifeq ($(ext),IRI) +F90C=f90 +FFLAGS = -Ofast -mp -n32 -LANG:recursive=ON +LFLAGS= +endif -#Set FISHER=Y to compile bispectrum fisher matrix code -FISHER= +ifeq ($(ext),Lin) +F90C=gfortran +FFLAGS= -O -fopenmp +LFLAGS= +endif -#Edit for your compiler -#Note there are many old ifc versions, some of which behave oddly +ifeq ($(ext),OSF) +F90C=f90 +FFLAGS= -omp -O -arch host -math_library fast -tune host -fpe1 +LFLAGS= +endif +ifeq ($(ext),Sun) +F90C=f90 +FFLAGS= -O4 -xarch=native64 -openmp -ftrap=%none +LFLAGS= +endif -#Intel , -openmp toggles mutli-processor: -#note version 10.0 gives wrong result for lensed when compiled with -openmp [fixed in 10.1] -F90C = ifort -FFLAGS = -openmp -fast -W0 -WB -fpp2 -vec_report0 -ifneq ($(FISHER),) -FFLAGS += -mkl +ifeq ($(ext),AIX) +F90C=xlf90_r +FFLAGS= -O4 -qsmp=omp -qmaxmem=-1 -qstrict -qfree=f90 -qsuffix=f=f90:cpp=F90 +LFLAGS= endif -#Gfortran compiler: -#The options here work in v4.5, delete from RHS in earlier versions (15% slower) -#if pre v4.3 add -D__GFORTRAN__ -#With v4.6+ try -Ofast -march=native -fopenmp -#On my machine v4.5 is about 20% slower than ifort -#F90C = gfortran -#FFLAGS = -O3 -fopenmp -ffast-math -march=native -funroll-loops - - -#Old Intel ifc, add -openmp for multi-processor (some have bugs): -#F90C = ifc -#FFLAGS = -O2 -Vaxlib -ip -W0 -WB -quiet -fpp2 -#some systems can can also add e.g. -tpp7 -xW - -#G95 compiler -#F90C = g95 -#FFLAGS = -O2 - -#SGI, -mp toggles multi-processor. Use -O2 if -Ofast gives problems. -#F90C = f90 -#FFLAGS = -Ofast -mp - -#Digital/Compaq fortran, -omp toggles multi-processor -#F90C = f90 -#FFLAGS = -omp -O4 -arch host -math_library fast -tune host -fpe1 - -#Absoft ProFortran, single processor: -#F90C = f95 -#FFLAGS = -O2 -cpu:athlon -s -lU77 -w -YEXT_NAMES="LCS" -YEXT_SFX="_" - -#NAGF95, single processor: -#F90C = f95 -#FFLAGS = -DNAGF95 -O3 - -#PGF90 -#F90C = pgf90 -#FFLAGS = -O2 -DESCAPEBACKSLASH -Mpreprocess - -#Sun V880 -#F90C = mpf90 -#FFLAGS = -O4 -openmp -ftrap=%none -dalign -DMPI - -#Sun parallel enterprise: -#F90C = f95 -#FFLAGS = -O2 -xarch=native64 -openmp -ftrap=%none -#try removing -openmp if get bus errors. -03, -04 etc are dodgy. - -#IBM XL Fortran, multi-processor (run gmake) -#F90C = xlf90_r -#FFLAGS = -DESCAPEBACKSLASH -DIBMXL -qsmp=omp -qsuffix=f=f90:cpp=F90 -O3 -qstrict -qarch=pwr3 -qtune=pwr3 - -#Settings for building camb_fits -#Location of FITSIO and name of library -FITSDIR = /home/cpac/cpac-tools/lib -FITSLIB = cfitsio -#Location of HEALPIX for building camb_fits -HEALPIXDIR = /home/cpac/cpac-tools/healpix - -ifneq ($(FISHER),) + +#Files containing evolution equations initial power spectrum module +EQUATIONS = equations +POWERSPECTRUM = power_sr +REIONIZATION = reionization +RECOMBINATION = recfast +BISPECTRUM = SeparableBispectrum +DENABLE_FISHER= + +#Module doing non-linear scaling +NONLINEAR = halofit + +#Driver program +DRIVER = inidriver.F90 +#DRIVER = sigma8.f90 +#DRIVER = tester.f90 + + +CAMBLIB = libcamb.a + +ifneq ($(DENABLE_FISHER),) FFLAGS += -DFISHER +LFLAGS += -llapack -lblas EXTCAMBFILES = Matrix_utils.o else EXTCAMBFILES = endif -include ./Makefile_main +HEALPIXDIR= + +#Shouldn't need to change anything else... + +F90FLAGS = $(FFLAGS) +HEALPIXLD = -L$(HEALPIXDIR)/lib -lhealpix + +CAMBOBJ = constants.o utils.o $(EXTCAMBFILES) subroutines.o inifile.o $(POWERSPECTRUM).o $(RECOMBINATION).o $(REIONIZATION).o modules.o \ + bessels.o $(EQUATIONS).o $(NONLINEAR).o lensing.o $(BISPECTRUM).o \ + cmbmain.o camb.o + +default: camb.$(ext) + +all: camb.$(ext) $(CAMBLIB) + +subroutines.o: constants.o utils.o +$(POWERSPECTRUM): subroutines.o inifile.o +$(RECOMBINATION).o: subroutines.o inifile.o +$(REIONIZATION).o: constants.o inifile.o +modules.o: $(REIONIZATION).o $(POWERSPECTRUM).o $(RECOMBINATION).o +bessels.o: modules.o +$(EQUATIONS): bessels.o +$(NONLINEAR).o: modules.o +lensing.o: bessels.o +$(BISPECTRUM).o: lensing.o modules.o +cmbmain.o: lensing.o $(NONLINEAR).o $(EQUATIONS).o $(BISPECTRUM).o +camb.o: cmbmain.o +Matrix_utils.o: utils.o + +camb.$(ext): $(CAMBOBJ) $(DRIVER) + $(F90C) $(F90FLAGS) $(CAMBOBJ) $(DRIVER) $(LFLAGS) -o $@ + +$(CAMBLIB): $(CAMBOBJ) + ar -r $@ $? + +camb_fits: writefits.f90 $(CAMBOBJ) $(DRIVER) + $(F90C) $(F90FLAGS) -I$(HEALPIXDIR)/include tils.o $(CAMBOBJ) writefits.f90 $(DRIVER) $(HEALPIXLD) $(LFLAGS) -DWRITE_FITS -o $@ + +%.o: %.f90 + $(F90C) $(F90FLAGS) -c $*.f90 + +utils.o: + $(F90C) $(F90FLAGS) -c utils.F90 + +$(BISPECTRUM).o: + $(F90C) $(F90FLAGS) -c $(BISPECTRUM).F90 + +Matrix_utils.o: + $(F90C) $(F90FLAGS) -c Matrix_utils.F90 + +clean: + -rm -f *.o *.a *.d core *.mod + diff --git a/camb/cmbmain.f90 b/camb/cmbmain.f90 index 90ce721..821b73c 100644 --- a/camb/cmbmain.f90 +++ b/camb/cmbmain.f90 @@ -657,17 +657,24 @@ contains SourceAccuracyBoost = AccuracyBoost if (CP%WantScalars .and. CP%Reion%Reionization .and. CP%AccuratePolarization) then - dlnk0=2._dl/10/SourceAccuracyBoost +!wig +! dlnk0=2._dl/10/SourceAccuracyBoost + dlnk0=2._dl/10/SourceAccuracyBoost/kSrcSampleBoost !Need this to get accurate low l polarization else - dlnk0=5._dl/10/SourceAccuracyBoost +! dlnk0=5._dl/10/SourceAccuracyBoost + dlnk0=5._dl/10/SourceAccuracyBoost/kSrcSampleBoost +!end wig if (CP%closed) dlnk0=dlnk0/2 end if if (CP%AccurateReionization) dlnk0 = dlnk0/2 - - dkn1=0.6_dl/taurst/SourceAccuracyBoost - dkn2=0.9_dl/taurst/SourceAccuracyBoost +!wig +! dkn1=0.6_dl/taurst/SourceAccuracyBoost +! dkn2=0.9_dl/taurst/SourceAccuracyBoost + dkn1=0.6_dl/taurst/SourceAccuracyBoost/kSrcSampleBoost + dkn2=0.9_dl/taurst/SourceAccuracyBoost/kSrcSampleBoost +!end wig if (HighAccuracyDefault) dkn2=dkn2/1.2 if (CP%WantTensors .or. CP%WantVectors) then dkn1=dkn1 *0.8_dl @@ -1082,8 +1089,10 @@ contains qmax_int = min(qmax,max_bessels_etak/CP%tau0) - - IntSampleBoost=AccuracyBoost +!wig +! IntSampleBoost=AccuracyBoost + IntSampleBoost=AccuracyBoost*kIntSampleBoost +!end wig if (do_bispectrum) then IntSampleBoost = IntSampleBoost * 2 if (hard_bispectrum) IntSampleBoost = IntSampleBoost * 2 diff --git a/camb/inidriver.F90 b/camb/inidriver.F90 index b355761..73102b4 100644 --- a/camb/inidriver.F90 +++ b/camb/inidriver.F90 @@ -282,6 +282,20 @@ end if end if +!wig + + kSrcSampleBoost= Ini_Read_Double('k_src_sample_boost',1.d0) + kIntSampleBoost= Ini_Read_Double('k_int_sample_boost',1.d0) + DoClsInterpol = Ini_Read_Logical('do_cls_interpol',.true.) + ComputeEveryl = Ini_Read_Logical('compute_every_l',.true.) + + if (FeedbackLevel > 0) then + write (*,*)'do_cls_interpol is',DoClsInterpol + write (*,*)'compute_every_l is',ComputeEveryl + write (*,*)'primordial power spectra is: ',Power_Name + endif +!end wig + call Ini_Close if (.not. CAMB_ValidateParams(P)) stop 'Stopped due to parameter error' @@ -293,6 +307,23 @@ if (FeedbackLevel > 0) then Age = CAMB_GetAge(P) write (*,'("Age of universe/GYr = ",f7.3)') Age +!wig + if ((Power_Name.eq.'power_sr').or.(Power_Name.eq.'power_dbi')) then + + do i=1,P%InitPower%nn + write (*,'("Slow-roll epsilon1 = ",f12.8)') P%InitPower%eps1(i) + write (*,'("Slow-roll epsilon2 = ",f12.8)') P%InitPower%eps2(i) + write (*,'("Slow-roll epsilon3 = ",f12.8)') P%InitPower%eps3(i) + write (*,'("DBI delta1 = ",f12.8)') P%InitPower%del1(i) + write (*,'("DBI gamma = ",f12.8)') P%InitPower%gam(i) + write (*,'("TPL amplitude = ",f12.8)') P%InitPower%Xsigma0(i) + write (*,'("TPL scale = ",f12.8)') P%InitPower%loginvsigma0(i) + write (*,'("TPL phase = ",f12.8)') P%InitPower%varphi(i) + enddo + else + stop 'Power_Name unknown' + endif +!end wig end if if (global_error_flag==0) call CAMB_GetResults(P) diff --git a/camb/modules.f90 b/camb/modules.f90 index 6053b4c..7764896 100644 --- a/camb/modules.f90 +++ b/camb/modules.f90 @@ -38,8 +38,10 @@ character(LEN=*), parameter :: version = 'Jan_12' integer :: FeedbackLevel = 0 !if >0 print out useful information about the model - - logical, parameter :: DebugMsgs=.false. !Set to true to view progress and timing +!wig +! logical, parameter :: DebugMsgs=.false. !Set to true to view progress and timing + logical, parameter :: DebugMsgs=.false. +!end wig logical, parameter :: DebugEvolution = .false. !Set to true to do all the evolution for all k @@ -187,15 +189,17 @@ real(dl) :: lSampleBoost=1._dl !Increase lSampleBoost to increase sampling in lSamp%l for Cl interpolation - - real(dl) :: AccuracyBoost =1._dl +!wig + real(dl) :: AccuracyBoost =8._dl +!end wig !Decrease step sizes, etc. by this parameter. Useful for checking accuracy. !Can also be used to improve speed significantly if less accuracy is required. !or improving accuracy for extreme models. !Note this does not increase lSamp%l sampling or massive neutrino q-sampling - - real(sp) :: lAccuracyBoost=1. +!wig + real(sp) :: lAccuracyBoost=2._dl +!end wig !Boost number of multipoles integrated in Boltzman heirarchy integer, parameter :: lmin = 2 @@ -212,6 +216,13 @@ ! lmax is max possible number of l's evaluated integer, parameter :: lmax_arr = l0max + +!wig + real(dl) :: kSrcSampleBoost=1._dl + real(dl) :: kIntSampleBoost=1._dl + logical :: DoClsInterpol=.true. + logical :: ComputeEveryl=.true. +!end wig character(LEN=1024) :: highL_unlensed_cl_template = 'HighLExtrapTemplate_lenspotentialCls.dat' !fiducial high-accuracy high-L C_L used for making small cosmology-independent numerical corrections @@ -581,6 +592,25 @@ integer lind, lvar, step,top,bot,ls(lmax_arr) real(dl) AScale +!wig + if (ComputeEveryl) then + write(*,*) 'Warning: computing every l' + if (max_l.gt.lmax_arr) then + write(*,*)'lmax_arr maxl',lmax_arr,max_l + stop + endif + lind=0 + do lvar=2, max_l + lind=lind+1 + ls(lind)=lvar + end do + + lSet%l0=lind + lSet%l = ls(1:lind) + + else +!end wig + Ascale=scale/lSampleBoost if (lSampleBoost >=50) then @@ -720,6 +750,9 @@ end if lSet%l0=lind lSet%l(1:lind) = ls(1:lind) +!wig + endif +!end wig end subroutine initlval @@ -738,8 +771,10 @@ if (max_ind > lSet%l0) stop 'Wrong max_ind in InterpolateClArr' xl = real(lSet%l(1:lSet%l0),dl) - call spline(xl,iCL(1),max_ind,cllo,clhi,ddCl(1)) - + call spline(xl,iCl(1),max_ind,cllo,clhi,ddCl(1)) +!wig + if (DoClsInterpol) then + llo=1 do il=lmin,lSet%l(max_ind) xi=il @@ -755,7 +790,22 @@ +(b0**3-b0)*ddCl(lhi))*ho**2/6 end do - + + else + llo=1 + do il=2,lSet%l(max_ind) + xi=il + if ((xi > lSet%l(llo+1)).and.(llo < max_ind)) then + llo=llo+1 + end if + + all_Cl(il) = iCl(llo+1) + + end do + all_Cl(2)=iCl(1) + endif +!end wig + end subroutine InterpolateClArr subroutine InterpolateClArrTemplated(lSet,iCl, all_Cl, max_ind, template_index) @@ -2041,6 +2091,10 @@ dlnkh = 0.02 points = log(MTrans%TransferData(Transfer_kh,MTrans%num_q_trans,itf)/minkh)/dlnkh+1 ! dlnkh = log(MTrans%TransferData(Transfer_kh,MTrans%num_q_trans,itf)/minkh)/(points-0.999) +!wig + write(*,*)'Transfert_SaveToFiles: dlnkh= points= ',dlnkh,points + write(*,*)'You should check with transfer_k_per_logint > 0' +!endwig allocate(outpower(points,CP%InitPower%nn)) do in = 1, CP%InitPower%nn call Transfer_GetMatterPower(MTrans,outpower(1,in), itf, in, minkh,dlnkh, points) diff --git a/camb/params.ini b/camb/params.ini index cf5b821..54fdde6 100644 --- a/camb/params.ini +++ b/camb/params.ini @@ -1,12 +1,12 @@ #Parameters for CAMB #output_root is prefixed to output file names -output_root = test +output_root = tpl #What to do get_scalar_cls = T get_vector_cls = F -get_tensor_cls = F +get_tensor_cls = T get_transfer = F #if do_lensing then scalar_output_file contains additional columns of l^4 C_l^{pp} and l^3 C_l^{pT} @@ -65,18 +65,27 @@ nu_mass_degeneracies = 0 #Fraction of total omega_nu h^2 accounted for by each eigenstate, eg. 0.5 0.5 nu_mass_fractions = 1 -#Initial power spectrum, amplitude, spectral index and running. Pivot k in Mpc^{-1}. -initial_power_num = 1 -pivot_scalar = 0.05 -pivot_tensor = 0.05 -scalar_amp(1) = 2.1e-9 -scalar_spectral_index(1) = 0.96 -scalar_nrun(1) = 0 -tensor_spectral_index(1) = 0 -#ratio is that of the initial tens/scal power spectrum amplitudes -initial_ratio(1) = 1 -#note vector modes use the scalar settings above +#Type of initial power spectrum: power_sr or power_dbi +initial_power_num = 1 +power_name = power_sr + +#slow-roll parameters +ln[10^{10} P_scalar] = 0.303161E+01 +eps1(1) = 0.210391E-02 +eps2(1) = 0.421584E-01 +eps3(1) = 0 +#dbi parameter +gam(1) = 1. +del1(1) = 0. + +#wiggles parameters +logeps1invsigma0(1) = 1.5 +Xsigma0(1) = 0.1 +varphi(1) = 0 + +#tpl pivot +k_star = 0.05 #Reionization, ignored unless reionization = T, re_redshift measures where x_e=0.5 reionization = T @@ -98,6 +107,7 @@ RECFAST_fudge_He = 0.86 RECFAST_Heswitch = 6 RECFAST_Hswitch = T + #Initial scalar perturbation mode (adiabatic=1, CDM iso=2, Baryon iso=3, # neutrino density iso =4, neutrino velocity iso = 5) initial_condition = 1 @@ -121,7 +131,7 @@ CMB_outputscale = 7.4311e12 #transfer_k_per_logint=5 samples fixed spacing in log-k #transfer_interp_matterpower =T produces matter power in regular interpolated grid in log k; # use transfer_interp_matterpower =F to output calculated values (e.g. for later interpolation) -transfer_high_precision = F +transfer_high_precision = T transfer_kmax = 2 transfer_k_per_logint = 0 transfer_num_redshifts = 1 @@ -208,7 +218,7 @@ high_accuracy_default=F #Increase accuracy_boost to decrease time steps, use more k values, etc. #Decrease to speed up at cost of worse accuracy. Suggest 0.8 to 3. -accuracy_boost = 1 +accuracy_boost = 4 #Larger to keep more terms in the hierarchy evolution. l_accuracy_boost = 1 @@ -218,3 +228,15 @@ l_accuracy_boost = 1 #interpolation errors may be up to 3% #Decrease to speed up non-flat models a bit l_sample_boost = 1 + +# Boost the sampling on k sources only, in unit of AccuracyBoost +k_src_sample_boost = 1. + +# Boost the sampling on k line of sight integral only, in unit of AccuracyBoost +k_int_sample_boost = 1. + +#Use spline interpolation for the Cls +do_cls_interpol = F + +# Remove undersampling of the Cls: if true lSampleBoost has no effect +compute_every_l = T diff --git a/camb/power_sr.f90 b/camb/power_sr.f90 new file mode 100644 index 0000000..5a518d2 --- /dev/null +++ b/camb/power_sr.f90 @@ -0,0 +1,544 @@ +!This module provides the initial power spectra, parameterized as an expansion in ln k +! +! ln P = ln A_s + (n_s -1)*ln(k/k_0) + n_{run}/2 * ln(k/k_0)^2 +! +! so if n_run = 0 +! +! P = A_s (k/k_0_scalar)^(n_s-1) +! +!for the scalar spectrum, when an(in) is the in'th spectral index. k_0_scalar +!is a pivot scale, fixed here to 0.05/Mpc (change it below as desired). +! +!This module uses the same inputs an(in), ant(in) and rat(in) as CMBFAST, however here +!rat(in) is used to set the ratio of the initial power spectra, so here +! +!** rat(in) is not the Cl quadrupole ratio *** +! +!in general models the quadrupole ratio depends in a complicated way on the ratio of the initial +!power spectra + +!The absolute normalization of the Cls is unimportant here, but the relative ratio +!of the tensor and scalar Cls generated with this module will be correct for general models + + +!The OutputNormalization parameter controls the final output +!Absolute Cls can be obtained by setting OuputNormalization=outNone, otherwise the overall normalization +!of the power spectra doesn't matter + +!This version December 2003 - changed default tensor pivot to 0.05 (consistent with CMBFAST 4.5) + + module InitialPower + use Precision + implicit none + + private + +!wig character(LEN=*), parameter :: Power_Name = 'power_tilt' + character(80) :: Power_Name = 'power_sr' + logical, parameter :: useTplPhase = .false. + logical, save :: useSecondOrder = .false. +!end wig + integer, parameter :: nnmax= 5 + !Maximum possible number of different power spectra to use + + Type InitialPowerParams + + integer nn !Must have this variable + !The actual number of power spectra to use + + !For the default implementation return power spectra based on spectral indices + real(dl) an(nnmax) !scalar spectral indices + real(dl) n_run(nnmax) !running of spectral index + real(dl) ant(nnmax) !tensor spectral indices + real(dl) rat(nnmax) !ratio of scalar to tensor initial power spectrum amplitudes + real(dl) k_0_scalar, k_0_tensor + real(dl) ScalarPowerAmp(nnmax) + +!wig +! slow-roll power + real(dl) P_scalar + real(dl) k_star + real(dl) eps1(nnmax) + real(dl) eps2(nnmax) + real(dl) eps3(nnmax) +!dbi params + real(dl) gam(nnmax) + real(dl) del1(nnmax) +! tpl power + real(dl) invsigma0(nnmax) + real(dl) loginvsigma0(nnmax) + real(dl) logeps1invsigma0(nnmax) + real(dl) Xsigma0(nnmax) + real(dl) varphi(nnmax) +!Internal parameters + real(dl) ks_pivot !scalar pivot wave number + real(dl) kt_pivot !tensor pivot wave number + real(dl) ko_pivot !wiggles pivot wave number + real(dl) P_0_scalar(nnmax) !scalar normalisation during inflation + real(dl) P_0_tensor(nnmax) !tensor normalisation during inflation + real(dl) as0(nnmax) !scalar power spectrum coefficients: order 0 + real(dl) as1(nnmax) !order 1 + real(dl) as2(nnmax) !order 2 + real(dl) amps0(nnmax) !amplitudes of tpl effect in scalar power: order 0 + real(dl) amps1(nnmax) !order 1 + real(dl) amps2(nnmax) !order 2 + real(dl) args0(nnmax) !argument of the modulations: order 0 + real(dl) args1(nnmax) !order 1 + real(dl) args2(nnmax) !order 2 + real(dl) at0(nnmax) !Tensor power spectrum coefficients: same as scalars + real(dl) at1(nnmax) + real(dl) at2(nnmax) + real(dl) ampt0(nnmax) + real(dl) ampt1(nnmax) + real(dl) ampt2(nnmax) + real(dl) argt0(nnmax) + real(dl) argt1(nnmax) + real(dl) argt2(nnmax) +!end wig + + end Type InitialPowerParams + + real(dl) curv !Curvature contant, set in InitializePowers + + Type(InitialPowerParams) :: P + +!Make things visible as neccessary... + + public InitialPowerParams, InitialPower_ReadParams, InitializePowers, ScalarPower, TensorPower + public nnmax,Power_Descript, Power_Name, SetDefPowerParams +! public + contains + + + subroutine SetDefPowerParams(AP) + Type (InitialPowerParams) :: AP + + AP%nn = 1 !number of initial power spectra +!wig + + AP%an = 1 !scalar spectral index + AP%n_run = 0 !running of scalar spectral index + AP%ant = 0 !tensor spectra index + AP%rat = 1 + AP%k_0_scalar = 0.05 + AP%k_0_tensor = 0.05 + AP%ScalarPowerAmp = 1 + + AP%eps1 = 0.01 + AP%eps2 = 0.01 + AP%eps3 = 0. + AP%P_scalar=1.d-9 + AP%k_star = 0.05 + + AP%del1 = 0. + AP%gam = 1. + + + AP%Xsigma0 = 0. + AP%logeps1invsigma0 = 0. + AP%loginvsigma0 = AP%logeps1invsigma0 - log10(AP%eps1) + AP%invsigma0=10**(AP%loginvsigma0) + AP%varphi=0. + +!end wig + + end subroutine SetDefPowerParams + + subroutine InitializePowers(AParamSet,acurv) + Type (InitialPowerParams) :: AParamSet + !Called before computing final Cls in cmbmain.f90 + !Could read spectra from disk here, do other processing, etc. + + real(dl) acurv +!wig + integer :: i +!from Hankel functions + real(dl), parameter :: C_const = -0.7296 +!wkb or uniform approx gives a slightly different constant + real(dl), parameter :: D_const = -0.7652 +! real(dl), parameter :: C_const = D_const + + real(dl), parameter :: pi = 3.14159265 +!end wig + + if (AParamSet%nn > nnmax) then + write (*,*) 'To use ',AParamSet%nn,'power spectra you need to increase' + write (*,*) 'nnmax in power_tilt.f90, currently ',nnmax + end if + P = AParamSet + + curv=acurv + +!Write implementation specific code here... +!wig + + if ((Power_Name.eq.'power_sr').or.(Power_Name.eq.'power_dbi')) then + if(abs(curv) .gt. 1.d-3) then + write(*,*)'Non flat cosmology, curv = ',curv + write(*,*)'Cannot use SR or TPL spectrum.' + stop + endif + + P%ks_pivot = P%k_star + P%kt_pivot = P%k_star + P%ko_pivot = P%k_star + + + do i=1, P%nn + +!TODO: derive second order tpl corrections? + if (useSecondOrder) then + if (P%Xsigma0(i).ne.0.) then + useSecondOrder=.false. + write(*,*)'PowerIni: SR second order disabled: TPL not implemented! ' + endif + if ((P%del1(i).ne.0.).or.(P%gam(i).ne.1.)) then + useSecondorder = .false. + write(*,*)'PowerIni: SR second order disabled: DBI slow-roll not implemented!' + endif + endif + + if ((power_Name.eq.'power_sr').and. & + ((P%del1(i).ne.0.).or.(P%gam(i).ne.1.)) ) then + write(*,*)'power_sr with non-vanishing DBI params!' + write(*,*)'resetting to standard SR' + P%del1(i) = 0. + P%gam(i) = 1. + endif + + +!first order tpl corrections + + P%loginvsigma0(i) = P%logeps1invsigma0(i) - log10(P%eps1(i)) + P%invsigma0(i)=10**(P%loginvsigma0(i)) + + P%amps0(i) = 2.*P%Xsigma0(i) + P%amps1(i) = pi*P%Xsigma0(i)*(2.*P%eps1(i) + P%eps2(i)) + P%amps2(i) = 0. + if (useTplPhase) then + P%args0(i) = 2.*P%invsigma0(i)*(1. + P%eps1(i)) + P%varphi(i) + else + P%args0(i) = P%varphi(i) + endif + P%args1(i) = 2.*P%invsigma0(i)*P%eps1(i) + P%args2(i) = 0. + + P%ampt0(i) = P%amps0(i) + P%ampt1(i) = pi*P%Xsigma0(i)*2.*P%eps1(i) + P%ampt2(i) = 0. + P%argt0(i) = P%args0(i) + P%argt1(i) = P%args1(i) + P%argt2(i) = 0. + + +!first and second order slowroll + + P%as0(i) = 1. + P%eps1(i)*(-2.*(C_const+1.))-P%eps2(i)*C_const & + + (C_const + 2.)*P%del1(i) + P%as1(i) = -2.*P%eps1(i) - P%eps2(i) + P%del1(i) + P%as2(i) = 0. + + P%at0(i) = 1. + P%eps1(i)*(-2.*(C_const+1.)) + P%at1(i) = -2.*P%eps1(i) + P%at2(i) = 0. + + if (useSecondOrder) then + + + P%as0(i) = P%as0(i) & + + P%eps1(i)*P%eps1(i) * (2.*C_const**2 + 2.*C_const + pi**2/2. - 5.) & + + P%eps1(i)*P%eps2(i) * (C_const**2 - C_const + 7.*pi**2/12. - 7.) & + + P%eps2(i)*P%eps2(i) * (0.5*C_const**2 + pi**2/8. - 1.) & + + P%eps2(i)*P%eps3(i) * (-0.5*C_const**2 + pi**2/24.) + + P%as1(i) = P%as1(i) & + + P%eps1(i)*P%eps1(i) * 2.*(2.*C_const + 1.) & + + P%eps1(i)*P%eps2(i) * (2.*C_const - 1.) & + + P%eps2(i)*P%eps2(i) * C_const & + + P%eps2(i)*P%eps3(i) * (-C_const) + + P%as2(i) = P%as2(i) & + + P%eps1(i)*P%eps1(i) * 4. & + + P%eps1(i)*P%eps2(i) * 2. & + + P%eps2(i)*P%eps2(i) * 1. & + + P%eps2(i)*P%eps3(i) * (-1.) + + + P%at0(i) = P%at0(i) & + + P%eps1(i)*P%eps1(i) * (2.*C_const**2 + 2.*C_const + pi**2/2. - 5.) & + + P%eps1(i)*P%eps2(i) * (-C_const**2 - 2.*C_const + pi**2/12. - 2.) + + P%at1(i) = P%at1(i) & + + P%eps1(i)*P%eps1(i) * 2.*(2.*C_const + 1.) & + + P%eps1(i)*P%eps2(i) * 2.*(-C_const - 1.) + + P%at2(i) = P%at2(i) & + + P%eps1(i)*P%eps1(i) * 4. & + + P%eps1(i)*P%eps2(i) * (-2.) + + endif + + + !Relative normalisation of scalar and tensor spectra + + P%P_0_scalar(i)=P%P_scalar/(P%as0(i)) +! P%P_0_tensor(i)=P%P_scalar/(P%as0(i))*P%eps1(i)*16. + P%P_0_tensor(i)=P%P_scalar/(P%as0(i))*P%eps1(i)*16/P%gam(i) + + enddo + + else + stop 'No such initial power spectra' + endif +!end wig + + end subroutine InitializePowers + + + function ScalarPower(k,in) + + !"in" gives the index of the power to return for this k + !ScalarPower = const for scale invariant spectrum + !The normalization is defined so that for adiabatic perturbations the gradient of the 3-Ricci + !scalar on co-moving hypersurfaces receives power + ! < |D_a R^{(3)}|^2 > = int dk/k 16 k^6/S^6 (1-3K/k^2)^2 ScalarPower(k) + !In other words ScalarPower is the power spectrum of the conserved curvature perturbation given by + !-chi = \Phi + 2/3*\Omega^{-1} \frac{H^{-1}\Phi' - \Psi}{1+w} + !(w=p/\rho), so < |\chi(x)|^2 > = \int dk/k ScalarPower(k). + !Near the end of inflation chi is equal to 3/2 Psi. + !Here nu^2 = (k^2 + curv)/|curv| + + !This power spectrum is also used for isocurvature modes where + !< |\Delta(x)|^2 > = \int dk/k ScalarPower(k) + !For the isocurvture velocity mode ScalarPower is the power in the neutrino heat flux. + + real(dl) ScalarPower,k, lnrat + integer in + +!wig + real :: angle + +! lnrat = log(k/P%k_0_scalar) +! ScalarPower=P%ScalarPowerAmp(in)*exp((P%an(in)-1)*lnrat + P%n_run(in)/2*lnrat**2) + + lnrat = log(k/P%ks_pivot) + + + if ((Power_Name.eq.'power_sr').or.(Power_Name.eq.'power_dbi')) then + + angle = P%args0(in) + P%args1(in) * log(k/P%ko_pivot) + + ScalarPower = (P%as0(in) + P%as1(in)*lnrat + 0.5*P%as2(in)*lnrat**2) & + * (1. - P%amps0(in)*cos(angle)) - P%amps1(in)*sin(angle) + + else + stop 'No such initial power spectra' + endif + + ScalarPower = ScalarPower*P%P_0_scalar(in) + + if(ScalarPower .lt. 0.0) then + ScalarPower =0. + end if + + +!end wig + +! ScalarPower = ScalarPower * (1 + 0.1*cos( lnrat*30 ) ) + + end function ScalarPower + + + function TensorPower(k,in) + + !TensorPower= const for scale invariant spectrum + !The normalization is defined so that + ! < h_{ij}(x) h^{ij}(x) > = \sum_nu nu /(nu^2-1) (nu^2-4)/nu^2 TensorPower(k) + !for a closed model + ! < h_{ij}(x) h^{ij}(x) > = int d nu /(nu^2+1) (nu^2+4)/nu^2 TensorPower(k) + !for an open model + !"in" gives the index of the power spectrum to return + !Here nu^2 = (k^2 + 3*curv)/|curv| + + + real(dl) TensorPower,k + real(dl), parameter :: PiByTwo=3.14159265d0/2._dl + + integer in + +!wig + real :: lnrat,angle + +! TensorPower=P%rat(in)*P%ScalarPowerAmp(in)*exp(P%ant(in)*log(k/P%k_0_tensor)) +! if (curv < 0) TensorPower=TensorPower*tanh(PiByTwo*sqrt(-k**2/curv-3)) + + lnrat = log(k/P%kt_pivot) + + if ((Power_Name.eq.'power_sr').or.(Power_Name.eq.'power_dbi')) then + + angle = P%argt0(in) + P%argt1(in) * log(k/P%ko_pivot) + + TensorPower = (P%at0(in) + P%at1(in)*lnrat + 0.5*P%at2(in)*lnrat**2) & + * (1. - P%ampt0(in)*cos(angle)) - P%ampt1(in)*sin(angle) + + else + stop 'No such initial power spectra' + endif + + TensorPower = TensorPower*P%P_0_tensor(in) + + if (curv < 0) TensorPower=TensorPower*tanh(PiByTwo*sqrt(-k**2/curv-3)) + + if(TensorPower .lt. 0.0) then + TensorPower=0. + endif + +!end wig + + + end function TensorPower + + !Get parameters describing parameterisation (for FITS file) + function Power_Descript(in, Scal, Tens, Keys, Vals) + character(LEN=8), intent(out) :: Keys(*) + real(dl), intent(out) :: Vals(*) + integer, intent(IN) :: in + logical, intent(IN) :: Scal, Tens + integer num, Power_Descript + num=0 +!wig + + + if (Power_Name.eq.'power_sr') then + + num=num+1 + Keys(num) = 'k_star' + Vals(num) = P%k_star + num=num+1 + Keys(num) = 'P_scalar' + Vals(num) = P%P_scalar + num=num+1 + Keys(num) = 'eps1' + Vals(num) = P%eps1(in) + num=num+1 + Keys(num) = 'eps2' + Vals(num) = P%eps2(in) + num=num+1 + Keys(num) = 'eps3' + Vals(num) = P%eps3(in) + num=num+1 + Keys(num) = 'isigma0' + Vals(num) = P%invsigma0(in) + num=num+1 + Keys(num) = 'lisigma0' + Vals(num) = P%loginvsigma0(in) + num=num+1 + Keys(num) = 'Xsigma0' + Vals(num) = P%Xsigma0(in) + num=num+1 + Keys(num) = 'varphi' + Vals(num) = P%varphi(in) + + elseif (Power_Name.eq.'power_dbi') then + + num=num+1 + Keys(num) = 'k_star' + Vals(num) = P%k_star + num=num+1 + Keys(num) = 'P_scalar' + Vals(num) = P%P_scalar + num=num+1 + Keys(num) = 'eps1' + Vals(num) = P%eps1(in) + num=num+1 + Keys(num) = 'eps2' + Vals(num) = P%eps2(in) + num=num+1 + Keys(num) = 'eps3' + Vals(num) = P%eps3(in) + num=num+1 + Keys(num) = 'gam' + Vals(num) = P%gam(in) + num=num+1 + Keys(num) = 'del1' + Vals(num) = P%del1(in) + num=num+1 + Keys(num) = 'isigma0' + Vals(num) = P%invsigma0(in) + num=num+1 + Keys(num) = 'lisigma0' + Vals(num) = P%loginvsigma0(in) + num=num+1 + Keys(num) = 'Xsigma0' + Vals(num) = P%Xsigma0(in) + num=num+1 + Keys(num) = 'varphi' + Vals(num) = P%varphi(in) + + else + stop 'No such initial power spectra' + endif +!end wig + Power_Descript = num + + end function Power_Descript + + subroutine InitialPower_ReadParams(InitPower, Ini, WantTensors) + use IniFile + Type(InitialPowerParams) :: InitPower + Type(TIniFile) :: Ini + logical, intent(in) :: WantTensors + integer i + + InitPower%k_0_scalar = Ini_Read_Double_File(Ini,'pivot_scalar',InitPower%k_0_scalar) + InitPower%k_0_tensor = Ini_Read_Double_File(Ini,'pivot_tensor',InitPower%k_0_tensor) + InitPower%nn = Ini_Read_Int('initial_power_num') + if (InitPower%nn>nnmax) stop 'Too many initial power spectra - increase nnmax in InitialPower' + +!wig + + Power_Name = Ini_Read_String('power_name') + + if ((Power_Name.eq.'power_sr').or.(Power_Name.eq.'power_dbi')) then + + do i=1, InitPower%nn + +! if (P%WantScalars .or. P%WantTransfer) then + + InitPower%eps1(i)= Ini_Read_Double_Array_File(Ini,'eps1',i,0.01d0) + InitPower%eps2(i)=Ini_Read_Double_Array_File(Ini,'eps2',i,0.01d0) + InitPower%eps3(i)=Ini_Read_Double_Array_File(Ini,'eps3',i,0.00d0) + InitPower%gam(i)=Ini_Read_Double_Array_File(Ini,'gam',i,1.00d0) + InitPower%del1(i)=Ini_Read_Double_Array_File(Ini,'del1',i,0.00d0) + InitPower%logeps1invsigma0(i) & + = Ini_Read_Double_Array_File(Ini,'logeps1invsigma0',i,1d0) + InitPower%loginvsigma0(i) = InitPower%logeps1invsigma0(i) & + - log10(InitPower%eps1(i)) + InitPower%invsigma0(i) = 10**(InitPower%loginvsigma0(i)) + InitPower%Xsigma0(i) = Ini_Read_Double_Array_File(Ini,'Xsigma0',i,0d0) + InitPower%varphi(i) = Ini_Read_Double_Array_file(Ini,'varphi',i,0.d0) + +! endif + + enddo + + InitPower%k_star= Ini_Read_Double('k_star',0.05d0) + InitPower%P_scalar= 1e-10*exp(Ini_Read_Double('ln[10^{10} P_scalar]',1d0)) + + else + + stop 'No such initial power spectrum' + + endif + +!end wig + + + + + + + end subroutine InitialPower_ReadParams + + + end module InitialPower diff --git a/camb/readme.html b/camb/readme.html index f5ada26..cff517b 100644 --- a/camb/readme.html +++ b/camb/readme.html @@ -112,7 +112,7 @@ where all quantities are in the synchronous gauge and evaluated at the requested

Version history

January 2012
December 2011
@@ -503,7 +503,7 @@ the Cls are normalized. Comments in the code explain this further. This file defines a module called Reionization that parameterizes the reionization history and supplies a function Reionization_xe that gives xe as a function of redshift. Optical depth input parameters are mapped into zre (defined as where xe is half its maximum (ex second He reionization)) using a binary search. See the CAMB notes for discussion. This module should be easily modifiable for alternative reionization models.

halofit.f90

-Implements the NonLinear module, to calculate non linear scalings of the matter power spectrum as a function of redshift. Uses HALOFIT (astro-ph/0207664, code thanks to Robert Smith. Note this is only reliable at the several percent level for standard ΛCDM models with power law initial power spectra. This module can be replaced to use a different non-linear fitting method. +Implements the NonLinear module, to calculate non linear scalings of the matter power spectrum as a function of redshift. Uses HALOFIT (astro-ph/0207664, code thanks to Robert Smith, with tweaks from arXiv:1109.4416. Note this is only reliable at the several percent level for standard ΛCDM models with power law initial power spectra. This module can be replaced to use a different non-linear fitting method.

@@ -563,7 +563,9 @@ Implements the NonLinear module, to calculate non linear scalings of the matter

Accuracy

-Scalar errors should rarely exceed 0.3% for min(2500, L well into the damping tail) at default accuracy setting, and 0.1% for 500<L<2000 with high_accuracy_default=T. Matter power spectrum errors are usually dominated by interpolation in the acoustic oscillations, with about 0.2% accuracy with high_accuracy_default (but much better rms accuracy). +Scalar numerical errors should rarely exceed 0.3% for min(2500, L well into the damping tail) at default accuracy setting, and 0.1% for 500<L<2000 with high_accuracy_default=T. Matter power spectrum errors are usually dominated by interpolation in the acoustic oscillations, with about 0.2% accuracy with high_accuracy_default (but much better rms accuracy). For a detailed study of numerical accuracy as of January 2012 see arXiv:1201.3654. + + See also comparison with CMBFAST. Accuracy of course assumes the model is correct, and is dependent on RECFAST being the correct ionization history. Lensed C_l TT, TE and EE are accurate at the same level (to within the approximation that the lensing potential is linear, or the accuracy of the the HALOFIT non-linear model).

Extreme models (e.g. scale > 4, h>1) may give errors of 5% or more. @@ -594,10 +596,15 @@ See also comparison with CMBFAST. Accuracy of course assu

REFERENCES

-Some notes and relevant Maple derivations are given here (see also the Appendix of astro-ph/0406096). The CAMB notes outline the equations and approximations used, and relation to standard synchronous-gauge and Newtonian-gauge variables. +Some notes and relevant Maple derivations are given here (see also the Appendix of astro-ph/0406096). The CAMB notes outline the equations and approximations used, and relation to standard synchronous-gauge and Newtonian-gauge variables; see also arXiv:1201.3654. There is a BibTex file of references (including CosmoMC).

+CMB power spectrum parameter degeneracies in the era of precision cosmology
+Cullan Howlett, Antony Lewis, Alex Hall, Anthony Challinor arXiv:1201.3654. +

Efficient computation of CMB anisotropies in closed FRW Models
Antony Lewis, Anthony Challinor and Anthony Lasenby astro-ph/9911177 Ap. J. 538:473-476, 2000. @@ -647,7 +654,8 @@ Antony Lewis, astro-ph/0403583HALOFIT

Stable clustering, the halo model and nonlinear cosmological power spectra
-Smith, R. E. and others, astro-ph/0207664 +Smith, R. E. and others, astro-ph/0207664. +

RECOMBINATION

@@ -696,13 +704,25 @@ The Cosmic Linear Anisotropy Solving System (CLASS) II: Blas, Diego and Lesgourgues, Julien and Tram, Thomas. arXiv:1104.2933

+Massive Neutrinos +

+CMB power spectrum parameter degeneracies in the era of precision cosmology +
+Cullan Howlett, Antony Lewis, Alex Hall, Anthony Challinor. +arXiv:1201.3654 +

+

Evolution of cosmological dark matter perturbations
+Antony Lewis and Anthony Challinor astro-ph/0203507 +Phys. Rev. D66, 023531 (2002) + +

Synchronous gauge theory and non-flat models

Complete treatment of CMB anisotropies in a FRW universe
Wayne Hu, Uros Seljak and Matias Zaldarriaga. Phys. Rev. D57:6, 3290-3301, 1998. astro-ph/9709066. -

WKB approx to hyperspherical Bessel functions @@ -719,10 +739,16 @@ Blas, Diego and Lesgourgues, Julien and Tram, Thomas. astro-ph/9603033 Ap.J. 469:2 437-444, 1996

+ +Integral solution for the microwave background + anisotropies in nonflat universes
+ Matias Zaldarriaga, Uros Seljak, Edmund Bertschinger. + ApJ. 494:491-501, 1998. astro-ph/9704265. + +

CMBFAST for spatially closed universes
Uros Seljak and Matias Zaldariaga, astro-ph/9911219 -

-See also the references on the CMBFAST home page. + diff --git a/distparams.ini b/distparams.ini index 35474e0..3f70f4d 100644 --- a/distparams.ini +++ b/distparams.ini @@ -54,6 +54,7 @@ line_labels= F plot_meanlikes = T shade_meanlikes = T +td_meanlikes = F # if non-zero, output _thin file, thinned by thin_factor thin_factor = 0 @@ -99,7 +100,7 @@ matlab_colscheme = #Parameters to use. If zero use all parameters which have lables. plotparams_num = 0 -plotparams = omegabh2 omegadmh2 tau ns +plotparams = omegabh2 omegadmh2 tau ns eps1 eps2 eps3 logomega amp psi logsig #Get set label to empty to not include a parameter in the parameter_names file #lab[asz]= diff --git a/params.ini b/params.ini index 6e524dd..5c3ff64 100644 --- a/params.ini +++ b/params.ini @@ -72,18 +72,18 @@ directional_grid_steps = 20 #use fast-slow parameter distinctions to speed up #(note for basic models WMAP3 code is only ~3x as fast as CAMB) -use_fast_slow = F +use_fast_slow = T #Can use covariance matrix for proposal density, otherwise use settings below #Covariance matrix can be produced using "getdist" program. -propose_matrix = params_CMB.covmat +#propose_matrix = params_CMB.covmat #If propose_matrix is blank (first run), can try to use numerical Hessian to #estimate a good propose matrix. As a byproduct you also get an approx best fit point estimate_propose_matrix = F #Tolerance on log likelihood to use when estimating best fit point -delta_loglike = 2 +delta_loglike = 0.1 #Scale of proposal relative to covariance; 2.4 is recommended by astro-ph/0405462 for Gaussians #If propose_matrix is much broader than the new distribution, make proportionately smaller @@ -150,7 +150,7 @@ highL_unlensed_cl_template = ./camb/HighLExtrapTemplate_lenspotentialCls.dat #CAMB parameters #If we are including tensors -compute_tensors = F +compute_tensors = T #Initial power spectrum amplitude point (Mpc^{-1}) #Note if you change this, may need new .covmat as degeneracy directions change pivot_k = 0.05 @@ -166,7 +166,12 @@ high_accuracy_default = F #increase accuracy_level to run CAMB on higher accuracy #(default is about 0.3%, 0.1% if high_accuracy_default =T) #accuracy_level can be increased to check stability/higher accuracy, but inefficient -accuracy_level = 1 +accuracy_level = 4 +l_accuracy_level = 2 +l_sample_level = 1 +k_src_sample_level = 1 +k_int_sample_level = 1 +camb_compute_every_l = T #If action = 1 redo_likelihoods = T @@ -192,12 +197,19 @@ param[omegak] = 0 0 0 0 0 param[fnu] = 0 0 0 0 0 param[w] = -1 -1 -1 0 0 -param[ns] = 0.95 0.5 1.5 0.02 0.01 -param[nt] = 0 0 0 0 0 -param[nrun] = 0 0 0 0 0 +param[sr1] = 0.1 0.00001 0.2 0.01 0.02 +param[sr2] = -0.04 -0.2 0.2 0.01 0.02 +param[sr3] = -0.04 -0.2 0.2 0.01 0.02 +param[unused] = 0 0 0 0 0 -#log[10^10 A_s] -param[logA] = 3 2.7 4 0.01 0.01 -param[r] = 0 0 0 0 0 +param[logfreq] = 1.5 1 2.5 0.5 0.2 +param[amposc] = 0.05 0 0.45 0.1 0.2 +param[phase] = 0 0 6.28 0.4 0.1 + +param[logsnd] = 0 0 0 0 0 + +#ln[10^10 A_s] +param[lnA] = 3 2.7 4 0.01 0.01 +param[r] = 1 1 1 0 0 #SZ amplitude, as in WMAP analysis param[asz]= 1 0 2 0.4 0.4 diff --git a/params_CMB.paramnames b/params_CMB.paramnames index 0c08e67..8fa3488 100644 --- a/params_CMB.paramnames +++ b/params_CMB.paramnames @@ -1,20 +1,25 @@ -omegabh2 \Omega_b h^2 #physical baryon density -omegadmh2 \Omega_{DM} h^2 #physical dark matter density, including CDM and massive neutrinos -theta \theta #100 times the ratio of the angular diameter distance to the LSS sound horizon (approx) +omegabh2 \Omega_b h^2 #physical baryon density +omegadmh2 \Omega_{DM} h^2 #physical dark matter density, including CDM and massive neutrinos +theta \theta #100 times the ratio of the angular diameter distance to the LSS sound horizon (approx) tau \tau omegak \Omega_K -fnu f_\nu #neutrino energy density as fraction of omegadmh2 -w w #constant equation of state parameter for scalar field dark energy -ns n_s #beware that pivot scale can change in .ini file -nt n_t -nrun n_{run} -logA log[10^{10} A_s] -r r #ratio of tensor to scalar primordial amplitudes at pivot scale -asz A_{SZ} #SZ template amplitude, as in WMAP +fnu f_\nu #neutrino energy density as fraction of omegadmh2 +w w #constant equation of state parameter for scalar field dark energy +sr1 \epsilon_1-\delta_1 #first slow roll parameter - first sound flow param +sr2 \epsilon_2+\delta_1 #second slow roll parameter + first sound flow param +sr3 \epsilon_3 #third slow roll parameter +unused log(\epsilon_1) #unused +logfreq log(\epsilon_1/\sigma_0) #superimposed oscillation frequency +amposc x\sigma_0 #superimposed oscillation amplitude +phase \psi #superimposed oscillation phase shift +logsnd log(\epsilon_1 c_s) #speed of sound +lnA ln[10^{10} P_*] +r r #ratio of tensor to scalar primordial amplitudes at pivot scale +asz A_{SZ} #SZ template amplitude, as in WMAP omegal* \Omega_\Lambda age* Age/Gyr omegam* \Omega_m sigma8* \sigma_8 zrei* z_{re} -r10* r_{10} #tensor-scalar C_l amplitude at l=10 -H0* H_0 #hubble parameter is H0 km/s/Mpc \ No newline at end of file +r10* r_{10} #tensor-scalar C_l amplitude at l=10 +H0* H_0 #hubble parameter is H0 km/s/Mpc diff --git a/readme.html b/readme.html index 9c8fd78..69a24f9 100644 --- a/readme.html +++ b/readme.html @@ -171,8 +171,12 @@ If you don't have a propose_matrix, set estimate_propose_matrix = T to au

  • Sampling method
    Set sampling_method=1 to use the default Metropolis algorithm (in the optimal fast/slow subspaces if use_fast_slow is true). Other options are slice sampling (2), slice sampling fast parameters (3), and directional gridding (4). These methods should work fine for most simple distributions; the temperature input parameter can be increased to probe further into the tails if required (e.g. to get better high-confidence limits). Further sampling_method options that you can try for nastier (e.g. multi-modal distributions) are multicanonical sampling (5) and Wang-Landau-like sampling (6). These latter methods could be modified to calculate the evidence, but at the moment are only implemented to sample nastier distributions via importance sampling. Multicanonical sampling probes into the tails a distance proportional to the running time, and all samples can be kept if the distribution turns out to be unimodal. The Wang-Landau-like sampling probes the full likelihood range from the word go, and produces no samples for the first 10,000 or so steps (thereafter all samples are strictly Markovian may be kept without burn in). Methods 5 and 6 can be used with MPI, but MPI stopping should probably be turned off; MPI proposal learning may work with method 6, though this is not extensively tested. Both are likely to require samples = 100000 or larger, hence significantly slow than a basic MCMC run for simple distributions. Settings for methods 5 and 6 can be edited at the top of MCMC.f90. See the notes for a more detailed explanation of sampling methods. +

    +
  • CAMB's numerical accuracy
    +The default accuracy is aimed at WMAP. Set high_accuracy_default = T to target ~0.1% accuracy on the CMB power spectra from CAMB. This is sufficient for numerical biases to be small compared to the error bar for Planck; for an accuracy study see arXiv:1201.3654. You can also use the accuracy_level parameter to increase accuracy further, however this is numerically inefficient and is best used mainly for numerical stability checks using importance sampling (action=1).

    +
  • Best-fit point
    Set action =2 to just calculate the best-fit point and stop. Set delta_loglike to the tolerance on the log likelihood for finding the best fit. The values are output to a file called file_root.minimum. Note that this function does not always work very reliably. @@ -185,7 +189,6 @@ to zero and remove entries when you come to process the output samples.

    -
  • Threads and run-time adjustment
    The num_threads parameter will determine the number of openMP threads (in MPI runs, usually set to the number of CPUs on each node). Scaling is linear up to about 8 processors on most systems, then falls off slowly. It is probably best to run several chains on 8 or fewer @@ -523,8 +526,8 @@ weighting function or parameter mapping.
  • January 2012
    diff --git a/source/CMB_Cls_simple.f90 b/source/CMB_Cls_simple.f90 index db9d905..9b9617d 100644 --- a/source/CMB_Cls_simple.f90 +++ b/source/CMB_Cls_simple.f90 @@ -404,6 +404,28 @@ contains integer zix real redshifts(matter_power_lnzsteps) +!wig + AccuracyBoost = AccuracyLevel + lAccuracyBoost = lAccuracyLevel + lSampleBoost = lSampleLevel + kSrcSampleBoost = kSrcSampleLevel + kIntSampleBoost = kIntSampleLevel + ComputeEveryl = CambComputeEveryl + + write(*,*) + write(*,*) '>>>>>>>>>>>>>> accuracy settings <<<<<<<<<<<<<<' + write(*,*) 'Power_Name = ',Power_Name + write(*,*) 'write_all_cls = ',write_all_cls + write(*,*) 'AccuracyBoost = ',AccuracyBoost + write(*,*) 'lAccuracyBoost = ',lAccuracyBoost + write(*,*) 'DoClsInterpol = ',DoClsInterpol + write(*,*) 'ComputeEveryl = ',ComputeEveryl + write(*,*) '>>>>>>>>>>>>>> accuracy settings <<<<<<<<<<<<<<' + write(*,*) + + if ((Power_Name == 'power_sr').and.(.not.compute_tensors)) stop 'Slow-Roll requires tensors' +!end wig + Threadnum =num_threads w_lam = -1 call CAMB_SetDefParams(P) @@ -436,9 +458,11 @@ contains P%Transfer%high_precision=.true. P%Transfer%kmax=P%Transfer%kmax + 0.2 end if - AccuracyBoost = AccuracyLevel - lAccuracyBoost = AccuracyLevel - lSampleBoost = AccuracyLevel +!wig +! AccuracyBoost = AccuracyLevel +! lAccuracyBoost = AccuracyLevel +! lSampleBoost = AccuracyLevel +!end wig P%AccurateReionization = .true. end if diff --git a/source/GetDist.f90 b/source/GetDist.f90 index 65cbc71..e72ad74 100644 --- a/source/GetDist.f90 +++ b/source/GetDist.f90 @@ -105,6 +105,32 @@ module MCSamples integer Num_ComparePlots logical :: prob_label = .false. +!addon +!exact field model + logical, parameter :: findMaxValues = .false. + logical, parameter :: mapRhoReh = .false. + logical, parameter :: mapLFPowerToWreh = .true. + +!slow-roll analysis + logical, parameter :: adjEps1ToLog = .false. + logical, parameter :: extractHinf = .false. + logical, parameter :: mapnslogr = .true. + logical, parameter :: use2ndSR=.true. + real, parameter :: ln10p5HinfCut=-10 + real, parameter :: logrCut=-3.7 + +#ifdef SRREHEAT + logical, parameter :: doReheatFromSr = .false. + character(len=*), parameter :: reheatModel='smallf' +#endif + + real, parameter :: wfix = -0.3 + real, parameter :: pmax = 15, pmin = 0 + real, parameter :: logmumax=10 +!end addon + + + contains subroutine AdjustPriors @@ -112,8 +138,20 @@ contains !Be careful as this code is parameterisation dependent integer i real ombh2, chisq +!addon + integer, parameter :: dp = kind(1._8) + real(dp) :: eps1,logeps1,eps2,eps3, lnR,lnRrad + real(dp) :: wreh,p,logmu,lnRhoReh,lnRhoEnd + + real :: ln10p5HinfMin=1e30 + real :: ln10p5HinfMax=-1e30 + + real :: lnRhoRehmax = -1e30, lnRhoEndMax=-1e30, lnRmax = -1e30, lnRradmax = -1e30 + + real, parameter :: lnRhoNuc = -187.747 - stop 'You need to write the AdjustPriors subroutine in GetDist.f90 first!' +! stop 'You need to write the AdjustPriors subroutine in GetDist.f90 first!' +!end addon ! write (*,*) 'Adjusting priors' ! do i=0, nrows-1 @@ -125,6 +163,214 @@ contains ! ! end do +!addon + wreh = wfix + + + if (findMaxValues) then + print *,'-->computing maximum values in chains' + + do i=0,nrows-1 + lnR = coldata(14+2,i) + lnRhoEnd = coldata(2+30,i) + lnRrad = coldata(2+28,i) + lnRhoEndMax = max(lnRHoEndMax,lnRhoEnd) + lnRmax = max(lnRmax, lnR) + lnRradmax = max(lnRradmax, lnRrad) + enddo + + print *,'lnRhoEndMax= ',lnRhoEndMax + print *,'lnRradMax= lnRmax= ',lnRradMax, lnRmax + endif + + if (mapRhoReh) then + print *,'--> rejecting all p <',pmin + print *,'--> rejecting all logmu >',logmumax + print *,'--> rejecting all lnRhoNuc > lnRhoReh > lnRhoEnd' + + do i=0,nrows-1 + p = coldata(2+10,i) + logmu = coldata(11+2,i) + + if (mapLFPowerToWreh) wreh = (p-2._dp)/(p+2._dp) + + lnR = coldata(14+2,i) + lnRrad = coldata(2+28,i) + lnRhoReh = coldata(27+2,i) + lnRhoEnd = coldata(2+30,i) + + if ((wreh.lt.-1._dp/3._dp) & + .or.(p.lt.pmin).or.(logmu.gt.logmumax) & + .or.(lnRhoReh.gt.lnRhoEnd) & + .or.(lnRhoReh.lt.lnRhoNuc)) then + chisq = 1e30 + coldata(1,i) = 0. + coldata(2,i) = coldata(2,i) + chisq/2 + else + lnRhoRehMax = max(lnRHoRehMax,lnRhoReh) + lnRhoEndMax = max(lnRHoEndMax,lnRhoEnd) + lnRmax = max(lnRmax, lnR) + lnRradmax = max(lnRradmax, lnRrad) + endif + + enddo + print *,'Max values after applying prior!' + print *,'lnRhoEndMax= ',lnRhoEndMax + print *,'lnRhoRehMax= ',lnRhoRehMax + print *,'lnRradMax= lnRmax= ',lnRradMax, lnRmax + endif + + if (adjEps1ToLog) then + print *,'----> assumed mapped log(eps1) in col 8' + print *,'----> from flat eps1 to flat log(eps1)' + read(*,*) + do i=0, nrows-1 + + logeps1 = coldata(8+2,i) + eps1 = 10.**logeps1 + + coldata(1,i) = coldata(1,i)/(eps1*log(10.)) + coldata(2,i) = coldata(2,i) + + end do + endif + + + + if (extractHinf) then + + write (*,*)'Computing and adjusting bounds for ln(10^5 Hinf/mpl) in 17!' + + do i=0,nrows-1 + ln10p5HinfMax = max(ln10p5HinfMax,coldata(2+17,i)) + ln10p5HinfMin = min(ln10p5HinfMin,coldata(2+17,i)) + + if (coldata(2+17,i).lt.ln10p5HinfCut) then + chisq = 1e30 + coldata(1,i) = 0. + coldata(2,i) = coldata(2,i) + chisq/2 + endif + + enddo + + write(*,*)'----> ln(10^5 Hinf) Min= Max= ',ln10p5HinfMin,ln10p5HinfMax + + endif + + if (mapnslogr) then + + write (*,*)'Adjusting bounds for log(r) in 14!' + + do i=0,nrows-1 + + if (coldata(2+14,i).lt.logrCut) then + chisq = 1e30 + coldata(1,i) = 0. + coldata(2,i) = coldata(2,i) + chisq/2 + endif + + enddo + endif + +#ifdef SRREHEAT + if (doReheatFromSR) then + write(*,*)'Adjusting prior for reheat SR analysis' + + select case (reheatModel) + + case ('largef') + + + do i=0,nrows-1 + + eps1 = 10**(coldata(2+8,i)) + eps2 = coldata(2+9,i) + p = coldata(2+14,i) + wreh = (p-2.)/(p+2.) + +!large field models cannot produce eps2<=0 (out of prior) + if ((eps2.le.0.).or.(p.gt.pmax).or.(p.lt.pmin)) then + chisq = 1e30 + coldata(1,i) = 0. + coldata(2,i) = coldata(2,i) + chisq/2 + endif + +!reheating is after end of inflation + lnRhoReh = coldata(11+2,i) + lnRhoEnd = coldata(12+2,i) + +! print *,'lnrhoreh end=',lnRhoReh,lnRhoEnd,eps1,eps2 + + if (wreh.lt.-1./3.) then + chisq = 1e30 + coldata(1,i) = 0. + coldata(2,i) = coldata(2,i) + chisq/2 + endif + if (lnRhoReh.gt.lnRhoEnd) then + chisq = 1e30 + coldata(1,i) = 0. + coldata(2,i) = coldata(2,i) + chisq/2 + endif + if (lnRhoReh.lt.lnRhoNuc) then + chisq = 1e30 + coldata(1,i) = 0. + coldata(2,i) = coldata(2,i) + chisq/2 + endif + enddo + + case ('smallf') + + wreh = wfix + + do i=0,nrows-1 + + eps1 = 10**(coldata(2+8,i)) + eps2 = coldata(2+9,i) + eps3 = coldata(2+10,i) + + p = coldata(2+14,i) + + +!small field models (out of prior) + if (.not.(sr_checkobs_sf(eps1,eps2,eps3)) & + .or.(p.gt.pmax).or.(p.lt.pmin)) then + chisq = 1e30 + coldata(1,i) = 0. + coldata(2,i) = coldata(2,i) + chisq/2 + endif + + +!reheating is after end of inflation + lnRhoReh = coldata(11+2,i) + lnRhoEnd = coldata(12+2,i) + + if (wreh.lt.-1./3.) then + chisq = 1e30 + coldata(1,i) = 0. + coldata(2,i) = coldata(2,i) + chisq/2 + endif + if (lnRhoReh.gt.lnRhoEnd) then + chisq = 1e30 + coldata(1,i) = 0. + coldata(2,i) = coldata(2,i) + chisq/2 + endif + if (lnRhoReh.lt.lnRhoNuc) then + chisq = 1e30 + coldata(1,i) = 0. + coldata(2,i) = coldata(2,i) + chisq/2 + endif + enddo + + end select + + + endif +#endif +!end addon + + + + end subroutine AdjustPriors subroutine MapParameters(invars) @@ -132,7 +378,194 @@ contains ! map parameters in invars: eg. invars(3)=invars(3)*invars(4) ! invars(2+13)=invars(17+2)*exp(-invars(2+4)) - stop 'Need to write MapParameters routine first' + + +!addon + integer, save :: counter = 0 + real :: lnR, lnRhoReh,lnRhoEnd, ns,logr + real :: cseps1,eps1mdel1,eps2pdel1,del1 + real :: as0,ln10p5Hinf + real, parameter :: C_const = -0.7296 + + integer, parameter :: dp=kind(1._8) + real(dp) :: p, mu, wreh, Pstar,bfold,eps1,eps2,eps3 + real(dp), dimension(2) :: buffer2D + real(dp), dimension(3) :: buffer3D + + + counter = counter + 1 + + if (mapRhoReh) then + + if (mapLFPowerToWreh) then + + p = invars(2+10) + wreh = (p-2._dp)/(p+2._dp) + + if (counter.eq.1) print *,'reheat marginalisation on p in col10' + else + + wreh = wfix + if (counter.eq.1) print *,'marginalising at fixed wreh= ',wreh + + endif + + if (wreh.eq.1./3.) stop 'Treh inversion is singular!' + + if (counter.eq.1) print *,'mapping lnR -> lnRhoReh in col27' + + lnR = invars(2+14) + + lnRhoEnd = invars(2+30) + + lnRhoReh = 4._dp*((1._dp+wreh)/(1._dp/3._dp-wreh))*lnR & + - 4._dp*(1._dp/3._dp + wreh)/(1._dp/3._dp - wreh)*lnRhoEnd/2._dp + + invars(27+2)=lnRhoReh + + end if + + + if (extractHinf) then + del1 = 0. + + if (counter.eq.1) print *,'assumiing 10**(col10) is eps1' + + eps1mdel1 = 10**invars(2+8) + eps2pdel1 = invars(2+9) + eps1 = eps1mdel1 + eps2 = eps2pdel1 + eps3 = invars(2+10) + + + cseps1=eps1mdel1 +! cseps1 = eps1/invars(2+15) +! cseps1= 10.**invars(2+15) + +! eps1 = 10.**invars(2+8) +! eps2 = invars(2+9) +! eps3 = invars(2+10) + + + + as0 = 1. + eps1mdel1*(-2.*(C_const+1.)) - eps2pdel1*C_const + + if (use2ndSR) then + if (counter.eq.1) print *,'using 2nd order slow-roll!' + as0 = as0 + eps1*eps1 * (2.*C_const**2 + 2.*C_const + pi**2/2. - 5.) & + + eps1*eps2 * (C_const**2 - C_const + 7.*pi**2/12. - 7.) & + + eps2*eps2 * (0.5*C_const**2 + pi**2/8. - 1.) & + + eps2*eps3 * (-0.5*C_const**2 + pi**2/24.) + endif + + ln10p5Hinf = + 0.5*invars(2+16) & + + 0.5*log(3.141592653589793238*cseps1) - 0.5*log(as0) + + + invars(2+17) = ln10p5Hinf + endif + + if (mapnslogr) then + + if (counter.eq.1) then + print *,'assumiing 10**(col10) is eps1' + print *,'Mapping ns --> col15 & logr --> col14' + endif + + eps1 = 10**invars(2+8) + eps2 = invars(2+9) + eps3 = invars(2+10) + + ns = 1 - 2.*eps1 - eps2 + logr = log10(16*eps1) + + if (use2ndSR) then + if (counter.eq.1) print *,'using 2nd order slow-roll!' + ns = ns - 2*eps1**2 - (2*C_const + 3)*eps1*eps2 - C_const*eps2*eps3 + logr = logr + log10(1 + C_const*eps2) + endif + + invars(2+14) = logr + invars(2+15) = ns + + endif + + +#ifdef SRREHEAT + if (doReheatFromSR) then + + if (counter.eq.1) then + write(*,*)'************************************************' + write(*,*)'Model is: ',reheatModel + write(*,*)'Doing Reheating analysis from SR:' + write(*,*)'Mapping lnRhoReh -> col11 & lnRhoEnd -> col12' + write(*,*)'& -N* -> col13 & p -> col14 & mu -> col15' + write(*,*)'************************************************' + endif + + Pstar = exp(invars(2+16))*1.e-10 + + eps1 = 10**(invars(2+8)) + eps2 = invars(2+9) + eps3 = invars(2+10) + + lnRhoReh = 0. + lnRhoEnd = 0. + bfold=0. + + select case(reheatModel) + + case ('largef') + +!large field models cannot produce eps2<=0 (out of prior) + if (eps2.le.0.) return + + buffer2d = sr_matter_obs_lf(eps1,eps2) + + p = buffer2d(1) + invars(15+2)=p + + if ((p.gt.pmax).or.(p.lt.pmin)) return + + lnRhoReh = sr_lnrhoreh_lf(wreh,eps1,eps2,Pstar,bfold) + lnRhoEnd = sr_lnrhoend_lf(p,Pstar) + + invars(11+2)=lnRhoReh + invars(12+2)=lnRhoEnd + invars(13+2)=bfold + + case ('smallf') + +!small field models cannot produce any eps1,eps2,eps3<=0 (out of prior) + if (.not.sr_checkobs_sf(eps1,eps2,eps3)) return + + buffer3d = sr_matterovermu_obs_sf(eps1,eps2,eps3) + + p = buffer3d(1) + mu = buffer3d(2) + wreh = wfix + + invars(14+2)=p + invars(15+2)=mu + + if ((p.gt.pmax).or.(p.lt.pmin)) return + + lnRhoReh = sr_lnrhoreh_sf(wreh,eps1,eps2,eps3,Pstar,bfold) + lnRhoEnd = sr_lnrhoend_sf(p,mu,Pstar) + + invars(11+2) = lnRhoReh + invars(12+2) = lnRhoEnd + invars(13+2) = bfold + + case default + + stop 'model not implemented!' + + end select + + endif +#endif +!end addon end subroutine MapParameters diff --git a/source/MCMC.f90 b/source/MCMC.f90 index 9c3285d..c731cb4 100644 --- a/source/MCMC.f90 +++ b/source/MCMC.f90 @@ -215,7 +215,7 @@ contains Like = GetLogLike(grid(r)) - if (Feedback > 1) write (*,*) r, 'Likelihood: ', Like, 'Current Like:', CurLike + if (Feedback > 1) write (*,*) r, 'Likelihood: ', Like, 'Current Like:', CurLike, 'Rank', MPIRank if ((Like /= logZero) .and. (CurLike > Like .or. randexp1() > Like - CurLike)) then !Accept @@ -501,7 +501,7 @@ function WL_Weight(L) result (W) output_lines = output_lines +1 call WriteParams(CurParams, real(mult), CurLike) end if - if (Feedback > 1) write (*,*) instance, 'Slicing, Current Like:', CurLike + if (Feedback > 1) write (*,*) instance, 'Slicing, Current Like:', CurLike, 'Rank', MPIRank mult = 1 if (num_slow /=0) call SliceSampleSlowParam(CurParams, CurLike) if (num_fast /=0) call SliceSampleFastParams(CurParams, CurLike) @@ -536,7 +536,7 @@ function WL_Weight(L) result (W) Like = GetLogLike(Trial) - if (Feedback > 1) write (*,*) 'Likelihood: ', Like, 'Current Like:', CurLike + if (Feedback > 1) write (*,*) 'Likelihood: ', Like, 'Current Like:', CurLike, 'Rank', MPIRank !!! accpt = (Like /= logZero) .and. (CurLike > Like .or. randexp1() > Like - CurLike) !Include the min() so that compilers not doing optimal compilation don't complain diff --git a/source/Makefile b/source/Makefile index 9f2ad5a..e694dfb 100644 --- a/source/Makefile +++ b/source/Makefile @@ -1,123 +1,84 @@ -#You may need to edit the library paths for MKL for Intel -#Beware of using optimizations that lose accuracy - may give errors when running +# >>> DESIGNED FOR GMAKE <<< -##Uncomment the next line to include dr7 LRG -EXTDATA = -#EXTDATA = LRG +# Unified Systems makefile for COSMOMC +# Add FLAGS -DMPI for using MPI + +ext=$(shell uname | cut -c1-3) + +CC=cc + +ifeq ($(ext),IRI) +F90C= f90 +FFLAGS= -Ofast -mp -n32 -LANG:recursive=ON -lmpi -DMPI +WMAPFLAGS= $(FFLAGS) +LAPACKL = -lcomplib.sgimath_mp +INCLUDE = -I../camb +CFITSIODIR = +GSLDIR = +endif + +ifeq ($(ext),Lin) +F90C=gfortran +FFLAGS= -O -cpp -fopenmp +WMAPFLAGS= -O +LAPACKL = -llapack -lblas +INCLUDE = -I../camb +CFITSIODIR = +GSLDIR = /usr +endif + +ifeq ($(ext),OSF) +F90C=f90 +FFLAGS= -omp -O -arch host -math_library fast -tune host -fpe1 +WMAPFLAGS= $(FFLAGS) +LAPACKL = -lcxml +INCLUDE = -I../camb +CFITSIODIR = +GSLDIR = +endif -#set WMAP empty not to compile with WMAP -WMAP = /home/aml1005/WMAP7/likelihood_v4 - -#Only needed for WMAP -cfitsio = /usr/local/cfitsio/intel10/64/3.040 - -#GSL only needed for DR7 LRG -GSLPATH = /home/aml1005/libs/gsl - -IFLAG = -I -INCLUDE= - -#Intel MPI (assuming mpif77 set to point to ifort) -#these settings for ifort 11.1 and higher; may need to add explicit link directory otherwise -#Can add -xHost if your cluster is uniform, or specify specific processor optimizations -x... -#If getdist gives segfaults remove openmp when compiling getdist -F90C = mpif90 -FFLAGS = -O3 -W0 -WB -openmp -fpp -DMPI -vec_report0 -mkl=parallel -LAPACKL = -lmpi - -#GFortran: defaults for v4.5; if pre v4.3 add -D__GFORTRAN__ -#F90C = gfortran -#in earlier versions use FFLAGS = -O2 -ffree-form -x f95-cpp-input -D__GFORTRAN__ -#FFLAGS = -O3 -fopenmp -ffree-form -x f95-cpp-input -ffast-math -march=native -funroll-loops -#LAPACKL = -Wl,-framework -Wl,accelerate -#commented above is (I think) for Mac; this is standard linux (sudo apt-get install liblapack-dev) -#LAPACKL = -lblas -llapack - -#HPCF settings. Use Intel 9 or 10.1+, not 10.0 -#F90C = mpif90 -#FFLAGS = -O2 -Vaxlib -W0 -WB -openmp -fpp -DMPI -vec_report0 -#LAPACKL = -L/usr/local/Cluster-Apps/intel/mkl/10.2.2.025/lib/em64t -lmkl_lapack -lmkl -lguide -lpthread -#GSLPATH = /usr/local/Cluster-Apps/gsl/1.9 -#cfitsio = /usr/local/Cluster-Users/cpac/cmb/2.1.0/cfitsio -#INCLUDE= - -#COSMOS: use "module load cosmolib latest" -#use "runCosmomc" (globally installed) to run, defining required memory usage -ifeq ($(COSMOHOST),cosmos) -F90C = ifort -FFLAGS = -openmp -fast -w -fpp2 -DMPI -LAPACKL = -mkl=sequential -lmkl_lapack -lmpi -cfitsio = $(CFITSIO) -WMAP = $(COSMOLIB)/WMAP7 -GSLPATH = $(GSL_ROOT) +ifeq ($(ext),Sun) +F90C=f90 +FFLAGS= -O4 -xarch=native64 -openmp -ftrap=%none +WMAPFLAGS= $(FFLAGS) +LAPACKL = -lsunperf -lfsu +INCLUDE = -I../camb -M../camb +CFITSIODIR = +GSLDIR = endif -#Intel fortran 8, check you have the latest update from the Intel web pages -#See Makefile_intel for ifc 7.1 or lower (some versions have problems) -#F90C = ifort -#FFLAGS = -O2 -Vaxlib -ip -W0 -WB -openmp -fpp -#LAPACKL = -L/opt/intel/mkl61/lib/32 -lmkl_lapack -lmkl_ia32 -lguide -lpthread - -#G95; make sure LAPACK and MPI libs also compiled with g95 -#F90C = mpif90 -#FFLAGS = -O2 -cpp -DMPI -#LAPACKL = /LAPACK/lapack_LINUX.a /LAPACK/blas_LINUX.a - -#SGI, -mp toggles multi-processor. Use -O2 if -Ofast gives problems. -#Not various versions of the compiler are buggy giving erroneous seg faults with -mp. -#Version 7.3 is OK, currently version 7.4 is bugged, as are some earlier versions. -#F90C = f90 -#LAPACKL = -lcomplib.sgimath -#FFLAGS = -Ofast -mp - -#Digital/Compaq fortran, -omp toggles multi-processor -#F90C = f90 -#FFLAGS = -omp -O4 -arch host -math_library fast -tune host -fpe1 -#LAPACKL = -lcxml - -#Absoft ProFortran, single processor, set -cpu:[type] for your local system -#F90C = f95 -#FFLAGS = -O2 -s -cpu:athlon -lU77 -w -YEXT_NAMES="LCS" -YEXT_SFX="_" -#LAPACKL = -llapack -lblas -lg2c -#IFLAG = -p - -#NAGF95, single processor: -#F90C = f95 -#FFLAGS = -DNAGF95 -O3 -#LAPACKL = -llapack -lblas -lg2c - -#PGF90 -#F90C = pgf90 -#FFLAGS = -O2 -DESCAPEBACKSLASH -Mpreprocess -#LAPACKL = -llapack -lblas - - -#Sun, single processor: -#F90C = f90 -#FFLAGS = -fast -ftrap=%none -#LAPACKL = -lsunperf -lfsu -#LAPACKL = -lsunperf -lfsu -lsocket -lm -#IFLAG = -M - -#Sun MPI -#F90C = mpf90 -#FFLAGS = -O4 -openmp -ftrap=%none -dalign -DMPI -#LAPACKL = -lsunperf -lfsu -lmpi_mt -#IFLAG = -M - -#Sun parallel enterprise: -#F90C = f95 -#FFLAGS = -O4 -xarch=native64 -openmp -ftrap=%none -#LAPACKL = -lsunperf -lfsu -#IFLAG = -M - - -#IBM XL Fortran, multi-processor (run "module load lapack" then run "gmake") -# See also http://cosmocoffee.info/viewtopic.php?t=326 -#F90C = xlf90_r $(LAPACK) -#FFLAGS = -WF,-DIBMXL -qsmp=omp -qsuffix=f=f90:cpp=F90 -O3 -qstrict -qarch=pwr3 -qtune=pwr3 -#INCLUDE = -lessl -#LAPACKL = +ifeq ($(ext),AIX) +F90C = mpxlf90_r +FFLAGS = -O4 -WF,-DIBMXL,-DMPI -qstrict -qsmp=omp -qmaxmem=-1 -qsuffix=f=f90:cpp=F90 +WMAPFLAGS= $(FFLAGS) +LAPACKL = -lessl +INCLUDE = -I../camb +CFITSIODIR = +GSLDIR = +endif + +#EXTDATA = LRG +EXTDATA = LRG +EXTINCLUDE = -I$(GSLDIR)/include +EXTOBJS = bsplinepk.o + + +WMAPDIR = ../WMAP +WMAPINCLUDE = -I$(CFITSIODIR)/include +WMAPOBJS = read_archive_map.o \ + read_fits.o \ + healpix_types.o \ + br_mod_dist.o \ + WMAP_7yr_options.o \ + WMAP_7yr_util.o \ + WMAP_7yr_gibbs.o \ + WMAP_7yr_tt_pixlike.o \ + WMAP_7yr_tt_beam_ptsrc_chisq.o \ + WMAP_7yr_teeebb_pixlike.o \ + WMAP_7yr_tetbeebbeb_pixlike.o \ + WMAP_7yr_likelihood.o + + PROPOSE = propose.o CLSFILE = CMB_Cls_simple.o @@ -125,51 +86,44 @@ CLSFILE = CMB_Cls_simple.o #Can use params_H if you prefer more generic parameters PARAMETERIZATION = params_CMB.o -F90FLAGS = -DMATRIX_SINGLE $(FFLAGS) $(IFLAG)../camb $(INCLUDE) -LINKFLAGS = -L../camb -lcamb $(LAPACKL) +LINKFLAGS = -L../camb -lcamb $(LAPACKL) + +F90FLAGS = -DMATRIX_SINGLE $(FFLAGS) $(INCLUDE) -DISTFILES = ParamNames.o Matrix_utils.o settings.o IO.o GetDist.o +DISTFILES = utils.o ParamNames.o Matrix_utils.o settings.o IO.o GetDist.o -OBJFILES= ParamNames.o Matrix_utils.o settings.o IO.o cmbtypes.o Planck_like.o \ - cmbdata.o WeakLen.o bbn.o lrggettheory.o mpk.o bao.o supernovae.o HST.o SDSSLy-a-v3.o \ + + +OBJFILES = utils.o ParamNames.o Matrix_utils.o settings.o IO.o cmbtypes.o Planck_like.o \ + cmbdata.o WeakLen.o bbn.o lrggettheory.o mpk.o bao.o supernovae.o HST.o SDSSLy-a-v3.o\ $(CLSFILE) paramdef.o $(PROPOSE) $(PARAMETERIZATION) calclike.o \ conjgrad_wrapper.o EstCovmat.o postprocess.o MCMC.o driver.o - -ifeq ($(EXTDATA),LRG) +ifeq ($(EXTDATA),LRG) F90FLAGS += -DDR71RG -OBJFILES += bsplinepk.o -GSLINC = -I$(GSLPATH)/include -LINKFLAGS += -L$(GSLPATH)/lib -lgsl -lgslcblas -endif - -ifneq ($(WMAP),) -F90FLAGS += $(IFLAG)$(cfitsio)/include $(IFLAG)$(WMAP) -LINKFLAGS += -L$(cfitsio)/lib -L$(WMAP) -lcfitsio - -OBJFILES += $(WMAP)/read_archive_map.o \ - $(WMAP)/read_fits.o \ - $(WMAP)/healpix_types.o \ - $(WMAP)/WMAP_7yr_options.o \ - $(WMAP)/WMAP_7yr_util.o \ - $(WMAP)/WMAP_7yr_tt_pixlike.o \ - $(WMAP)/WMAP_7yr_teeebb_pixlike.o \ - $(WMAP)/WMAP_7yr_likelihood.o \ - $(WMAP)/WMAP_7yr_gibbs.o \ - $(WMAP)/WMAP_7yr_tt_beam_ptsrc_chisq.o \ - $(WMAP)/br_mod_dist.o +EXTINCLUDE = -I$(GSLDIR)/include +LINKFLAGS += -lgsl -lgslcblas +OBJFILES += $(EXTOBJS) +endif + +ifneq ($(WMAPDIR),) +F90FLAGS += $(WMAPINCLUDE) +LINKFLAGS += -lcfitsio +OBJFILES += $(WMAPOBJS) else F90FLAGS += -DNOWMAP endif -default: cosmomc -all : cosmomc getdist -settings.o: ../camb/libcamb.a +default: cosmomc.$(ext) + +all : cosmomc.$(ext) getdist.$(ext) + +settings.o: utils.o cmbtypes.o: settings.o Planck_like.o: cmbtypes.o -cmbdata.o: Planck_like.o $(WMAPOBJS) +cmbdata.o: Planck_like.o $(WMAPOBJS) WeakLen.o: cmbtypes.o bbn.o: settings.o mpk.o: cmbtypes.o lrggettheory.o @@ -188,14 +142,20 @@ postprocess.o: calclike.o MCMC.o: calclike.o driver.o: MCMC.o -export FFLAGS -export F90C - .f.o: f77 $(F90FLAGS) -c $< %.o: %.c - $(CC) $(GSLINC) -c $*.c + $(CC) $(EXTINCLUDE) -c $*.c + +%.o: $(WMAPDIR)/%.f90 + $(F90C) $(WMAPFLAGS) $(WMAPINCLUDE) -c $< + +%.o: $(WMAPDIR)/%.F90 + $(F90C) $(WMAPFLAGS) $(WMAPINCLUDE) -c $< + +utils.o: ../camb/utils.F90 + $(F90C) $(F90FLAGS) -c $< %.o: %.f90 $(F90C) $(F90FLAGS) -c $*.f90 @@ -204,19 +164,12 @@ export F90C $(F90C) $(F90FLAGS) -c $*.F90 -cosmomc: camb $(OBJFILES) - $(F90C) -o ../cosmomc $(OBJFILES) $(LINKFLAGS) $(F90FLAGS) +cosmomc.$(ext): $(OBJFILES) ../camb/libcamb.a + $(F90C) -o ../$@ $(OBJFILES) $(LINKFLAGS) $(F90FLAGS) -clean: cleancosmomc - rm -f ../camb/*.o ../camb/*.obj ../camb/*.mod - -cleancosmomc: +clean: rm -f *.o *.mod *.d *.pc *.obj ../core - -getdist: camb $(DISTFILES) - $(F90C) -o ../getdist $(DISTFILES) $(LINKFLAGS) $(F90FLAGS) - -camb: - cd ../camb && $(MAKE) --file=Makefile_main libcamb.a +getdist.$(ext): $(DISTFILES) + $(F90C) -o ../$@ $(DISTFILES) $(LINKFLAGS) $(F90FLAGS) diff --git a/source/driver.F90 b/source/driver.F90 index 8188b94..a423f68 100644 --- a/source/driver.F90 +++ b/source/driver.F90 @@ -59,11 +59,17 @@ program SolveCosmology #ifdef MPI - if (instance /= 0) call DoAbort('With MPI should not have second parameter') +!wig if (instance /= 0) call DoAbort('With MPI should not have second parameter') + instance_shift=instance + write(*,*)'MPI file names shifted by: ',instance_shift +!end wig call mpi_comm_rank(mpi_comm_world,MPIrank,ierror) - + instance = MPIrank +1 !start at 1 for chains - write (numstr,*) instance +!wig write (numstr,*) instance + instance_shift = MPIrank + instance_shift + write (numstr,*) instance_shift +!end wig rand_inst = instance if (ierror/=MPI_SUCCESS) call DoAbort('MPI fail') @@ -104,6 +110,13 @@ program SolveCosmology HighAccuracyDefault = Ini_Read_Logical('high_accuracy_default',.false.) AccuracyLevel = Ini_Read_Real('accuracy_level',1.) +!wig + lAccuracyLevel = Ini_Read_Real('l_accuracy_level',1.) + lSampleLevel = Ini_Read_Real('l_sample_level',1.) + kSrcSampleLevel = Ini_Read_Real('k_src_sample_level',1.) + kIntSampleLevel = Ini_Read_Real('k_int_sample_level',1.) + CambComputeEveryl = Ini_Read_Logical('camb_compute_every_l',.false.) +!end wig if (Ini_HasKey('highL_unlensed_cl_template')) & highL_unlensed_cl_template= ReadIniFilename(DefIni,'highL_unlensed_cl_template') diff --git a/source/params_CMB.f90 b/source/params_CMB.f90 index a453def..0634a8a 100644 --- a/source/params_CMB.f90 +++ b/source/params_CMB.f90 @@ -57,24 +57,46 @@ integer, intent(in) :: in - if (Power_Name == 'power_tilt') then - - P%InitPower%k_0_scalar = pivot_k - P%InitPower%k_0_tensor = pivot_k +!wig + if (Power_Name == 'power_sr') then + + P%InitPower%P_scalar = cl_norm*CMB%norm(norm_As) + P%InitPower%eps1(in) = CMB%InitPower(1) + P%InitPower%eps2(in) = CMB%InitPower(2) + P%InitPower%eps3(in) = CMB%InitPower(3) + + P%InitPower%logeps1invsigma0(in) = CMB%InitPower(5) + P%InitPower%loginvsigma0(in) = P%InitPower%logeps1invsigma0(in) & + - log10(P%InitPower%eps1(in)) + P%InitPower%invsigma0(in) = 10**(P%InitPower%loginvsigma0(in)) + P%InitPower%Xsigma0(in) = CMB%InitPower(6) + P%InitPower%varphi(in) = CMB%InitPower(7) + !CMB%InitPower(7) unused. + CMB%InitPower(7) = 1./P%InitPower%invsigma0(in) + + elseif (Power_Name == 'power_dbi') then + + P%InitPower%P_scalar = cl_norm*CMB%norm(norm_As) + P%InitPower%eps1(in) = 10.**(CMB%InitPower(4)) + P%InitPower%del1(in) = P%InitPower%eps1(in) - CMB%InitPower(1) + P%InitPower%eps2(in) = CMB%InitPower(2) - P%InitPower%del1(in) + P%InitPower%eps3(in) = CMB%InitPower(3) + + P%InitPower%logeps1invsigma0(in) = CMB%InitPower(5) + P%InitPower%loginvsigma0(in) = P%InitPower%logeps1invsigma0(in) & + - log10(P%InitPower%eps1(in)) + P%InitPower%invsigma0(in) = 10**(P%InitPower%loginvsigma0(in)) + P%InitPower%Xsigma0(in) = CMB%InitPower(6) + P%InitPower%varphi(in) = CMB%InitPower(7) + !CMB%InitPower(7) unused. + CMB%InitPower(7) = 1./P%InitPower%invsigma0(in) + + P%InitPower%gam(in) = 10**(-CMB%InitPower(8))*P%InitPower%eps1(in) - P%InitPower%ScalarPowerAmp(in) = cl_norm*CMB%norm(norm_As) - P%InitPower%rat(in) = CMB%norm(norm_amp_ratio) - - P%InitPower%an(in) = CMB%InitPower(1) - P%InitPower%ant(in) = CMB%InitPower(2) - P%InitPower%n_run(in) = CMB%InitPower(3) - if (inflation_consistency) then - P%InitPower%ant(in) = - CMB%norm(norm_amp_ratio)/8. - !note input n_T is ignored, so should be fixed (to anything) - end if else stop 'params_CMB:Wrong initial power spectrum' end if +!end wig end subroutine SetCAMBInitPower diff --git a/source/settings.f90 b/source/settings.f90 index 90d84d0..0998c93 100644 --- a/source/settings.f90 +++ b/source/settings.f90 @@ -6,6 +6,13 @@ module settings implicit none real :: AccuracyLevel = 1. +!wig + real :: lAccuracyLevel = 1. + real :: lSampleLevel = 1. + real :: kSrcSampleLevel = 1. + real :: kIntSampleLevel = 1. + logical :: CambComputeEveryl = .false. +!end wig !Set to >1 to use CAMB etc on higher accuracy settings. !Does not affect MCMC (except making it all slower) @@ -19,7 +26,10 @@ module settings ! (e.g. beam uncertainty modes, etc, specific to dataset) integer, parameter :: num_hard =7 - integer, parameter :: num_initpower = 3 +!wig +! integer, parameter :: num_initpower = 3 + integer, parameter :: num_initpower = 8 +!end wig integer, parameter :: num_freq_params = 1 integer, parameter :: num_norm = 2 + num_freq_params integer, parameter :: num_nuisance_params= 0 @@ -66,6 +76,9 @@ module settings integer :: num_threads = 0 integer :: instance = 0 +!wig + integer :: instance_shift = 0 +!end wig integer :: MPIchains = 1, MPIrank = 0 logical :: Use_LSS = .true.