Initial commit of xcorr (#22)

* initial commit of xcorr * checks against CCL, fixed units in calculation of mag kernel * added sacc file handling * exposed defaults in xcorr initialize * removed dummy_pk * updated code style * added proposal for xcorr b parm * added docstring, made sacc load default behaviour * whitespace * added other parms to docstring, smaller b1 proposal * added s1 proposal * removed redundant chi setup option * changed to use soliket constants * changed theory to provider * rename cosmo to provider for clarity * changed chi result to dict * removed zeff and autoCMB for now * tidy up codestyle Co-authored-by: Ian Harrison <[email protected]>
simonsobs · Oct 25, 2021 · 1d0a7fc · 1d0a7fc
1 parent 17a02da
commit 1d0a7fc
Show file tree

Hide file tree

Showing 14 changed files with 2,463 additions and 2 deletions.
diff --git a/notebooks/dev/cross_correlation_simulatedata.ipynb b/notebooks/dev/cross_correlation_simulatedata.ipynb
diff --git a/soliket/__init__.py b/soliket/__init__.py
@@ -3,6 +3,8 @@
 from .ps import PSLikelihood, BinnedPSLikelihood  # noqa: F401
 from .clusters import ClusterLikelihood  # noqa: F401
 from .mflike import MFLike  # noqa: F401
+from .xcorr import XcorrLikelihood  # noqa: F401
+
 try:
     import pyccl as ccl  # noqa: F401
     from .ccl import CCL  # noqa: F401

diff --git a/soliket/constants.py b/soliket/constants.py
@@ -2,6 +2,7 @@
 
 C_M_S = constants.c
 C_KM_S = constants.c * 1.e-3
+C_HMPC = constants.c * 1.e-5
 h_Planck = constants.h
 k_Boltzmann = constants.k
 elementary_charge = constants.e

diff --git a/soliket/tests/data/unwise_g-so_kappa.sim.sacc.fits b/soliket/tests/data/unwise_g-so_kappa.sim.sacc.fits
diff --git a/soliket/tests/test_runs.py b/soliket/tests/test_runs.py
@@ -10,7 +10,8 @@
                           "lensing",
                           "lensing_lite",
                           "multi",
-                          "cross_correlation"])
+                          "cross_correlation",
+                          "xcorr"])
 def test_evaluate(lhood):
     info = yaml_load(pkgutil.get_data("soliket", f"tests/test_{lhood}.yaml"))
     info["force"] = True
@@ -24,7 +25,8 @@ def test_evaluate(lhood):
                           "lensing",
                           "lensing_lite",
                           "multi",
-                          "cross_correlation"])
+                          "cross_correlation",
+                          "xcorr"])
 def test_mcmc(lhood):
     info = yaml_load(pkgutil.get_data("soliket", f"tests/test_{lhood}.yaml"))
     info["force"] = True

