Merge pull request #6 from joezuntz/v2-demos

Fix the demos and their ingredients

.github/run-demos.sh

@@ -0,0 +1,71 @@
+#!/usr/bin/env bash
+
+set -e
+
+# Bash is insane. If you don't do this next bit then it
+# decides that you wanted to pass all the arguments from this
+# script on to cosmosis-configure. Madness.
+function DoSource() { source cosmosis-configure ; } ; DoSource
+
+
+# for demo 11 first run
+export HALOFIT=takahashi
+
+if [ "$1" == "" ]; then
+    demos=$(seq 1 19)
+else
+    demos=$1
+fi
+
+echo "Running demos: " $demos
+
+for i in $demos
+do
+    if [ "$i" == "3" ]; then
+        args="-p grid.nsample_dimension=10"
+    elif [ "$i" == "4" ]; then
+        args="-v cosmological_parameters.n_s=0.962 cosmological_parameters.h0=0.680 -p maxlike.tolerance=1.0"
+    elif [ "$i" == "5" ]; then
+        args="-p emcee.samples=30 emcee.nsteps=10 emcee.walkers=20"
+    elif [ "$i" == "9" ]; then
+        args="-p multinest.live_points=50"
+    elif [ "$i" == "10" ]; then
+        args="-p grid.nsample_dimension=10"
+    elif [ "$i" == "11" ]; then
+        args="-p grid.nsample_dimension=10"
+    elif [ "$i" == "13" ]; then
+        args="-p snake.threshold=3"
+    elif [ "$i" == "14" ]; then
+        args="-p zeus.samples=15 zeus.nsteps=5"
+    elif [ "$i" == "17" ]; then
+        args="-v cosmological_parameters.n_s=0.96 cosmological_parameters.omega_b=0.0468 cosmological_parameters.h0=0.6881"
+    elif [ "$i" == "18" ]; then
+        echo "Skipping demo 18"
+        continue
+    elif [ "$i" == "19" ]; then
+        args="-p grid.nsample_dimension=20"
+    else
+        args=
+    fi
+    echo Running demo $i:  cosmosis demos/demo${i}.ini $args
+    time cosmosis demos/demo${i}.ini $args
+    echo Done demo $i
+    echo ""
+    echo ""
+    echo "Postprocessing demo $i" 
+    if [ "$i" == "10" ]; then
+        # Run the second part of demo 10
+        export HALOFIT=mead2020
+        echo "Running demo 10 (part 2):  cosmosis demos/demo${i}.ini $args"
+        time cosmosis demos/demo${i}.ini $args
+
+        # Postprocess demo10
+        cosmosis-postprocess output/demo10_mead2020.txt output/demo10_takahashi.txt -o output/plots/demo10
+    else
+        cosmosis-postprocess demos/demo${i}.ini  -o output/plots/demo${i}
+    fi
+    echo ""
+    echo ""
+    echo ""
+    echo ""
+done

.github/workflows/ci.yml

@@ -69,6 +69,14 @@ jobs:
           cosmosis examples/des-y3.ini | tee output/des-y3.log
           grep 'Likelihood =  6043.23' output/des-y3.log
 
+      - name: Run Demos
+        shell: bash -l {0}
+        run: |
+          .github/run-demos.sh
+
+
+
+
   apt-get-test:
     runs-on: ubuntu-latest
     strategy:

boltzmann/camb/camb_interface.py

@@ -451,7 +451,7 @@ def compute_growth_rates(r, block, P_tot, k, z, more_config):
     kmin = k.min()
     P_kmin = P_tot(z, kmin)
 
-    D = np.sqrt(P_kmin / P_kmin[0])
+    D = np.sqrt(P_kmin / P_kmin[0]).squeeze()
     a = 1/(1+z)
 
     loga = np.log(a[::-1])
@@ -461,7 +461,7 @@ def compute_growth_rates(r, block, P_tot, k, z, more_config):
     # since z is increasing so a is decreasing.
     logD_spline = InterpolatedUnivariateSpline(loga, np.log(D[::-1]) )
     f_spline = logD_spline.derivative()
-    f = f_spline(loga)[::-1]
+    f = f_spline(loga)[::-1].squeeze()
 
     return f, D, a
 

