diff -r -c -b -B cosmomc/params.ini cosmomc_sdss/params.ini *** cosmomc/params.ini 2004-06-19 21:12:18.000000000 +0200 --- cosmomc_sdss/params.ini 2004-06-19 22:26:02.000000000 +0200 *************** *** 32,37 **** --- 32,38 ---- use_CMB = T use_HST = T use_2dF = F + use_sdss = F use_clusters = F use_BBN = F use_Age_Tophat_Prior = T diff -r -c -b -B cosmomc/source/cmbtypes.f90 cosmomc_sdss/source/cmbtypes.f90 *** cosmomc/source/cmbtypes.f90 2004-06-19 21:12:18.000000000 +0200 --- cosmomc_sdss/source/cmbtypes.f90 2004-06-19 22:27:59.000000000 +0200 *************** *** 17,24 **** integer, parameter :: num_matter_power = 49 !number of points computed in matter power spectrum real, parameter :: matter_power_minkh = 3.16227766E-03 !minimum value of k/h to store real, parameter :: matter_power_dlnkh = 0.143911568 !log spacing in k/h ! real, parameter :: matter_power_maxz = 0. !6.0 ! integer, parameter :: matter_power_lnzsteps = 1 !20 !Note these are the interpolated/extrapolated values. The k at which matter power is computed up to --- 17,24 ---- integer, parameter :: num_matter_power = 49 !number of points computed in matter power spectrum real, parameter :: matter_power_minkh = 3.16227766E-03 !minimum value of k/h to store real, parameter :: matter_power_dlnkh = 0.143911568 !log spacing in k/h ! real, parameter :: matter_power_maxz = 0.17 !6.0 ! integer, parameter :: matter_power_lnzsteps = 2 !20 !Note these are the interpolated/extrapolated values. The k at which matter power is computed up to *************** *** 323,335 **** integer i x = log(kh/matter_power_minkh) / matter_power_dlnkh if (x < 0 .or. x >= num_matter_power-1) then ! write (*,*) ' k/h out of bounds in MatterPowerAt (',kh,')' ! stop ! end if i = int(x) d = x - i MatterPowerAt = exp(log(T%matter_power(i+1,1))*(1-d) + log(T%matter_power(i+2,1))*d) !Just do linear interpolation in logs for now.. !(since we already cublic-spline interpolated to get the stored values) --- 323,342 ---- integer i x = log(kh/matter_power_minkh) / matter_power_dlnkh + + ! SDSS window function need a wide range. Setting nul power spectrum + !from extremal values does not change the chi2 significantly. + ! if (dz...) to respect array bound + if (x < 0 .or. x >= num_matter_power-1) then ! ! write (*,*) ' k/h out of bounds in MatterPowerAt (',kh,')' ! ! stop ! MatterPowerAt=0. ! else i = int(x) d = x - i MatterPowerAt = exp(log(T%matter_power(i+1,1))*(1-d) + log(T%matter_power(i+2,1))*d) + endif !Just do linear interpolation in logs for now.. !(since we already cublic-spline interpolated to get the stored values) *************** *** 355,364 **** stop end if x = log(kh/matter_power_minkh) / matter_power_dlnkh if (x < 0 .or. x >= num_matter_power-1) then ! write (*,*) ' k/h out of bounds in MatterPowerAt_Z (',kh,')' ! stop ! end if iz = int(y) dz = y - iz --- 362,377 ---- stop end if x = log(kh/matter_power_minkh) / matter_power_dlnkh + + ! SDSS window function need a wide range. Setting nul power spectrum + !from extremal values does not change the chi2 significantly. + ! if (dz...) to respect array bound + if (x < 0 .or. x >= num_matter_power-1) then ! ! write (*,*) ' k/h out of bounds in MatterPowerAt_Z (',kh,')' ! ! stop ! MatterPowerAt_Z=0. ! else iz = int(y) dz = y - iz *************** *** 366,376 **** i = int(x) d = x - i ! mup = log(T%matter_power(i+1,iz+2))*(1-d) + log(T%matter_power(i+2,iz+2))*d ! mdn = log(T%matter_power(i+1,iz+1))*(1-d) + log(T%matter_power(i+2,iz+1))*d MatterPowerAt_Z = exp(mdn*(1-dz) + mup*dz) end function MatterPowerAt_Z --- 379,393 ---- i = int(x) d = x - i ! mup=0. ! mdn=0. ! if (dz.ne.0.) mup = log(T%matter_power(i+1,iz+2))*(1-d) + log(T%matter_power(i+2,iz+2))*d ! if (dz.ne.1.) mdn = log(T%matter_power(i+1,iz+1))*(1-d) + log(T%matter_power(i+2,iz+1))*d MatterPowerAt_Z = exp(mdn*(1-dz) + mup*dz) + endif + end function MatterPowerAt_Z diff -r -c -b -B cosmomc/source/driver.F90 cosmomc_sdss/source/driver.F90 *** cosmomc/source/driver.F90 2004-06-19 21:12:18.000000000 +0200 --- cosmomc_sdss/source/driver.F90 2004-06-19 22:28:52.000000000 +0200 *************** *** 133,139 **** Use_2dF = Ini_Read_Logical('use_2dF',.false.) Use_Clusters = Ini_Read_Logical('use_clusters',.false.) ! Use_HST = Ini_Read_Logical('use_HST',.true.) Use_BBN = Ini_Read_Logical('use_BBN',.false.) Use_Age_Tophat_Prior= Ini_Read_Logical('use_Age_Tophat_Prior',.true.) --- 133,139 ---- Use_2dF = Ini_Read_Logical('use_2dF',.false.) Use_Clusters = Ini_Read_Logical('use_clusters',.false.) ! Use_sdss = Ini_Read_Logical('use_sdss',.false.) Use_HST = Ini_Read_Logical('use_HST',.true.) Use_BBN = Ini_Read_Logical('use_BBN',.false.) Use_Age_Tophat_Prior= Ini_Read_Logical('use_Age_Tophat_Prior',.true.) *************** *** 142,148 **** Use_WeakLen = Ini_Read_Logical('use_WeakLen',.false.) Use_min_zre = Ini_Read_Double('use_min_zre',0.d0) ! use_LSS = Use_2dF .or. Use_Clusters .or. Use_WeakLen Temperature = Ini_Read_Real('temperature',1.) --- 142,148 ---- Use_WeakLen = Ini_Read_Logical('use_WeakLen',.false.) Use_min_zre = Ini_Read_Double('use_min_zre',0.d0) ! use_LSS = Use_2dF .or. Use_Clusters .or. Use_WeakLen .or. Use_sdss Temperature = Ini_Read_Real('temperature',1.) diff -r -c -b -B cosmomc/source/twodf.f90 cosmomc_sdss/source/twodf.f90 *** cosmomc/source/twodf.f90 2004-06-19 21:12:18.000000000 +0200 --- cosmomc_sdss/source/twodf.f90 2004-02-12 16:06:45.000000000 +0100 *************** *** 1,5 **** ! !Module for using 147k 2dF data !See http://www.roe.ac.uk/~wjp/2dFGRS/2dF_Pk.html - adapted for f90 module twodf use settings --- 1,9 ---- ! !Module for using 147k 2dF data and SDSS data !See http://www.roe.ac.uk/~wjp/2dFGRS/2dF_Pk.html - adapted for f90 + !See http://www.hep.upenn.edu/~max/sdss.html - adapted for f90 + + ! version february 2004. + ! TODO: introduce bias as fast parameter; implement corrections from HaloFit. module twodf use settings *************** *** 14,19 **** --- 18,54 ---- real :: twodf_P(twodf_points), twodf_k(twodf_points), twodf_Ninv(twodf_points,twodf_points) logical :: do_twodf_init = .true. + real, parameter :: twodf_redshift = 0.17 + + !SDSS declaration + + logical :: use_sdss = .false. + + !To use hyper parameters between twodf and sdss + logical :: use_lss_hyper=.false. + + integer, parameter :: sdss_bands_full = 22 + integer, parameter :: sdss_kbands_full = 97 + + !The window functions are defined up to kh=100 for kbands=97. The + !observed power spectrum is measured up to kh=0.3 for bands=22. Linear + !regime is trusted up to k/h=0.15 or bands=17 (k/h=0.2 for 19). The + !extrapolated power spectrum is arbitrarily set to zero for k/h>3 (see + !cmbtype.f90) + + integer, parameter :: sdss_bands = 17 + integer, parameter :: sdss_kbands = 97 + real, parameter :: sdss_redshift = 0.1 + real, dimension(sdss_bands) :: sdss_P + real, dimension(sdss_kbands) :: sdss_k + real, dimension(sdss_bands,sdss_kbands) :: sdss_W + logical :: do_sdss_init = .true. + + !Ln(Gamma(ndata/2 + 1)) for hyperparameters + real :: Twodfo2p1,Sdsso2p1 + real :: LnGammaTwodf,LnGammaSdss + + contains subroutine twodf_init *************** *** 28,36 **** end do close(tmp_file_unit) call ReadMatrix('data/twodf_june02_matrix.dat',twodf_Ninv,twodf_points,twodf_points) ! do_twodf_init = .false. ! read (tmp_file_unit,*, end = 300, err=300) model 200 stop 'Error reading 2dF file' --- 63,72 ---- end do close(tmp_file_unit) call ReadMatrix('data/twodf_june02_matrix.dat',twodf_Ninv,twodf_points,twodf_points) ! Twodfo2p1 = real(twodf_points)/2. + 1. ! LnGammaTwodf = gammln(Twodfo2p1) do_twodf_init = .false. ! if (Feedback > 2) write(*,*) 'Ln[Gamma(1+twodf_points/2)]= ',LnGammaTwodf read (tmp_file_unit,*, end = 300, err=300) model 200 stop 'Error reading 2dF file' *************** *** 76,82 **** do imu = 1, nint2 rmu = -1 + (imu-0.5)*dmu yk = sqrt(rk2 + xk2 - 2*xk*rk*rmu) ! sp = sp + MatterPowerAt(Theory, yk) end do s = s + sp*xk2*twodf_Wksq(xk) xk = rkmax; xk2= xk**2 --- 112,119 ---- do imu = 1, nint2 rmu = -1 + (imu-0.5)*dmu yk = sqrt(rk2 + xk2 - 2*xk*rk*rmu) ! ! sp = sp + MatterPowerAt(Theory, yk) ! sp = sp + MatterPowerAt_Z(Theory, yk,twodf_redshift) end do s = s + sp*xk2*twodf_Wksq(xk) xk = rkmax; xk2= xk**2 *************** *** 99,105 **** do imu = 1, nint2 rmu = -1 + (imu-0.5)*dmu yk = sqrt(rk2+xk2-2*xk*rk*rmu) ! sp = sp + MatterPowerAt(Theory, yk) end do s = s+ sp*xk2*twodf_Wksq(xk) end do --- 136,143 ---- do imu = 1, nint2 rmu = -1 + (imu-0.5)*dmu yk = sqrt(rk2+xk2-2*xk*rk*rmu) ! ! sp = sp + MatterPowerAt(Theory, yk) ! sp = sp + MatterPowerAt_Z(Theory, yk,twodf_redshift) end do s = s+ sp*xk2*twodf_Wksq(xk) end do *************** *** 115,121 **** function LSS_2dflike(CMB, Theory) Type (CMBParams) CMB Type (CosmoTheory) Theory ! real LSS_2dflike real th(twodf_points), sum1, chisq real Ninvth(twodf_points), Ninvd(twodf_points) integer i --- 153,159 ---- function LSS_2dflike(CMB, Theory) Type (CMBParams) CMB Type (CosmoTheory) Theory ! real, dimension(2) :: LSS_2dflike real th(twodf_points), sum1, chisq real Ninvth(twodf_points), Ninvd(twodf_points) integer i *************** *** 141,165 **** chisq = sum(twodf_P*Ninvd) - sum(th*Ninvd)**2/sum1 + log(sum1) ! LSS_2dflike = chisq/2 if (Feedback>1) write(*,*) '2df chi-sq: ', chisq end function LSS_2dflike ! function LSSLnLike(CMB, Theory) Type (CMBParams) CMB Type (CosmoTheory) Theory real LSSLnLike ! if (use_2dF) then ! LSSLnLike = LSS_2dflike(CMB,Theory) ! else ! LSSLnLike = 0 ! end if ! end function LSSLnLike end module twodf --- 179,336 ---- chisq = sum(twodf_P*Ninvd) - sum(th*Ninvd)**2/sum1 + log(sum1) ! LSS_2dflike(1) = chisq/2. ! LSS_2dflike(2) = chisq/2. - log(sum1)/2. if (Feedback>1) write(*,*) '2df chi-sq: ', chisq end function LSS_2dflike ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! The SDSS stuff ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! subroutine sdss_init ! ! integer i,j,iopb ! real keff,klo,khi,beff ! real, dimension(sdss_bands_full,sdss_kbands_full) :: sdss_Wfull ! real, dimension(sdss_kbands_full) :: sdss_kfull ! real, dimension(sdss_bands) :: sdss_sdev ! character(80) :: dummychar ! ! call ReadVector('data/sdss_kbands.txt',sdss_kfull,sdss_kbands_full) ! sdss_k(1:sdss_kbands)=sdss_kfull(1:sdss_kbands) ! ! if (Feedback > 0) then ! write(*,*) 'reading: sdss data' ! write(*,*) 'kbands Windows truncated at k/h= ',sdss_k(sdss_kbands) ! endif ! call OpenTxtFile('data/sdss_measurements.txt', tmp_file_unit) ! sdss_P=0. ! read (tmp_file_unit,*)dummychar ! read (tmp_file_unit,*)dummychar ! do i =1, sdss_bands ! read (tmp_file_unit,*, iostat=iopb) keff,klo,khi,sdss_P(i),sdss_sdev(i),beff ! end do ! close(tmp_file_unit) ! ! if (Feedback > 0) write(*,*) 'bands truncated at keff= ',keff ! ! call ReadMatrix('data/sdss_windows.txt',sdss_Wfull,sdss_bands_full,sdss_kbands_full) ! sdss_W(1:sdss_bands,1:sdss_kbands)=sdss_Wfull(1:sdss_bands,1:sdss_kbands) ! ! ! ! Prewhiten so that WLOG N=I: ! do i=1,sdss_bands ! sdss_P(i) = sdss_P(i)/sdss_sdev(i) ! do j=1,sdss_kbands ! sdss_W(i,j) = sdss_W(i,j)/sdss_sdev(i) ! end do ! end do ! ! Now chi2 = (obs-W p)^t (obs-W p) = |obs-Wp|^2. ! Sdsso2p1 = real(sdss_bands)/2. + 1. ! LnGammaSdss = gammln(Sdsso2p1) ! do_sdss_init = .false. ! if (Feedback > 2) write(*,*) 'Ln[Gamma(1+sdss_bands/2)]= ',LnGammaSdss ! if (iopb.ne.0) then ! stop 'Error reading sdss file' ! endif ! ! end subroutine sdss_init ! ! ! function LSS_sdsslike(CMB, Theory) ! Type (CMBParams) CMB ! Type (CosmoTheory) Theory ! real, dimension(sdss_kbands) :: sdss_Pth ! real, dimension(sdss_bands) :: sdss_WPth ! real :: chisq,normV ! real, dimension(2) :: LSS_sdsslike ! integer :: i ! ! if (do_sdss_init) call sdss_init ! ! do i=1, sdss_kbands ! sdss_Pth(i)=MatterPowerAt_Z(Theory,sdss_k(i),sdss_redshift) ! end do ! ! sdss_WPth = matmul(sdss_W,sdss_Pth) ! ! !with analytic marginalization over normalization nuisance (flat prior) ! normV = sum(sdss_WPth*sdss_WPth) ! chisq = sum(sdss_P*sdss_P) - sum(sdss_WPth*sdss_P)**2/normV + log(normV) ! ! !without analytic marginalization ! ! chisq = sum((sdss_P(:) - sdss_WPth(:))**2) ! ! LSS_sdsslike(1) = chisq/2. ! LSS_sdsslike(2) = chisq/2. - log(normV)/2. ! ! if (Feedback>1) write(*,*) 'sdss chi-sq:', chisq ! ! end function LSS_sdsslike ! ! ! function gammln(xx) ! ! from numrec ! REAL gammln,xx ! INTEGER j ! DOUBLE PRECISION ser,stp,tmp,x,y,cof(6) ! SAVE cof,stp ! DATA cof,stp/76.18009172947146d0,-86.50532032941677d0, & ! 24.01409824083091d0,-1.231739572450155d0,.1208650973866179d-2, & ! -.5395239384953d-5,2.5066282746310005d0/ ! x=xx ! y=x ! tmp=x+5.5d0 ! tmp=(x+0.5d0)*log(tmp)-tmp ! ser=1.000000000190015d0 ! do j=1,6 ! y=y+1.d0 ! ser=ser+cof(j)/y ! enddo ! gammln=tmp+log(stp*ser/x) ! return ! END function gammln ! ! ! function LSSLnLike(CMB, Theory) Type (CMBParams) CMB Type (CosmoTheory) Theory real LSSLnLike + real, dimension(2) :: twodflike,sdsslike ! LSSLnLike = 0. ! if (use_lss_hyper.and.use_2dF.and.use_sdss) then ! twodflike = LSS_2dflike(CMB,Theory) ! sdsslike = LSS_sdsslike(CMB,Theory) ! ! LSSLnLike = Twodfo2p1*log(1.+twodflike(2)) + twodflike(1)-twodflike(2) - LnGammaTwodf & ! + Sdsso2p1*log(1.+sdsslike(2)) + sdsslike(1)-sdsslike(2) - LnGammaSdss ! ! if (Feedback>1) then ! write(*,*) 'Effective hyperparam 2dF: ',2.*real(twodf_points)/twodflike(2) ! write(*,*) 'Effective hyperparam sdss: ',2.*real(sdss_bands)/sdsslike(2) ! write(*,*) 'LSS marginalized -ln(like):',LSSLnLike ! endif ! ! if (Feedback>2) then ! write(*,*) 'Ln[Evidence ratio for hyperparams]: ',-LSSLnLike + twodflike(1) + sdsslike(1) ! endif + else + if (use_2dF) then + twodflike = LSS_2dflike(CMB,Theory) + LSSLnLike = twodflike(1) + LSSLnLike + endif + if (use_sdss) then + sdsslike = LSS_sdsslike(CMB,Theory) + LSSLnLike = sdsslike(1) + LSSLnLike + endif + endif + end function LSSLnLike end module twodf