diff --git a/soliket/tests/test_xcorr.py b/soliket/tests/test_xcorr.py
@@ -0,0 +1,171 @@
+# pytest -k xcorr -v --pdb .
+
+import pytest
+import numpy as np
+
+from cobaya.yaml import yaml_load
+from cobaya.model import get_model
+
+import os
+import pdb
+
+
+def get_demo_xcorr_model(theory):
+    if theory == "camb":
+        info_yaml = r"""
+        likelihood:
+            soliket.XcorrLikelihood:
+                stop_at_error: True
+                datapath: soliket/tests/data/unwise_g-so_kappa.sim.sacc.fits
+                k_tracer_name: ck_so
+                gc_tracer_name: gc_unwise
+
+        theory:
+            camb:
+                extra_args:
+                  lens_potential_accuracy: 1
+
+        params:
+            tau: 0.05
+            mnu: 0.0
+            nnu: 3.046
+            b1:
+                prior:
+                    min: 0.
+                    max: 10.
+                ref:
+                    min: 1.
+                    max: 4.
+                proposal: 0.1
+            s1:
+                prior:
+                    min: 0.1
+                    max: 1.0
+                proposal: 0.1
+        """
+    elif theory == "classy":
+        info_yaml = r"""
+        likelihood:
+            soliket.XcorrLikelihood:
+                stop_at_error: True
+                datapath: soliket/tests/data/unwise_g-so_kappa.sim.sacc.fits
+                k_tracer_name: ck_so
+                gc_tracer_name: gc_unwise
+
+        theory:
+            classy:
+                extra_args:
+                    output: lCl, tCl
+                path: global
+
+        params:
+            b1:
+                prior:
+                    min: 0.
+                    max: 10.
+                ref:
+                    min: 1.
+                    max: 4.
+                proposal: 0.1
+            s1:
+                prior:
+                    min: 0.1
+                    max: 1.0
+                proposal: 0.1
+
+        """
+
+    info = yaml_load(info_yaml)
+    model = get_model(info)
+    return model
+
+
+@pytest.mark.parametrize("theory", ["camb"])#, "classy"])
+def test_xcorr(theory):
+
+    params = {'b1': 1.0, 's1': 0.4}
+
+    model = get_demo_xcorr_model(theory)
+
+    lnl = model.loglike(params)[0]
+    assert np.isfinite(lnl)
+
+    xcorr_lhood = model.likelihood['soliket.XcorrLikelihood']
+
+    setup_chi_out = xcorr_lhood._setup_chi()
+
+    Pk_interpolator = xcorr_lhood.theory.get_Pk_interpolator(("delta_nonu", "delta_nonu"),
+                                                             extrap_kmax=1.e8,
+                                                             nonlinear=False).P
+
+    from soliket.xcorr.limber import do_limber
+
+    cl_gg, cl_kappag = do_limber(xcorr_lhood.ell_range,
+                                   xcorr_lhood.provider,
+                                   xcorr_lhood.dndz,
+                                   xcorr_lhood.dndz,
+                                   params['s1'],
+                                   params['s1'],
+                                   Pk_interpolator,
+                                   params['b1'],
+                                   params['b1'],
+                                   xcorr_lhood.alpha_auto,
+                                   xcorr_lhood.alpha_cross,
+                                   setup_chi_out,
+                                   Nchi=xcorr_lhood.Nchi,
+                                   dndz1_mag=xcorr_lhood.dndz,
+                                   dndz2_mag=xcorr_lhood.dndz)
+
+    ell_load = xcorr_lhood.data.x
+    cl_load = xcorr_lhood.data.y
+    # cov_load = xcorr_lhood.data.cov
+    # cl_err_load = np.sqrt(np.diag(cov_load))
+    n_ell = len(ell_load) // 2
+
+    ell_obs_gg = ell_load[n_ell:]
+    ell_obs_kappag = ell_load[:n_ell]
+
+    cl_obs_gg = cl_load[:n_ell]
+    cl_obs_kappag = cl_load[n_ell:]
+
+    # Nell_unwise_g = np.ones_like(cl_gg) \
+    #                         / (xcorr_lhood.ngal * (60 * 180 / np.pi)**2)
+    Nell_obs_unwise_g = np.ones_like(cl_obs_gg) \
+                            / (xcorr_lhood.ngal * (60 * 180 / np.pi)**2)
+
+    import pyccl as ccl
+    h2 = (xcorr_lhood.provider.get_param('H0') / 100)**2
+
+    cosmo = ccl.Cosmology(Omega_c=xcorr_lhood.provider.get_param('omch2') / h2,
+                          Omega_b=xcorr_lhood.provider.get_param('ombh2') / h2,
+                          h=xcorr_lhood.provider.get_param('H0') / 100,
+                          n_s=xcorr_lhood.provider.get_param('ns'),
+                          A_s=xcorr_lhood.provider.get_param('As'),
+                          Omega_k=xcorr_lhood.provider.get_param('omk'),
+                          Neff=xcorr_lhood.provider.get_param('nnu'),
+                          matter_power_spectrum='linear')
+
+    g_bias_zbz = (xcorr_lhood.dndz[:, 0],
+                  params['b1'] * np.ones(len(xcorr_lhood.dndz[:, 0])))
+    mag_bias_zbz = (xcorr_lhood.dndz[:, 0],
+                    params['s1'] * np.ones(len(xcorr_lhood.dndz[:, 0])))
+
+    tracer_g = ccl.NumberCountsTracer(cosmo,
+                                      has_rsd=False,
+                                      dndz=xcorr_lhood.dndz.T,
+                                      bias=g_bias_zbz,
+                                      mag_bias=mag_bias_zbz)
+
+    tracer_k = ccl.CMBLensingTracer(cosmo, z_source=1100)
+
+    cl_gg_ccl = ccl.cls.angular_cl(cosmo, tracer_g, tracer_g, xcorr_lhood.ell_range)
+    cl_kappag_ccl = ccl.cls.angular_cl(cosmo, tracer_k, tracer_g, xcorr_lhood.ell_range)
+
+    assert np.allclose(cl_gg_ccl, cl_gg)
+    assert np.allclose(cl_kappag_ccl, cl_kappag)
+
+    cl_obs_gg_ccl = ccl.cls.angular_cl(cosmo, tracer_g, tracer_g, ell_obs_gg)
+    cl_obs_kappag_ccl = ccl.cls.angular_cl(cosmo, tracer_k, tracer_g, ell_obs_kappag)
+
+    assert np.allclose(cl_obs_gg_ccl + Nell_obs_unwise_g, cl_obs_gg)
+    assert np.allclose(cl_obs_kappag_ccl, cl_obs_kappag)
diff --git a/soliket/tests/test_xcorr.yaml b/soliket/tests/test_xcorr.yaml
@@ -0,0 +1,33 @@
+# A simple cobaya likelihood for SO
+
+debug: True
+
+likelihood:
+  soliket.XcorrLikelihood:
+      stop_at_error: True
+
+params:
+  b1: 
+    prior:
+      min: 0.
+      max: 10.
+    ref:
+      min: 1.
+      max: 4.
+    proposal: 0.1
+    latex: b_1
+  s1: 
+    value: 0.4
+    latex: s_1
+
+theory:
+  camb:
+    stop_at_error: False
+    extra_args:
+      lens_potential_accuracy: 1
+
+sampler:
+  evaluate:
+
+
+output: chains/test_xcorr
diff --git a/soliket/xcorr/XcorrLikelihood.yaml b/soliket/xcorr/XcorrLikelihood.yaml
@@ -0,0 +1,41 @@
+auto_file:
+cross_file:
+dndz_file:
+
+datapath: soliket/tests/data/unwise_g-so_kappa.sim.sacc.fits
+k_tracer_name: ck_so
+gc_tracer_name: gc_unwise
+# tracer_combinations: gc_unwise-ck_so, gc_unwise-gc_unwise
+
+high_ell: 600
+nz: 149
+Nchi: 149
+Nchi_mag: 149
+
+Pk_interp_kmax: 10.0
+
+params:
+  b1: 
+    prior:
+      min: 0.
+      max: 10.
+    ref:
+      min: 1.
+      max: 4.
+    proposal: 0.1
+  s1: 
+    value: 0.4
+  H0:
+    value: 70.
+  omegam:
+    value: 0.3
+  ombh2:
+    value: 0.0245
+  omch2:
+    value: 0.1225
+  omk:
+    value: 0.0
+  ns:
+    value: 0.965
+  As:
+    value: 2.11e-9
diff --git a/soliket/xcorr/__init__.py b/soliket/xcorr/__init__.py
@@ -0,0 +1 @@
+from .xcorr import XcorrLikelihood
diff --git a/soliket/xcorr/data/clgg.txt b/soliket/xcorr/data/clgg.txt
@@ -0,0 +1,3 @@
+3.550000000000000000e+01 7.650000000000000000e+01 1.265000000000000000e+02 1.765000000000000000e+02 2.265000000000000000e+02 2.765000000000000000e+02 3.265000000000000000e+02 3.765000000000000000e+02 4.265000000000000000e+02 4.765000000000000000e+02 5.265000000000000000e+02 5.765000000000000000e+02
+1.171259386682930336e-07 1.074925971098873865e-07 7.605910992260147274e-08 5.561247886809266997e-08 4.345656443980961222e-08 3.471274000295990805e-08 2.776643347807964378e-08 2.248227907027175782e-08 1.867015082977402487e-08 1.590664246359614681e-08 1.377366655041607371e-08 1.201950138681249373e-08
+1.171386555992278504e-08 1.074818774537805139e-08 7.605695767423656527e-09 5.562085700838220625e-09 4.345907898321461207e-09 3.470700836078987746e-09 2.777317313635696080e-09 2.248324340497516305e-09 1.866778760685389341e-09 1.590802127607215577e-09 1.377460168904915491e-09 1.201961522334897607e-09
diff --git a/soliket/xcorr/data/clkg.txt b/soliket/xcorr/data/clkg.txt
@@ -0,0 +1,3 @@
+3.550000000000000000e+01 7.650000000000000000e+01 1.265000000000000000e+02 1.765000000000000000e+02 2.265000000000000000e+02 2.765000000000000000e+02 3.265000000000000000e+02 3.765000000000000000e+02 4.265000000000000000e+02 4.765000000000000000e+02 5.265000000000000000e+02 5.765000000000000000e+02
+1.425194692235914338e-07 1.234653760415391458e-07 8.611126066544479263e-08 6.248167114060786722e-08 4.838074702560995877e-08 3.853147673160977486e-08 3.090768970482272427e-08 2.506062365326854316e-08 2.078165440429319254e-08 1.768399884093924078e-08 1.532187351346465496e-08 1.340877933987993602e-08
+1.424982171873059966e-08 1.234732647730792068e-08 8.611508401982062142e-09 6.248438515784782318e-09 4.837007692483073648e-09 3.852305256406458467e-09 3.090458671774709973e-09 2.505965588344152520e-09 2.078012214338600843e-09 1.768136302352318975e-09 1.532330752048001851e-09 1.340720288202413845e-09