boltzmann/extrapolate/extrapolate_power.py

@@ -31,14 +31,14 @@ def linear_extend(x, y, xmin, xmax, nmin, nmax, nfit):
     return x, y
 
 
-def extrapolate_section(block, section, kmin, kmax, nmin, nmax, npoint):
+def extrapolate_section(block, section, kmin, kmax, nmin, nmax, npoint, key):
     # load current values
     k = block[section, "k_h"]
     z = block[section, "z"]
     nk = len(k)
     nz = len(z)
     # load other current values
-    k, z, P = block.get_grid(section, "k_h", "z", "p_k")
+    k, z, P = block.get_grid(section, "k_h", "z", key)
     # extrapolate
     P_out = []
     for i in range(nz):
@@ -50,7 +50,7 @@ def extrapolate_section(block, section, kmin, kmax, nmin, nmax, npoint):
     k = exp(logk)
     P_out = np.dstack(P_out).squeeze()
 
-    block.replace_grid(section, "z", z, "k_h", k, "P_k", P_out.T)
+    block.replace_grid(section, "z", z, "k_h", k, key, P_out.T)
 
 
 def setup(options):
@@ -76,6 +76,10 @@ def execute(block, config):
 
     # extrapolate non-linear power
     for section in [names.matter_power_nl, names.matter_power_lin]+extrapk:
+        if '.' in section:
+            section, key = section.split('.')
+        else:
+            key = 'p_k'
         if block.has_section(section):
-            extrapolate_section(block, section, kmin, kmax, nmin, nmax, npoint)
+            extrapolate_section(block, section, kmin, kmax, nmin, nmax, npoint, key)
     return 0

boltzmann/isitgr/Makefile

@@ -1,7 +1,7 @@
 include ${COSMOSIS_SRC_DIR}/config/compilers.mk
 
 USER_FFLAGS=-I.
-USER_LDFLAGS=-lcosmosis
+USER_LDFLAGS=-lcosmosis -Wl,-rpath,$(PWD)/camb_Jan12_isitgr
 
 export EQUATIONS ?= equations
 export RECOMBINATION ?= recfast

boltzmann/isitgr/camb_interface.F90

@@ -457,26 +457,31 @@ function camb_interface_save_sigma8(block) result(status)
 	end function
 
 	function camb_interface_save_transfer(block) result(status)
+		use nonlinear
 		integer (cosmosis_block) :: block
 		integer (cosmosis_status) :: status
-		Type(MatterPowerData) :: PK
+		Type(MatterPowerData) :: PK, PK_nl
 		integer nz, nk, iz, ik
 		external dtauda
 		real(8) dtauda
 		real(8), allocatable, dimension(:) :: k, z
-		real(8), allocatable, dimension(:,:) :: P, T
+		real(8), allocatable, dimension(:,:) :: P, P_nl, T
 		real(8), allocatable, dimension(:,:) :: ModifiedGravity_D, ModifiedGravity_D_dot
 		real(8), allocatable, dimension(:,:) :: ModifiedGravity_Q, ModifiedGravity_Q_dot
 		real(8) a, adotoa, TGR_D_T, TGR_D_T_dot, TGR_Q_T, TGR_Q_T_dot
 
 		call Transfer_GetMatterPowerData(MT, PK, 1)
+		call Transfer_GetMatterPowerData(MT, PK_nl, 1)
+		call NonLinear_GetNonLinRatios(PK_nl)
+		call MatterPowerdata_MakeNonlinear(PK_nl)
 
 		nz = CP%Transfer%num_redshifts
 		nk = MT%num_q_trans
 
 		allocate(k(nk))
 		allocate(z(nz))
 		allocate(P(nk,nz))
+		allocate(P_nl(nk,nz))
 		allocate(ModifiedGravity_D(nk,nz))
 		allocate(ModifiedGravity_Q(nk,nz))
 		allocate(ModifiedGravity_D_dot(nk,nz))
@@ -494,6 +499,7 @@ function camb_interface_save_transfer(block) result(status)
 		do ik=1,nk
 			do iz=1,nz
 				P(ik,iz) = MatterPowerData_k(PK, k(ik), nz-iz+1)
