Adds footprint option to cluster likelihoods

.gitignore

@@ -102,3 +102,19 @@ venv.bak/
 
 # mypy
 .mypy_cache/
+
+
+# clusters
+soliket/clusters/data/selFn*
+soliket/clusters/data/*zip
+soliket/clusters/data/*fits
+soliket/clusters/chains
+soliket/clusters/notebooks/figures
+soliket/clusters/data/act/*
+soliket/clusters/data/advact/*
+.DS_Store
+soliket/ymap
+soliket/cosmopower
+soliket/clusters/checks/*
+soliket/clusters/notebooks/figs/*
+soliket/clusters/notebooks/ongoing/*

requirements.txt

@@ -16,3 +16,4 @@ pyccl
 sacc
 fgspectra>=1.1.0
 syslibrary
+nemo-sz>=0.6.1

setup.cfg

@@ -25,6 +25,7 @@ install_requires =
     sacc
     fgspectra>=1.1.0
     syslibrary
+    nemo-sz>=0.6.1
 
 [options.package_data]
 soliket = *.yaml,*.bibtex,clusters/data/*,clusters/data/selFn_equD56/*,lensing/data/*.txt,lensing/data/*.fits,mflike/*.yaml,tests/*.yaml,data/xcorr_simulated/*.txt,data/CosmoPower/CP_paper/CMB/*.pkl,tests/data/test_bandpass/*,cross_correlation/*.yaml

soliket/__init__.py

@@ -2,6 +2,7 @@
 from .gaussian import GaussianLikelihood, MultiGaussianLikelihood  # noqa: F401
 # from .studentst import StudentstLikelihood  # noqa: F401
 from .ps import PSLikelihood, BinnedPSLikelihood  # noqa: F401
+from .clusters import BinnedClusterLikelihood, UnbinnedClusterLikelihood  # noqa: F401
 from .mflike import MFLike  # noqa: F401
 from .mflike import TheoryForge_MFLike
 from .cross_correlation import GalaxyKappaLikelihood, ShearKappaLikelihood  # noqa: F401, E501
@@ -12,7 +13,7 @@
 from .ccl import CCL  # noqa: F401
 
 try:
-    from .clusters import ClusterLikelihood  # noqa: F401
+    from .clusters import BinnedClusterLikelihood, UnbinnedClusterLikelihood  # noqa: F401
 except ImportError:
     print('Skipping cluster likelihood (is pyCCL installed?)')
-    pass
+    pass

soliket/cash.py

@@ -1,31 +1,26 @@
-import numpy as np
-from typing import Optional
-
 from cobaya.likelihood import Likelihood
 from .cash_data import CashCData
 
 
 class CashCLikelihood(Likelihood):
     name: str = "Cash-C"
-    datapath = Optional[str]
 
     def initialize(self):
-
-        x, N = self._get_data()
+        N = self._get_data()
         self.data = CashCData(self.name, N)
 
     def _get_data(self):
-        data = np.loadtxt(self.datapath, unpack=False)
-        N = data[:, -1] # assume data stored like column_stack([z, q, N])
-        x = data[:, :-1]
-        return x, N
+        raise NotImplementedError
 
     def _get_theory(self, **kwargs):
         if ("cash_test_logp" in kwargs):
             return np.arange(kwargs["cash_test_logp"])
         else:
             raise NotImplementedError
 
-    def logp(self, **params_values):
-        theory = self._get_theory(**params_values)
+    def logp(self, **kwargs):
+        pk_intp = self.provider.get_Pk_interpolator()
+        theory = self._get_theory(pk_intp, **kwargs)
         return self.data.loglike(theory)
+
+

soliket/cash_data.py

@@ -1,39 +1,62 @@
 import numpy as np
 from scipy.special import factorial
-import math as m
 
 
-def cash_c_logpdf(theory, data, usestirling=True):
+def cash_c_logpdf(theory, data, usestirling=True, name="binned"):
 
-    data = np.asarray(data, dtype=int)
+    __, __, delN2Dcat, zcut = data
 
-    ln_fac = np.zeros_like(data, dtype=float)
+    obs = np.asarray(delN2Dcat, dtype=int)
+    ln_fac = np.zeros_like(obs, dtype=float)
+    zcut_arr = np.arange(zcut)
+
+    # print("WARNING: Zapped some S/N bins as a test")
+    # theory=theory[:, 1:]
+    # obs=obs[:, 1:]
+    # ln_fac=ln_fac[:, 1:]
+
+    if zcut > 0:
+        theory = np.delete(theory, zcut_arr, 0)
+        obs = np.delete(obs, zcut_arr, 0)
+        ln_fac = np.delete(ln_fac, zcut_arr, 0)
+        print("\r ::: Excluding first {} redshift bins in likelihood.".format(zcut))
+
+        for i in range(theory.shape[0]):
+            print('\r Number of clusters in redshift bin {}: {}.'.format(i, theory[i,:].sum()))
+        print('------------')
+        for i in range(theory.shape[1]):
+            print('\r Number of clusters in SNR bin {}: {}.'.format(i, theory[:,i].sum()))
+        print('------------')
+        print('\r Total predicted N = {}'.format(theory.sum()))
+        print('\r Total observed N = {}'.format(obs.sum()))
 
     if usestirling: # use Stirling's approximation for N > 10
-        ln_fac[data > 10] = 0.918939 + (data[data > 10] + 0.5) \
-                                    * np.log(data[data > 10]) - data[data > 10]
-        ln_fac[data <= 10] = np.log(factorial(data[data <= 10]))
-    else:
-        ln_fac[data > 0] = np.log(factorial(data[data > 0]))
-    ln_fac[data == 0] = 0.
+        ln_fac[obs > 10] = 0.918939 + (obs[obs > 10] + 0.5) * np.log(obs[obs > 10]) - obs[obs > 10]
+        ln_fac[obs <= 10] = np.log(factorial(obs[obs <= 10]))
+    else: # direct compuation of factorial
+        ln_fac[obs > 0] = np.log(factorial(obs[obs > 0]))
+    ln_fac[obs == 0] = 0.
+
+    log_theory = np.zeros_like(theory, dtype=float)
+    log_theory[theory > 0] = np.log(theory[theory > 0])
+    log_theory[theory == 0] = 0.
 
-    loglike = data * np.log(theory) - theory - ln_fac
+    loglike = obs * log_theory - theory - ln_fac
+
+    print("\r ::: 2D ln likelihood = ", np.nansum(loglike[np.isfinite(loglike)]))
 
     return np.nansum(loglike[np.isfinite(loglike)])
 
 
+
 class CashCData:
     """Named multi-dimensional Cash-C distributed data
     """
 
     def __init__(self, name, N, usestirling=True):
-
         self.name = str(name)
         self.data = N
         self.usestirling = usestirling
 
-    def __len__(self):
-        return len(self.data)
-
     def loglike(self, theory):
-        return cash_c_logpdf(theory, self.data)
+        return cash_c_logpdf(theory, self.data, name=self.name)

soliket/clusters/__init__.py

@@ -1 +1,2 @@
-from .clusters import ClusterLikelihood  # noqa: F401
+from .clusters import BinnedClusterLikelihood, UnbinnedClusterLikelihood  # noqa: F401
+from .ccl_th import CCL