diff --git a/README b/README new file mode 100644 index 0000000..6c059c2 --- /dev/null +++ b/README @@ -0,0 +1,67 @@ +********************************************************************** +Adds fingerprints features and modulations to power law primordial +power spectra as in arXigv:1106.1635 +********************************************************************** + +-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 "feat" 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_tilt.f90 +Implementation of fingerprints'features. + +- 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/driver.F90 +2nd instance in mpirun shifts the file names. + +- cosmomc/params.ini +New keys for the above material. + + + diff --git a/camb/Makefile b/camb/Makefile index 5c09a67..856b7ee 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_tilt +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..2fb45ac 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 +!feat +! 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 feat 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 +!feat +! 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 feat if (HighAccuracyDefault) dkn2=dkn2/1.2 if (CP%WantTensors .or. CP%WantVectors) then dkn1=dkn1 *0.8_dl @@ -1082,8 +1089,11 @@ contains qmax_int = min(qmax,max_bessels_etak/CP%tau0) +!feat +! IntSampleBoost=AccuracyBoost + IntSampleBoost=AccuracyBoost*kIntSampleBoost +!end feat - IntSampleBoost=AccuracyBoost 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..db94141 100644 --- a/camb/inidriver.F90 +++ b/camb/inidriver.F90 @@ -281,6 +281,19 @@ write(*,*) 'Output _params.ini not created as would overwrite input' end if end if +!feat + + 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 feat call Ini_Close @@ -295,11 +308,28 @@ write (*,'("Age of universe/GYr = ",f7.3)') Age end if +!feat + do i=1,P%InitPower%nn + write (*,'("Scalar tilt = ",f12.8)') P%InitPower%an(i) + write (*,'("Scalar running = ",f12.8)') P%InitPower%n_run(i) + if (P%WantTensors) then + write (*,'("primordial T/S = ",f12.8)') P%InitPower%rat(i) + write (*,'("tensor tilt = ",f12.8)') P%InitPower%ant(i) + endif + write (*,'("Feature pivot = ",f12.8)') P%InitPower%krhalf(i) + write (*,'("Feature size = ",f12.8)') P%InitPower%size(i) + write (*,'("Feature frequency = ",f14.8)') P%InitPower%omega(i) + write (*,'("Feature phase = ",f12.8)') P%InitPower%phase(i) + write (*,'("Feature power = ",f12.8)') P%InitPower%p(i) + enddo +!end feat + if (global_error_flag==0) call CAMB_GetResults(P) if (global_error_flag/=0) then write (*,*) 'Error result '//trim(global_error_message) stop endif + if (P%WantTransfer .and. .not. (P%NonLinear==NonLinear_lens .and. P%DoLensing)) then call Transfer_SaveToFiles(MT,TransferFileNames) diff --git a/camb/modules.f90 b/camb/modules.f90 index 6053b4c..e79258e 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 +!feat +! logical, parameter :: DebugMsgs=.false. !Set to true to view progress and timing + logical, parameter :: DebugMsgs=.false. +!end feat 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 +!feat + real(dl) :: AccuracyBoost =8._dl +!end feat !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. +!feat + real(sp) :: lAccuracyBoost=2._dl +!end feat !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 + +!feat + real(dl) :: kSrcSampleBoost=1._dl + real(dl) :: kIntSampleBoost=1._dl + logical :: DoClsInterpol=.true. + logical :: ComputeEveryl=.true. +!end feat 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 +!feat + 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 feat + Ascale=scale/lSampleBoost if (lSampleBoost >=50) then @@ -720,6 +750,9 @@ end if lSet%l0=lind lSet%l(1:lind) = ls(1:lind) +!feat + endif +!end feat 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)) +!feat + 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 feat + end subroutine InterpolateClArr subroutine InterpolateClArrTemplated(lSet,iCl, all_Cl, max_ind, template_index) @@ -2038,9 +2088,16 @@ minkh = 1e-4 - dlnkh = 0.02 +!feat +! dlnkh = 0.02 + dlnkh = 1e-3 +!end feat 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) +!feat + write(*,*)'Transfert_SaveToFiles: dlnkh= points= ',dlnkh,points + write(*,*)'You should check with transfer_k_per_logint > 0' +!end feat 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..833657d 100644 --- a/camb/params.ini +++ b/camb/params.ini @@ -1,13 +1,13 @@ #Parameters for CAMB #output_root is prefixed to output file names -output_root = test +output_root = feature #What to do get_scalar_cls = T get_vector_cls = F get_tensor_cls = F -get_transfer = F +get_transfer = T #if do_lensing then scalar_output_file contains additional columns of l^4 C_l^{pp} and l^3 C_l^{pT} #where p is the projected potential. Output lensed CMB Cls (without tensors) are in lensed_output_file below. @@ -77,6 +77,12 @@ tensor_spectral_index(1) = 0 initial_ratio(1) = 1 #note vector modes use the scalar settings above +#Primordial features +feature_kscale(1) = 0.05 +feature_size(1) = 1. +feature_power(1) = 0.5 +feature_pulsation(1) = 50 +feature_phase(1) = 0 #Reionization, ignored unless reionization = T, re_redshift measures where x_e=0.5 reionization = T @@ -207,8 +213,8 @@ number_of_threads = 0 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 +#Decrease to speed up at cost of worse accuracy. Suggest 3 to 8 for features. +accuracy_boost = 2 #Larger to keep more terms in the hierarchy evolution. l_accuracy_boost = 1 @@ -218,3 +224,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_tilt.f90 b/camb/power_tilt.f90 index c23155f..d1ace2b 100644 --- a/camb/power_tilt.f90 +++ b/camb/power_tilt.f90 @@ -35,7 +35,9 @@ character(LEN=*), parameter :: Power_Name = 'power_tilt' - integer, parameter :: nnmax= 5 +!features + integer, parameter :: nnmax= 1 +!end features !Maximum possible number of different power spectra to use Type InitialPowerParams @@ -50,7 +52,16 @@ 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) - +!features + real(dl) :: krhalf(nnmax) !=kr/2 + real(dl) :: size(nnmax) !=beta x model dep factor + real(dl) :: phase(nnmax) != phi + p^2/(p-1) 2 M/H + real(dl) :: omega(nnmax) != 2 M/H + real(dl) :: p(nnmax) + + +!end features + end Type InitialPowerParams real(dl) curv !Curvature contant, set in InitializePowers @@ -76,6 +87,13 @@ AP%k_0_scalar = 0.05 AP%k_0_tensor = 0.05 AP%ScalarPowerAmp = 1 +!feat + AP%krhalf = 0.07 + AP%size = 0. + AP%phase = 0. + AP%omega = 10 + AP%p = 10 +!end feat end subroutine SetDefPowerParams @@ -116,14 +134,42 @@ !< |\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 + real(dl) ScalarPower,k, lnrat, lnf + real(dl) :: angle, amp + integer in - + lnrat = log(k/P%k_0_scalar) - ScalarPower=P%ScalarPowerAmp(in)*exp((P%an(in)-1)*lnrat + P%n_run(in)/2*lnrat**2) - -! ScalarPower = ScalarPower * (1 + 0.1*cos( lnrat*30 ) ) - +!features + lnf = log(k/P%krhalf(in)) + amp = 0._dl + + if ((P%p(in).eq.1._dl).or.(abs(lnf/P%p(in)).gt.log(huge(1._dl)))) then + + angle = 0._dl + + else + + angle = P%phase(in) + P%omega(in) * (P%p(in)/(P%p(in)-1._dl)) & + *( P%p(in)*exp(lnf/P%p(in))- P%p(in) ) + + if ((P%p(in).gt.0._dl).and.(P%p(in).lt.1._dl)) then + if (k.le.P%krhalf(in)) amp = exp(lnf*(-3._dl + 2.5_dl/P%p(in))) * P%size(in) + else + if (k.ge.P%krhalf(in)) amp = exp(lnf*(-3._dl + 2.5_dl/P%p(in))) * P%size(in) + endif + + endif + + ScalarPower=P%ScalarPowerAmp(in)*exp((P%an(in)-1)*lnrat + P%n_run(in)/2*lnrat**2) & + * (1._dl + amp*sin(angle)) + + + if(ScalarPower.lt.0.) then + ScalarPower=0._dl + end if +!end features + end function ScalarPower @@ -168,6 +214,23 @@ num=num+1 Keys(num) = 's_pivot' Vals(num) = P%k_0_scalar +!features + num=num+1 + Keys(num) = 'f_pivot' + Vals(num) = P%krhalf(in) + num=num+1 + Keys(num) = 'size' + Vals(num) = P%size(in) + num=num+1 + Keys(num) = 'omega' + Vals(num) = P%omega(in) + num=num+1 + Keys(num) = 'power' + Vals(num) = P%p(in) + num=num+1 + Keys(num) = 'phase' + Vals(num) = P%phase(in) +!end features end if if (Tens) then num=num+1 @@ -202,7 +265,13 @@ InitPower%an(i) = Ini_Read_Double_Array_File(Ini,'scalar_spectral_index', i) InitPower%n_run(i) = Ini_Read_Double_Array_File(Ini,'scalar_nrun',i,0._dl) - +!features + InitPower%krhalf(i) = Ini_Read_Double_Array_File(Ini,'feature_kscale',i,0.02_dl) + InitPower%size(i) = Ini_Read_Double_Array_File(Ini,'feature_size',i, 1._dl) + InitPower%omega(i) = Ini_Read_Double_Array_File(Ini,'feature_pulsation', i) + InitPower%phase(i) = Ini_Read_Double_Array_File(Ini,'feature_phase',i,0._dl) + InitPower%p(i) = Ini_Read_Double_Array_File(Ini,'feature_power', i) +!end features if (WantTensors) then InitPower%ant(i) = Ini_Read_Double_Array_File(Ini,'tensor_spectral_index',i) InitPower%rat(i) = Ini_Read_Double_Array_File(Ini,'initial_ratio',i) 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/params.ini b/params.ini index 6e524dd..887ee86 100644 --- a/params.ini +++ b/params.ini @@ -1,7 +1,7 @@ #Sample parameters for cosmomc in default parameterization #Root name for files produced -file_root = chains/test +file_root = chains/feat #action = 0: MCMC, action=1: postprocess .data file, action=2: find best fit point only action = 0 @@ -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 @@ -140,7 +140,7 @@ rand_seed = #If true, generate checkpoint files and terminated runs can be restarted using exactly the same command #and chains continued from where they stopped #With checkpoint=T note you must delete all chains/file_root.* files if you want new chains with an old file_root -checkpoint = F +checkpoint = T #whether to stop on CAMB error, or continue ignoring point stop_on_error= T @@ -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 = 2 +l_accuracy_level = 1 +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 @@ -196,8 +201,14 @@ 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 -#log[10^10 A_s] -param[logA] = 3 2.7 4 0.01 0.01 -param[r] = 0 0 0 0 0 +param[logkf] = -1.3 -1.3 -1.3 0 0 +param[size] = 0 0 0 0 0 +param[pow] = 0 0 0 0 0 +param[freq] = 0 0 0 0 0 +param[phase] = 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..895d4fa 100644 --- a/params_CMB.paramnames +++ b/params_CMB.paramnames @@ -8,7 +8,12 @@ w w #constant equation of state parameter for scalar field dark en ns n_s #beware that pivot scale can change in .ini file nt n_t nrun n_{run} -logA log[10^{10} A_s] +logkf log(k_r/2) #feature scale +size A_{\omega} #feature size +pow p #power dependency +freq \omega #oscillation frequency +phase \psi #phase shift +lnA ln[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 omegal* \Omega_\Lambda @@ -17,4 +22,4 @@ 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 +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..e3659ed 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) +!feat + 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 feat + 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 +!feat +! AccuracyBoost = AccuracyLevel +! lAccuracyBoost = AccuracyLevel +! lSampleBoost = AccuracyLevel +!end feat P%AccurateReionization = .true. end if 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..6f5f333 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,43 @@ 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) -DISTFILES = ParamNames.o Matrix_utils.o settings.o IO.o GetDist.o +F90FLAGS = -DMATRIX_SINGLE $(FFLAGS) $(INCLUDE) -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 \ +DISTFILES = utils.o ParamNames.o Matrix_utils.o settings.o IO.o GetDist.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 + +default: cosmomc.$(ext) + +all : cosmomc.$(ext) getdist.$(ext) settings.o: ../camb/libcamb.a 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 +141,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 +163,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..84b2de8 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') +!feat if (instance /= 0) call DoAbort('With MPI should not have second parameter') + instance_shift=instance + write(*,*)'MPI file names shifted by: ',instance_shift +!end feat call mpi_comm_rank(mpi_comm_world,MPIrank,ierror) - + instance = MPIrank +1 !start at 1 for chains - write (numstr,*) instance +!feat write (numstr,*) instance + instance_shift = MPIrank + instance_shift + write (numstr,*) instance_shift +!end feat 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.) +!feat + 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 feat 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..d90d3bd 100644 --- a/source/params_CMB.f90 +++ b/source/params_CMB.f90 @@ -68,6 +68,13 @@ P%InitPower%an(in) = CMB%InitPower(1) P%InitPower%ant(in) = CMB%InitPower(2) P%InitPower%n_run(in) = CMB%InitPower(3) +!feat + P%InitPower%krhalf(in) = 10._dl**CMB%InitPower(4) + P%InitPower%size(in) = CMB%InitPower(5) + P%InitPower%p(in) = CMB%InitPower(6) + P%InitPower%omega(in) = CMB%InitPower(7) + P%InitPower%phase(in) = CMB%InitPower(8) +!end feat 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) diff --git a/source/settings.f90 b/source/settings.f90 index 90d84d0..deab0c4 100644 --- a/source/settings.f90 +++ b/source/settings.f90 @@ -6,6 +6,13 @@ module settings implicit none real :: AccuracyLevel = 1. +!feat + real :: lAccuracyLevel = 1. + real :: lSampleLevel = 1. + real :: kSrcSampleLevel = 1. + real :: kIntSampleLevel = 1. + logical :: CambComputeEveryl = .false. +!end feat !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 +!feat +! integer, parameter :: num_initpower = 3 + integer, parameter :: num_initpower = 8 +!end feat 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 +!feat + integer :: instance_shift = 0 +!end feat integer :: MPIchains = 1, MPIrank = 0 logical :: Use_LSS = .true.