+				P_nl(ik,iz) = MatterPowerData_k(PK_nl, k(ik), nz-iz+1)
 				T(ik,iz) = MT%TransferData(Transfer_cdm,ik,nz-iz+1)
 				!Modifications for IsItGR
 				a = 1.0/(1+z(iz))
@@ -506,9 +512,16 @@ function camb_interface_save_transfer(block) result(status)
 			enddo
 		enddo
 
+		call MatterPowerdata_Free(PK)
+		call MatterPowerdata_Free(PK_nl)
+
 		status = datablock_put_double_grid(block, matter_power_lin_section, &
         	"k_h", k, "z", z, "P_k", P)
 
+
+		status = datablock_put_double_grid(block, matter_power_nl_section, &
+        	"k_h", k, "z", z, "P_k", P_NL)
+
 		if (status .ne. 0) then
 			write(*,*) "Failed to save matter power in CAMB."
 		endif
@@ -531,7 +544,7 @@ function camb_interface_save_transfer(block) result(status)
 		endif
 
 
-		deallocate(k, z, P, T)
+		deallocate(k, z, P, P_NL, T)
 	end function
 
 

boltzmann/mgcamb/Makefile

@@ -1,7 +1,7 @@
 include ${COSMOSIS_SRC_DIR}/config/compilers.mk
 
 USER_FFLAGS=-I.
-USER_LDFLAGS=-lcosmosis
+USER_LDFLAGS=-lcosmosis -Wl,-rpath,$(PWD)/camb_Jan12_mgcamb
 
 export EQUATIONS ?= equations
 export RECOMBINATION ?= recfast

demos/demo1.ini

@@ -6,7 +6,7 @@ root = ${PWD}
 [test]
 ; These are the parameters for this sampler.
 ; In this case there is just one parameter
-save_dir=demo_output_1
+save_dir=output/demo1
 
 ; If something goes wrong this will aid diagnosis
 fatal_errors=T

demos/demo10.ini

@@ -13,7 +13,7 @@ nsample_dimension = 60
 [pipeline]
 ; You can use environment variables in these ini files.
 ; they can be used anywhere, in section names, parameter names, or values.
-modules = consistency camb sigma8_rescale ${HALOFIT} load_nz shear_shear 2pt cfhtlens
+modules = consistency camb load_nz shear_shear 2pt cfhtlens
 values = demos/values10.ini
 
 likelihoods = cfhtlens
@@ -28,7 +28,7 @@ debug=F
 ;tells cosmosis how and where to save liklihood values
 [output]
 format=text
-filename=demo10_output_${HALOFIT}.txt
+filename=output/demo10_${HALOFIT}.txt
 
 [consistency]
 file = ./utility/consistency/consistency_interface.py
@@ -39,23 +39,17 @@ mode=all
 lmax=2500
 feedback=0
 
-; new simple module for using sigma8 as an input instead of A_s
-; the only parameter required is the file location
-[sigma8_rescale]
-file = ./utility/sample_sigma8/sigma8_rescale.py
 
-;the same halofit module as used in demo 6
-[halofit]
-file = ./boltzmann/halofit/halofit_module.so
+[camb]
+file = boltzmann/camb/camb_interface.py
+mode = cmb
+lmax = 2800          ;max ell to use for cmb calculation
+feedback=0         ;amount of output to print
+AccuracyBoost=1.1 ;CAMB accuracy boost parameter
+NonLinear = pk
+halofit_version = ${HALOFIT}
 
-; new module. works the same as halofit but uses the Takahashi non-linear correction
-; will be called when the variable HALOFIT = halofit_takahashi
-[halofit_takahashi]
-file = ./boltzmann/halofit_takahashi/halofit_interface.so
 
-; This is an example of a very simple module - it simple loads in 
-; a file once at the start when setting up and supplies fixed data
-; when executed.
 
 [load_nz]
 file = ./number_density/load_nz_fits/load_nz_fits.py
@@ -69,9 +63,9 @@ data_sets=SAMPLE
 
 [shear_shear]
 file = ./structure/projection/project_2d.py
-ell_min = 20.0
-ell_max = 10000.0
-n_ell = 150
+ell_min_logspaced = 20.0
+ell_max_logspaced = 10000.0
+n_ell_logspaced = 150
 shear-shear = sample-sample 
 verbose = F