diff --git a/soliket/__init__.py b/soliket/__init__.py index f9251ffa..44957e61 100644 --- a/soliket/__init__.py +++ b/soliket/__init__.py @@ -18,4 +18,5 @@ from .gaussian import GaussianLikelihood, MultiGaussianLikelihood from .lensing import LensingLikelihood, LensingLiteLikelihood from .ps import BinnedPSLikelihood, PSLikelihood +from .szlike import KSZLikelihood, TSZLikelihood from .xcorr import XcorrLikelihood diff --git a/soliket/constants.py b/soliket/constants.py index 45dd45c0..b8ed461b 100644 --- a/soliket/constants.py +++ b/soliket/constants.py @@ -7,10 +7,13 @@ k_Boltzmann = constants.k elementary_charge = constants.e electron_mass_kg = constants.m_e +proton_mass_kg = constants.m_p T_CMB = 2.72548 # in K MSUN_CGS = 1.98840987e33 G_CGS = constants.G * 1e3 MPC2CM = constants.parsec * 1e8 +ST_CGS = constants.value('Thomson cross section') * 1.e4 +hydrogen_fraction = 0.76 full_sky_area_sqdeg = 41252.9612 SO_sky_area_sqdeg = 987.5 diff --git a/soliket/szlike/KSZLikelihood.yaml b/soliket/szlike/KSZLikelihood.yaml new file mode 100644 index 00000000..87404b29 --- /dev/null +++ b/soliket/szlike/KSZLikelihood.yaml @@ -0,0 +1,56 @@ +sz_data_file: gnfw_test_projections.txt +beam_file: beam_f150_daynight.txt +transform_type: Hankel #FFT or Hankel +#optional beam response file +beam_response: False #False, or give file +cov_ksz_file: cov_diskring_ksz_varweight_bootstrap.txt +#optional 2-halo file +twohalo_term: twohalo_cmass_average.txt #False, or give file. Need to flesh out this option more, but the idea is that if False it will use two_halo.py to calculate a 2h term + +redshift: 0.55 +mass_halo_mean_Msol: 1.e13 +frequency_GHz: 150. +#input_model: "gnfw" + +params: + Omega_m: + value: 0.25 + Omega_b: + value: 0.044 + hh: + value: 0.7 + Omega_L: + value: 0.75 + rhoc_0: + value: 2.77525e2 + C_OVER_HUBBLE: + value: 2997.9 + XH: + value: 0.76 #hydrogen fraction + v_rms: + value: 1.06e-3 #v_rms/c + + gnfw_rho0: + prior: + min: 1.0 + max: 6.0 + ref: 3.1864410705359125 + latex: \log_{10}\rho_0 + gnfw_bt_ksz: + prior: + min: 1.0 + max: 10.0 + ref: 3.45493977635 + latex: \beta_k + gnfw_A2h_ksz: + prior: + min: 0.01 + max: 5.0 + ref: 1.0 + latex: A_{k2h} + + +#prior: #al 1.0: + ans = ( + Dh + / np.sqrt(O_tot - 1.0) + * np.sinh(np.sqrt(O_tot - 1.0) * quad(ComInt, 0, z, args=provider)[0]) + ) + return ans + + +def AngDist(z, provider): + """angular distance""" + ans = ComDist(z, provider) / (1.0 + z) + return ans diff --git a/soliket/szlike/cov_diskring_ksz_varweight_bootstrap.txt b/soliket/szlike/cov_diskring_ksz_varweight_bootstrap.txt new file mode 100644 index 00000000..fcb0dae2 --- /dev/null +++ b/soliket/szlike/cov_diskring_ksz_varweight_bootstrap.txt @@ -0,0 +1,9 @@ +3.114388966750203853e-17 1.030466021667764340e-17 1.989903589431795886e-17 3.990747139864290666e-17 6.004681351830791974e-17 9.119641331246055472e-17 1.271410703778790521e-16 1.683135917755793293e-16 2.073662546189791365e-16 +1.030466021667764340e-17 1.830150670513891649e-16 1.849107550524976019e-16 2.994324746968246459e-16 4.993569819628552851e-16 7.874829963378285997e-16 1.120465128842835845e-15 1.510519395238511552e-15 1.883330966121615058e-15 +1.989903589431795886e-17 1.849107550524976019e-16 6.638526085148106398e-16 1.010364944329085446e-15 1.557366682051142589e-15 2.494360208378304450e-15 3.533796615818870131e-15 4.785553145152298330e-15 5.985057482699239769e-15 +3.990747139864290666e-17 2.994324746968246459e-16 1.010364944329085446e-15 2.647182926194105083e-15 4.204831771129999671e-15 6.400698092240780311e-15 9.184964656162738387e-15 1.257311137225875672e-14 1.591892840180708825e-14 +6.004681351830791974e-17 4.993569819628552851e-16 1.557366682051142589e-15 4.204831771129999671e-15 8.058461500115115439e-15 1.278421375477372909e-14 1.808888093954592656e-14 2.489170429910855343e-14 3.189271851269197359e-14 +9.119641331246055472e-17 7.874829963378285997e-16 2.494360208378304450e-15 6.400698092240780311e-15 1.278421375477372909e-14 2.236055684952939516e-14 3.256821147544840557e-14 4.512584635638064550e-14 5.810045151907715322e-14 +1.271410703778790521e-16 1.120465128842835845e-15 3.533796615818870131e-15 9.184964656162738387e-15 1.808888093954592656e-14 3.256821147544840557e-14 4.993574276989252032e-14 7.080973726244963718e-14 9.233652831754360816e-14 +1.683135917755793293e-16 1.510519395238511552e-15 4.785553145152298330e-15 1.257311137225875672e-14 2.489170429910855343e-14 4.512584635638064550e-14 7.080973726244963718e-14 1.043057404038440461e-13 1.386543053592368014e-13 +2.073662546189791365e-16 1.883330966121615058e-15 5.985057482699239769e-15 1.591892840180708825e-14 3.189271851269197359e-14 5.810045151907715322e-14 9.233652831754360816e-14 1.386543053592368014e-13 1.892093946543024183e-13 diff --git a/soliket/szlike/cov_diskring_tsz_varweight_bootstrap.txt b/soliket/szlike/cov_diskring_tsz_varweight_bootstrap.txt new file mode 100644 index 00000000..09d705f2 --- /dev/null +++ b/soliket/szlike/cov_diskring_tsz_varweight_bootstrap.txt @@ -0,0 +1,9 @@ +1.412192710030839058e-17 5.036712064238222188e-18 1.011790444493683631e-17 2.006345091227972931e-17 3.102982687531567003e-17 4.853033233404774898e-17 6.871439670268841824e-17 9.083432492648790024e-17 1.137987339383397044e-16 +5.036712064238222188e-18 8.280234483335549423e-17 8.159173019447813252e-17 1.333888523460657689e-16 2.319605378036357951e-16 3.646456945963971215e-16 5.107305946529757585e-16 6.788065465095132587e-16 8.457304964526317898e-16 +1.011790444493683631e-17 8.159173019447813252e-17 3.091414429180411255e-16 4.818530988493671526e-16 7.610721920794591973e-16 1.212627962072382284e-15 1.717740090712217118e-15 2.331312408656376973e-15 2.944011971508447919e-15 +2.006345091227972931e-17 1.333888523460657689e-16 4.818530988493671526e-16 1.266489127949988896e-15 2.039896793352412389e-15 3.110133008105503366e-15 4.457807519531365550e-15 6.118055926071910318e-15 7.778094826101351645e-15 +3.102982687531567003e-17 2.319605378036357951e-16 7.610721920794591973e-16 2.039896793352412389e-15 3.952118167130136203e-15 6.284179981077994492e-15 8.901258432542199811e-15 1.227987444500612835e-14 1.577478824267336246e-14 +4.853033233404774898e-17 3.646456945963971215e-16 1.212627962072382284e-15 3.110133008105503366e-15 6.284179981077994492e-15 1.096969697544288357e-14 1.598379082541362185e-14 2.218305091936667372e-14 2.863610761587509863e-14 +6.871439670268841824e-17 5.107305946529757585e-16 1.717740090712217118e-15 4.457807519531365550e-15 8.901258432542199811e-15 1.598379082541362185e-14 2.447308891192420811e-14 3.476542149692437069e-14 4.545693473121323415e-14 +9.083432492648790024e-17 6.788065465095132587e-16 2.331312408656376973e-15 6.118055926071910318e-15 1.227987444500612835e-14 2.218305091936667372e-14 3.476542149692437069e-14 5.129687488210757834e-14 6.835878336187716378e-14 +1.137987339383397044e-16 8.457304964526317898e-16 2.944011971508447919e-15 7.778094826101351645e-15 1.577478824267336246e-14 2.863610761587509863e-14 4.545693473121323415e-14 6.835878336187716378e-14 9.346118486458904190e-14 diff --git a/soliket/szlike/flat_map.py b/soliket/szlike/flat_map.py new file mode 100644 index 00000000..feea2bc3 --- /dev/null +++ b/soliket/szlike/flat_map.py @@ -0,0 +1,4599 @@ +import numpy as np +import matplotlib.pyplot as plt +import scipy +from scipy import stats, special, integrate +import sys +import os +from astropy.io import fits +from scipy.interpolate import interp1d, UnivariateSpline, RectBivariateSpline + +################################################################### + + +class FlatMap(object): + def __init__( + self, + nX=256, + nY=256, + sizeX=5.0 * np.pi / 180.0, + sizeY=5.0 * np.pi / 180.0, + name="test", + ): + """n is number of pixels on one side + size is the angular size of the side, in radians + """ + self.name = name + self.nX = nX + self.sizeX = sizeX + self.dX = float(sizeX) / (nX - 1) + x = self.dX * np.arange(nX) # the x value corresponds to the center of the cell + # + self.nY = nY + self.sizeY = sizeY + self.dY = float(sizeY) / (nY - 1) + y = self.dY * np.arange(nY) # the y value corresponds to the center of the cell + # + self.x, self.y = np.meshgrid(x, y, indexing="ij") + # + self.fSky = self.sizeX * self.sizeY / (4.0 * np.pi) + + self.data = np.zeros((nX, nY)) + + lx = np.zeros(nX) + lx[: int(nX / 2 + 1)] = 2.0 * np.pi / sizeX * np.arange(nX // 2 + 1) + lx[int(nX / 2 + 1):] = 2.0 * np.pi / sizeX * np.arange(-nX // 2 + 1, 0, 1) + ly = 2.0 * np.pi / sizeY * np.arange(nY // 2 + 1) + self.lx, self.ly = np.meshgrid(lx, ly, indexing="ij") + + self.ell = np.sqrt(self.lx**2 + self.ly**2) + self.dataFourier = np.zeros((nX, nY // 2 + 1)) + + def copy(self): + newMap = FlatMap( + nX=self.nX, nY=self.nY, sizeX=self.sizeX, sizeY=self.sizeY, name=self.name + ) + newMap.data = self.data.copy() + newMap.dataFourier = self.dataFourier.copy() + return newMap + + def write(self, path=None): + # primary hdu: size of map + prihdr = fits.Header() + prihdr["name"] = self.name + prihdr["nX"] = self.nX + prihdr["sizeX"] = self.sizeX + prihdr["nY"] = self.nY + prihdr["sizeY"] = self.sizeY + hdu0 = fits.PrimaryHDU(header=prihdr) + # secondary hdu: real-space and Fourier maps + c1 = fits.Column(name="data", format="D", array=self.data.flatten()) + c2 = fits.Column( + name="dataFourier", format="M", array=self.dataFourier.flatten() + ) + hdu1 = fits.BinTableHDU.from_columns([c1, c2]) + # + hdulist = fits.HDUList([hdu0, hdu1]) + + if path is None: + path = "./output/lens_simulator/" + self.name + ".fits" + # print("writing to "+path) + hdulist.writeto(path, overwrite=True) + + def read(self, path=None): + if path is None: + path = "./output/lens_simulator/" + self.name + ".fits" + # print("reading from "+path) + hdulist = fits.open(path) + # + # self.name = hdulist[0].header['name'] + self.nX = hdulist[0].header["nX"] + self.sizeX = hdulist[0].header["sizeX"] + self.nY = hdulist[0].header["nY"] + self.sizeY = hdulist[0].header["sizeY"] + # + self.__init__( + nX=self.nX, nY=self.nY, sizeX=self.sizeX, sizeY=self.sizeY, name=self.name + ) + # + self.data = hdulist[1].data["data"].reshape((self.nX, self.nY)) + self.dataFourier = ( + hdulist[1] + .data["dataFourier"][: self.nX * (self.nY // 2 + 1)] + .reshape((self.nX, self.nY // 2 + 1)) + ) + # + hdulist.close() + + def saveDataFourier(self, dataFourier, path): + # print("saving Fourier map to", path) + # primary hdu: size of map + prihdr = fits.Header() + prihdr["nX"] = self.nX + prihdr["sizeX"] = self.sizeX + prihdr["nY"] = self.nY + prihdr["sizeY"] = self.sizeY + hdu0 = fits.PrimaryHDU(header=prihdr) + # secondary hdu: maps + c1 = fits.Column(name="dataFourier", format="M", array=dataFourier.flatten()) + hdu1 = fits.BinTableHDU.from_columns([c1]) + # + hdulist = fits.HDUList([hdu0, hdu1]) + # + hdulist.writeto(path, overwrite=True) + + def loadDataFourier(self, path): + # print("reading fourier map from", path) + hdulist = fits.open(path) + dataFourier = ( + hdulist[1] + .data["dataFourier"][: self.nX * (self.nY // 2 + 1)] + .reshape((self.nX, self.nY // 2 + 1)) + ) + hdulist.close() + return dataFourier + + ############################################################################### + # change resolution of map + + def downResolution(self, nXNew, nYNew, data=None, test=False): + """Expects nXNew < self.nX and nYNew < self.nY""" + if data is None: + data = self.data.copy() + # Fourier transform the map + dataFourier = self.fourier(data) + # truncate the fourier map + newMap = FlatMap(nXNew, nYNew, sizeX=self.sizeX, sizeY=self.sizeY) + IX = range(nXNew // 2 + 1) + range(-nXNew // 2 + 1, 0) + IY = range(nYNew // 2 + 1) + IX, IY = np.meshgrid(IX, IY, indexing="ij") + newMap.dataFourier = dataFourier[IX, IY] + # update real space map + newMap.data = newMap.inverseFourier() + if test: + self.plot(data) + newMap.plot(newMap.data) + return newMap.data + + def upResolution(self, nXNew, nYNew, data=None, test=False): + """Expects nXNew > self.nX and nYNew > self.nY""" + if data is None: + data = self.data.copy() + # Fourier transform the map + dataFourier = self.fourier(data) + # zero-pad the fourier map + newMap = FlatMap(nXNew, nYNew, sizeX=self.sizeX, sizeY=self.sizeY) + IX = range(self.nX // 2 + 1) + range(-self.nX // 2 + 1, 0) + IY = range(self.nY // 2 + 1) + newMap.dataFourier[IX, IY] = dataFourier + # update real space map + newMap.data = newMap.inverseFourier() + if test: + self.plot(data) + newMap.plot(newMap.data) + return newMap.data + + ############################################################################### + # plots + + def plot(self, data=None, save=False, path=None, cmap="viridis"): + if data is None: + data = self.data.copy() + # sigma = np.std(data.flatten()) + vmin = np.min(data.flatten()) + vmax = np.max(data.flatten()) + + fig = plt.figure(0) + ax = fig.add_subplot(111) + # + # pcolor wants x and y to be edges of cell, + # ie one more element, and offset by half a cell + x = self.dX * (np.arange(self.nX + 1) - 0.5) + y = self.dY * (np.arange(self.nY + 1) - 0.5) + x, y = np.meshgrid(x, y, indexing="ij") + # + cp = ax.pcolormesh( + x * 180.0 / np.pi, y * 180.0 / np.pi, data, linewidth=0, rasterized=True + ) + # + # choose color map: jet, summer, winter, Reds, gist_gray, YlOrRd, bwr, seismic + cp.set_cmap(cmap) + # cp.set_clim(0.,255.) + # cp.set_clim(-3.*sigma, 3.*sigma) + cp.set_clim(vmin, vmax) + fig.colorbar(cp) + # + plt.axis("scaled") + ax.set_xlim(np.min(x) * 180.0 / np.pi, np.max(x) * 180.0 / np.pi) + ax.set_ylim(np.min(y) * 180.0 / np.pi, np.max(y) * 180.0 / np.pi) + ax.set_xlabel("$x$ [deg]") + ax.set_ylabel("$y$ [deg]") + # + if save == True: + if path is None: + path = "./figures/lens_simulator/" + self.name + ".pdf" + # print("saving plot to "+path) + fig.savefig(path, bbox_inches="tight") + fig.clf() + else: + plt.show() + + def plotFourier(self, dataFourier=None, save=False, name=None, cmap="viridis"): + if dataFourier is None: + dataFourier = self.dataFourier.copy() + dataFourier = np.real(dataFourier) + # sigma = np.std(dataFourier.flatten()) + vmin = np.min(dataFourier.flatten()) + vmax = np.max(dataFourier.flatten()) + + fig = plt.figure(0) + ax = fig.add_subplot(111) + # + # pcolor wants x and y to be edges of cell, + # ie one more element, and offset by half a cell + # + # left part of plot + lxLeft = 2.0 * np.pi / self.sizeX * (np.arange(-self.nX / 2 + 1, 1, 1) - 0.5) + ly = 2.0 * np.pi / self.sizeY * (np.arange(self.nY // 2 + 1 + 1) - 0.5) + lx, ly = np.meshgrid(lxLeft, ly, indexing="ij") + cp1 = ax.pcolormesh( + lx, ly, dataFourier[int(self.nX / 2 + 1):, :], linewidth=0, rasterized=True + ) + # + # right part of plot + lxRight = 2.0 * np.pi / self.sizeX * (np.arange(self.nX / 2 + 1 + 1) - 0.5) + ly = 2.0 * np.pi / self.sizeY * (np.arange(self.nY // 2 + 1 + 1) - 0.5) + lx, ly = np.meshgrid(lxRight, ly, indexing="ij") + cp2 = ax.pcolormesh( + lx, ly, dataFourier[: int(self.nX / 2 + 1), :], linewidth=0, rasterized=True + ) + # + # choose color map: jet, summer, winter, Reds, gist_gray, YlOrRd, bwr, seismic + cp1.set_cmap(cmap) + cp2.set_cmap(cmap) + # cp1.set_clim(0.,255.); cp2.set_clim(0.,255.) + # cp1.set_clim(-3.*sigma, 3.*sigma); cp2.set_clim(-3.*sigma, 3.*sigma) + cp1.set_clim(vmin, vmax) + cp2.set_clim(vmin, vmax) + # + fig.colorbar(cp1) + # + plt.axis("scaled") + ax.set_xlim(np.min(lxLeft), np.max(lxRight)) + ax.set_ylim(np.min(ly), np.max(ly)) + ax.set_xlabel(r"$\ell_x$") + ax.set_ylabel(r"$\ell_y$") + # + if save == True: + if name is None: + name = self.name + # print("saving plot to "+"./figures/lens_simulator/"+name+".pdf") + fig.savefig( + "./figures/lens_simulator/" + name + ".pdf", bbox_inches="tight" + ) + fig.clf() + else: + plt.show() + + def plotHistogram(self, data=None, save=False, name=None, nBins=100): + """plot a pixel histogram""" + if data is None: + data = self.data.copy() + + # data mean, var, skewness, kurtosis + mean = np.mean(data) + sigma = np.std(data) + # skewness = np.mean((data - mean) ** 3) / sigma**3 + # kurtosis = np.mean((data - mean) ** 4) / sigma**4 + + # print("mean =", mean) + # print("std. dev =", sigma) + # print("skewness =", skewness) + # print("kurtosis =", kurtosis) + + # data histogram, and error bars on it + pdf, binEdges = np.histogram(data, density=False, bins=nBins) + binWidth = (np.max(data) - np.min(data)) / nBins + spdf = np.sqrt(pdf) / (np.sum(pdf) * binWidth) # absolute 1-sigma error on pdf + pdf = pdf / (np.sum(pdf) * binWidth) + + # Gaussian histogram with same mean and var + samples = np.random.normal(loc=mean, scale=sigma, size=self.nX * self.nY) + pdfGauss, binEdgesGauss = np.histogram( + samples, density=True, bins=nBins, range=(binEdges[0], binEdges[-1]) + ) + + # Poisson histogram with same mean + if mean > 0.0: + samples = np.random.poisson(lam=mean, size=self.nX * self.nY) + pdfPoiss, binEdgesPoiss = np.histogram( + samples, density=True, bins=nBins, range=(binEdges[0], binEdges[-1]) + ) + + fig = plt.figure(0) + ax = fig.add_subplot(111) + # + # data histogram + ax.step(binEdges[:-1], pdf, where="post", color="b", lw=2, label=r"map") + # error band on data histogram + ax.fill_between( + np.linspace(binEdges[0], binEdges[-1], 100 * len(pdf)), + np.repeat(pdf - spdf, 100), + np.repeat(pdf + spdf, 100), + color="b", + ) + # + # Gaussian histogram + ax.step( + binEdgesGauss[:-1], + pdfGauss, + where="post", + color="r", + lw=2, + label=r"Gaussian", + ) + ax.axvline(mean, color="k") + ax.axvline(mean + 1.0 * sigma, color="gray") + ax.axvline(mean + 2.0 * sigma, color="gray") + ax.axvline(mean + 3.0 * sigma, color="gray") + ax.axvline(mean + 4.0 * sigma, color="gray") + ax.axvline(mean + 5.0 * sigma, color="gray") + ax.axvline(mean + -1.0 * sigma, color="gray") + ax.axvline(mean + -2.0 * sigma, color="gray") + ax.axvline(mean + -3.0 * sigma, color="gray") + ax.axvline(mean + -4.0 * sigma, color="gray") + ax.axvline(mean + -5.0 * sigma, color="gray") + # + # Poisson histogram + if mean > 0.0: + ax.step( + binEdgesPoiss[:-1], + pdfPoiss, + where="post", + color="g", + lw=2, + label=r"Poisson", + ) + # + ax.legend(loc=1) + # ax.set_xlim((-4.*sigma, 4.*sigma)) + ax.set_yscale("log") + # + if save == True: + if name is None: + name = self.name + # print("saving plot to "+"./figures/lens_simulator/histogram_"+name+".pdf") + fig.savefig( + "./figures/lens_simulator/histogram_" + name + ".pdf", + bbox_inches="tight", + ) + fig.clf() + else: + plt.show() + + # !!! the x and y axes might be reversed! + def plotDeflectionFieldLines(self, dx, dy): + d = np.sqrt(dx**2 + dy**2) + d /= np.max(d.flatten()) + + fig = plt.figure(0) + ax = fig.add_subplot(111) + # + strm = ax.streamplot( + self.y * 180.0 / np.pi, + self.x * 180.0 / np.pi, + dy, + dx, + color=d, + linewidth=3.0 * d, + cmap=plt.cm.autumn, + density=3.0, + ) + fig.colorbar(strm.lines) + # + plt.axis("scaled") + ax.set_xlim( + np.min(self.x.flatten()) * 180.0 / np.pi, + np.max(self.x.flatten()) * 180.0 / np.pi, + ) + ax.set_ylim( + np.min(self.y.flatten()) * 180.0 / np.pi, + np.max(self.y.flatten()) * 180.0 / np.pi, + ) + ax.set_xlabel("$x$ [deg]") + ax.set_ylabel("$y$ [deg]") + + plt.show() + + def plotDeflectionArrows(self, dx, dy, save=False, name=None): + # d = np.sqrt(dx**2 + dy**2) + + # Keep only some of the points + skip = (slice(None, None, self.nX / 25), slice(None, None, self.nX / 25)) + + fig = plt.figure(0) + ax = fig.add_subplot(111) + # + # + # pcolor wants x and y to be edges of cell, + # ie one more element, and offset by half a cell + x = self.dX * (np.arange(self.nX + 1) - 0.5) + y = self.dY * (np.arange(self.nY + 1) - 0.5) + x, y = np.meshgrid(x, y, indexing="ij") + # + # sigma = self.data.std() + cp = ax.pcolormesh( + x * 180.0 / np.pi, + y * 180.0 / np.pi, + self.data, + linewidth=0, + rasterized=True, + alpha=1, + cmap=plt.cm.jet, + ) + fig.colorbar(cp) + # + # + ax.quiver( + self.x[skip] * 180.0 / np.pi, + self.y[skip] * 180.0 / np.pi, + dx[skip] * 180.0 / np.pi, + dy[skip] * 180.0 / np.pi, + facecolor="k", + units="xy", + angles="xy", + scale_units="xy", + scale=1, + width=0.003 * self.sizeX * 180.0 / np.pi, + ) + # + # + plt.axis("scaled") + ax.set_xlim( + np.min(self.x.flatten()) * 180.0 / np.pi, + np.max(self.x.flatten()) * 180.0 / np.pi, + ) + ax.set_ylim( + np.min(self.y.flatten()) * 180.0 / np.pi, + np.max(self.y.flatten()) * 180.0 / np.pi, + ) + ax.set_xlabel("$x$ [deg]") + ax.set_ylabel("$y$ [deg]") + # + if save == True: + if name is None: + name = self.name + # print "saving plot to "+"./figures/lens_simulator/"+name+".pdf" + fig.savefig( + "./figures/lens_simulator/" + name + "_arrows.pdf", bbox_inches="tight" + ) + fig.clf() + else: + plt.show() + + ############################################################################### + # Fourier transforms, notmalized such that + # f(k) = int dx e-ikx f(x) + # f(x) = int dk/2pi eikx f(k) + + def fourier(self, data=None): + """Fourier transforms, notmalized such that + f(k) = int dx e-ikx f(x) + f(x) = int dk/2pi eikx f(k) + """ + if data is None: + data = self.data.copy() + result = np.fft.rfftn(data) + # # use pyfftw's fft. Make sure the real-space data has type np.float128 + # result = pyfftw.interfaces.numpy_fft.rfftn((np.float128)(data)) + result *= self.dX * self.dY # NOTE for precision: this leads to ~1e-3 % + # # differences in outer radial bins of projection + # print(self.dX,self.dY) + return result + + def inverseFourier(self, dataFourier=None): + """Fourier transforms, notmalized such that + f(k) = int dx e-ikx f(x) + f(x) = int dk/2pi eikx f(k) + """ + if dataFourier is None: + dataFourier = self.dataFourier.copy() + # use numpy's fft + result = np.fft.irfftn(dataFourier) + # # use pyfftw's fft. Make sure the Fourier data has type np.complex128 + # result = pyfftw.interfaces.numpy_fft.irfftn((np.complex128)(dataFourier)) + result /= self.dX * self.dY + return result + + ############################################################################### + # Measure power spectrum + + def crossPowerSpectrum( + self, + dataFourier1, + dataFourier2, + theory=[], + fsCl=None, + nBins=51, + lRange=None, + plot=False, + name="test", + save=False, + ): + + # define ell bins + ell = self.ell.flatten() + if lRange is None: + lEdges = np.logspace(np.log10(1.0), np.log10(np.max(ell)), nBins, 10.0) + else: + lEdges = np.logspace(np.log10(lRange[0]), np.log10(lRange[-1]), nBins, 10.0) + + # bin centers + lCen, lEdges, binIndices = stats.binned_statistic( + ell, ell, statistic="mean", bins=lEdges + ) + # when bin is empty, replace lCen by a naive expectation + lCenNaive = 0.5 * (lEdges[:-1] + lEdges[1:]) + lCen[np.where(np.isnan(lCen))] = lCenNaive[np.where(np.isnan(lCen))] + # number of modes + Nmodes, lEdges, binIndices = stats.binned_statistic( + ell, np.zeros_like(ell), statistic="count", bins=lEdges + ) + Nmodes = np.nan_to_num(Nmodes) + # power spectrum + power = (dataFourier1 * np.conj(dataFourier2)).flatten() + power = np.real( + power + ) # unnecessary in principle, but avoids binned_statistics to complain + Cl, lEdges, binIndices = stats.binned_statistic( + ell, power, statistic="mean", bins=lEdges + ) + Cl = np.nan_to_num(Cl) + # finite volume correction + Cl /= self.sizeX * self.sizeY + # 1sigma uncertainty on Cl + if fsCl is None: + sCl = Cl * np.sqrt(2) + else: + sCl = np.array(map(fsCl, lCen)) + # In case of a cross-correlation, Cl may be negative. + # the absolute value is then still some estimate of the error bar + sCl = np.abs(sCl) + sCl /= np.sqrt(Nmodes) + sCl[np.where(np.isfinite(sCl) == False)] = 0.0 + + if plot: + factor = 1.0 # lCen**2 + + fig = plt.figure(0) + ax = fig.add_subplot(111) + # + Ipos = np.where(Cl >= 0.0) + Ineg = np.where(Cl < 0.0) + ax.errorbar( + lCen[Ipos], factor * Cl[Ipos], yerr=factor * sCl[Ipos], c="b", fmt="." + ) + ax.errorbar( + lCen[Ineg], -factor * Cl[Ineg], yerr=factor * sCl[Ineg], c="r", fmt="." + ) + # + for f in theory: + L = np.logspace(np.log10(1.0), np.log10(np.max(ell)), 201, 10.0) + ClExpected = np.array(map(f, L)) + ax.plot(L, factor * ClExpected, "k") + # + # ax.axhline(0.) + ax.set_xscale("log", nonposx="clip") + ax.set_yscale("log", nonposy="clip") + # ax.set_xlim(1.e1, 4.e4) + # ax.set_ylim(1.e-5, 2.e5) + ax.set_xlabel(r"$\ell$") + # ax.set_ylabel(r'$\ell^2 C_\ell$') + ax.set_ylabel(r"$C_\ell$") + # + if save == True: + if name is None: + name = self.name + # print "saving plot to "+"./figures/lens_simulator/"+name+"_power.pdf" + fig.savefig( + "./figures/lens_simulator/" + name + "_power.pdf", + bbox_inches="tight", + ) + fig.clf() + else: + plt.show() + + return lCen, Cl, sCl + + def powerSpectrum( + self, + dataFourier=None, + theory=[], + fsCl=None, + nBins=51, + lRange=None, + plot=False, + name="test", + save=False, + ): + if dataFourier is None: + dataFourier = self.dataFourier.copy() + return self.crossPowerSpectrum( + dataFourier1=dataFourier, + dataFourier2=dataFourier, + theory=theory, + fsCl=fsCl, + nBins=nBins, + lRange=lRange, + plot=plot, + name=name, + save=save, + ) + + def binTheoryPowerSpectrum(self, fCl, nBins=17, lRange=None): + """Bin a theory power spectrum to allow to compare it + with the measured power spectrum of a map.""" + # define ell bins + ell = self.ell.flatten() + if lRange is None: + lEdges = np.logspace(np.log10(1.0), np.log10(np.max(ell)), nBins, 10.0) + else: + lEdges = np.logspace(np.log10(lRange[0]), np.log10(lRange[-1]), nBins, 10.0) + + # bin centers + lCen, lEdges, binIndices = stats.binned_statistic( + ell, ell, statistic="mean", bins=lEdges + ) + # when bin is empty, replace lCen by a naive expectation + lCenNaive = 0.5 * (lEdges[:-1] + lEdges[1:]) + lCen[np.where(np.isnan(lCen))] = lCenNaive[np.where(np.isnan(lCen))] + + # generate map with theory power spectrum + clmapFourier = self.filterFourierIsotropic( + fCl, dataFourier=np.ones_like(self.ell), test=False + ) + clmapFourier = np.real(clmapFourier.flatten()) + Cl, lEdges, binIndices = stats.binned_statistic( + ell, clmapFourier, statistic="mean", bins=lEdges + ) + Cl = np.nan_to_num(Cl) + return lCen, Cl + + ############################################################################### + # Gaussians and tests for Fourier transform conventions + + def genGaussian(self, meanX=0.0, meanY=0.0, sigma1d=1.0): + result = np.exp( + -0.5 * ((self.x - meanX) ** 2 + (self.y - meanY) ** 2) / sigma1d**2 + ) + result /= 2.0 * np.pi * sigma1d**2 + return result + + def genGaussianFourier(self, meanLX=0.0, meanLY=0.0, sigma1d=1.0): + result = np.exp( + -0.5 * ((self.lx - meanLX) ** 2 + (self.ly - meanLY) ** 2) / sigma1d**2 + ) + result /= 2.0 * np.pi * sigma1d**2 + return result + + def testFourierGaussian(self): + """tests that the FT of a Gaussian is a Gaussian, + with correct normalization and variance + """ + # generate a quarter of a Gaussian + sigma1d = self.sizeX / 10.0 + self.data = self.genGaussian(sigma1d=sigma1d) + # show it + self.plot() + + # fourier transform it + self.dataFourier = self.fourier() + self.plotFourier() + + # computed expected Gaussian + expectedFourier = self.genGaussianFourier(sigma1d=1.0 / sigma1d) + expectedFourier *= 2.0 * np.pi * (1.0 / sigma1d) ** 2 + expectedFourier /= 4.0 # because only one quadrant in real space + + # self.plotFourier(data=self.dataFourier/expectedFourier-1.) + + # compare along one axis + plt.plot(self.dataFourier[0, :], "k") + plt.plot(expectedFourier[0, :], "r") + plt.show() + + # compare along other axis + plt.plot(self.dataFourier[:, 0], "k") + plt.plot(expectedFourier[:, 0], "r") + plt.show() + + def testFourierCos(self): + """tests that the FT of cos(k*x) + peaks at the right k + """ + # generate a quarter of a Gaussian + ell = 100.0 + self.data = np.cos(ell * self.x) + np.cos(ell * self.y) + # show it + self.plot() + + # fourier transform it + self.dataFourier = self.fourier() + # self.plotFourier() + + self.plotFourier() + + def testInverseFourier(self): + """test that the inverse FT of the forward FT is the initial function""" + # generate a quarter of a Gaussian + sigma1d = 2.0 + self.data = self.genGaussian(sigma1d=sigma1d) + # show it + self.plot() + + # Fourier transform it + self.dataFourier = self.fourier() + # inverse Fourier transform it + expectedData = self.inverseFourier() + + # compare along each axis + plt.plot(self.y[0, :], self.data[0, :] / expectedData[0, :] - 1.0, "g") + plt.plot(self.x[:, 0], self.data[:, 0] / expectedData[:, 0] - 1.0, "b--") + plt.show() + + ############################################################################### + # generate Gaussian random field with any power spectrum + + def genGRF(self, fCl, test=False): + + # generate Gaussian white noise in real space + data = np.zeros_like(self.data) + data = np.random.normal( + loc=0.0, scale=1.0 / np.sqrt(self.dX * self.dY), size=len(self.x.flatten()) + ) + data = data.reshape(np.shape(self.x)) + + # Fourier transform + dataFourier = self.fourier(data) + if test: + # check that the power spectrum is Cl = 1 + self.powerSpectrum(dataFourier, theory=[lambda ell: 1.0], plot=True) + + # multiply by desired power spectrum + f = lambda ell: np.sqrt(fCl(ell)) + clFourier = np.array(map(f, self.ell.flatten())) + clFourier = np.nan_to_num(clFourier) + clFourier = clFourier.reshape(np.shape(self.ell)) + dataFourier *= clFourier + if test: + # check 0 mode + # print "l=0 mode is:", dataFourier[0,0] + # check that the power spectrum is the desired one + self.powerSpectrum(dataFourier, theory=[fCl], plot=True) + # show the fourier map + self.plotFourier(dataFourier) + # show the real space map + data = self.inverseFourier(dataFourier) + self.plot(data) + + return dataFourier + + def saveGRFMocks(self, fCl, nRand, directory=None, name=None): + """create nRand GRF mock maps""" + if directory is None: + directory = "./output/lens_simulator/mocks/" + self.name + if name is None: + name = "mock_" + # create folder if needed + if not os.path.exists(directory): + os.makedirs(directory) + + for iRand in range(nRand): + path = directory + "/" + name + str(iRand) + ".fits" + randomFourier = self.genGRF(fCl) + self.saveDataFourier(randomFourier, path) + + def loadAllGRFMocks(self, nRand, directory=None, name=None): + """DO NOT USE unless you have infinite RAM ;)""" + if directory is None: + directory = "./output/lens_simulator/mocks/" + self.name + if name is None: + name = "mock_" + + self.mockDataFourier = {} + for iRand in range(nRand): + path = directory + "/" + name + str(iRand) + ".fits" + self.mockDataFourier[iRand] = self.loadDataFourier(path) + + def genCorrGRF(self, fC11, fC22, fC12, test=False): + """Generate two correlated GRFs, with the correct auto and cross spectra.""" + # mao 1: generate GRF + data1Fourier = self.genGRF(fC11, test=False) + + # map 2: start with part correlated with map 1 + f = lambda ell: fC12(ell) / fC11(ell) + data2Fourier = self.filterFourierIsotropic( + f, dataFourier=data1Fourier, test=False + ) + # map 2: add uncorrelated part + f = lambda ell: fC22(ell) - (fC12(ell) * fC12(ell)) / (fC11(ell)) + data2Fourier += self.genGRF(f, test=False) + # avoid nan and inf + data1Fourier[np.where(np.isfinite(data1Fourier) == False)] = 0.0 + data2Fourier[np.where(np.isfinite(data2Fourier) == False)] = 0.0 + + if test: + self.powerSpectrum(data1Fourier, theory=[fC11], plot=True) + self.powerSpectrum(data2Fourier, theory=[fC22], plot=True) + self.crossPowerSpectrum( + data1Fourier, data2Fourier, theory=[fC12], plot=True + ) + + return data1Fourier, data2Fourier + + ############################################################################### + # generate map where T_ell has C_ell for modulus square, + # and a random phase. + # Useful to avoid sample variance + + def genMockIsotropicNoSampleVar(self, fCl, test=False, path=None): + # generate map with correct modulus + f = lambda ell: np.sqrt(fCl(ell)) + resultFourier = self.filterFourierIsotropic( + f, dataFourier=np.ones_like(self.dataFourier), test=test + ) + resultFourier *= np.sqrt(self.sizeX * self.sizeY) + + # generate random phases + f = lambda lx, ly: np.exp(1j * np.random.uniform(0.0, 2.0 * np.pi)) + resultFourier = self.filterFourier(f, dataFourier=resultFourier) + + # keep it real ;) + # if lx=ly=0, set to zero + resultFourier[0, 0] = 0.0 + # if ly=0, make sure T[lx, 0] = T[-lx, 0]^* + for i in range(self.nX // 2 + 1, self.nX): + resultFourier[i, 0] = np.conj(resultFourier[self.nX - i, 0]) + + # save if needed + if path is not None: + self.saveDataFourier(resultFourier, path) + + return resultFourier + + ############################################################################### + # Generate Poisson white noise map + + def genPoissonWhiteNoise(self, nbar, norm=False, test=False): + """Generate Poisson white noise. + Returns real space map. + nbar mean number density of objects, in 1/sr. + if norm=False, returns a map of N (number) + if norm=True, returns a map of delta = N/Nbar - 1 + """ + # number of objects per pixel + Ngal = nbar * self.dX * self.dY + if test: + print("generate Poisson white noise") + # print nbar+" objects per sr, i.e. "+Ngal+" objects per pixel" + data = np.random.poisson(lam=Ngal, size=len(self.x.flatten())) + if norm: + data = data / Ngal - 1.0 + data = data.reshape(np.shape(self.x)) + return data + + ############################################################################### + # filter map + + def filterFourierIsotropic(self, fW, dataFourier=None, test=False): + """the filter fW is assumed to be function of |ell|""" + if dataFourier is None: + dataFourier = self.dataFourier.copy() + W = np.array(map(fW, self.ell.flatten())) + W = W.reshape(self.ell.shape) + if test: + self.plotFourier(dataFourier=W) + # + plt.plot(self.ell.flatten(), W.flatten(), "b.") + plt.show() + + result = dataFourier * W + result = np.nan_to_num(result) + return result + + def filterFourier(self, fW, dataFourier=None, test=False): + """the filter fW is assumed to be function of lx, ly""" + if dataFourier is None: + dataFourier = self.dataFourier.copy() + + f = np.vectorize(fW) + W = f(self.lx, self.ly) + if test: + self.plotFourier(dataFourier=W) + # + plt.plot(self.ell.flatten(), W.flatten(), "b.") + plt.show() + + result = dataFourier * W + result = np.nan_to_num(result) + return result + + ############################################################################### + # Matched filter and point source mask + + def matchedFilterIsotropic(self, fCl, fprof=None, dataFourier=None, test=False): + if dataFourier is None: + dataFourier = self.dataFourier.copy() + if fprof is None: + fprof = lambda ell: 1.0 + + # filter function + def fW(ell): + result = fprof(ell) / fCl(ell) + if not np.isfinite(result): + result = 0.0 + return result + + # filter the map + resultFourier = self.filterFourierIsotropic( + fW, dataFourier=dataFourier, test=test + ) + + # normalization function + def fNorm(ell): + result = ell / (2.0 * np.pi) * fprof(ell) ** 2 / fCl(ell) + if not np.isfinite(result): + result = 0.0 + return result + + # compute normalization + normalization = integrate.quad( + fNorm, self.ell.min(), self.ell.max(), epsabs=0.0, epsrel=1.0e-3 + )[0] + + # normalize the filtered map + resultFourier /= normalization + if test: + result = self.inverseFourier(resultFourier) + self.plot(result) + return resultFourier + + def pointSourceMaskMatchedFilterIsotropic( + self, + fCl, + fluxCut, + fprof=None, + dataFourier=None, + maskPatchRadius=None, + test=False, + ): + """Returns the mask for point sources with flux above fluxCut. + If T(x) in Jy/sr, T_l in Jy, Cl in Jy^2/sr, fluxCut in Jy. + If T(x) in muK, T_l in muK*sr, Cl in (muK)^2*sr, fluxCut in muK*sr. + prof_l is the profile before beam convolution (i.e. 1 for point source). + The mask is 1 almost everywhere, and 0 on point sources. + Patch is patch radius in rad, to mask around point sources + """ + # matched-filter the map + filteredFourier = self.matchedFilterIsotropic( + fCl, fprof=fprof, dataFourier=dataFourier, test=test + ) + filtered = self.inverseFourier(filteredFourier) + + # threshold the map + mask = 1.0 * (np.abs(filtered) < fluxCut) + + if maskPatchRadius is not None: + # make a Gaussian such that the fwhm is twice the maskPatchRadius + fwhm = 2.0 * maskPatchRadius # 5. * np.pi/(180.*60.) + s = fwhm / np.sqrt(8.0 * np.log(2.0)) + f = lambda ell: np.exp(-0.5 * ell**2 * s**2) + gaussFourier = self.filterFourierIsotropic( + f, dataFourier=np.ones_like(self.dataFourier) + ) + + # normalize properly, accounting for Kronecker vs Dirac and finite pixel size + gaussFourier /= special.erf(self.dX / (np.sqrt(8.0) * s)) * special.erf( + self.dY / (np.sqrt(8.0) * s) + ) + + # smooth the mask + # Fourier transform 1-mask, to have + # zero everywhere except on the point sources + maskFourier = self.fourier(1.0 - mask) + mask = self.inverseFourier(maskFourier * gaussFourier) + + # threshold the smoothed mask at half-max + mask = 1.0 * (mask < 0.5) # *np.max(mask.flatten())) + + return mask + + ############################################################################### + + ############################################################################### + + def randomizePhases(self, dataFourier=None, test=False): + """Generate new map with same Fourier amplitudes, + but replacing the original phases + with iid uniformly distributed Fourier phases. + """ + if dataFourier is None: + dataFourier = self.dataFourier.copy() + + f = lambda z: np.abs(z) * np.exp(1j * np.random.uniform(0.0, 2.0 * np.pi)) + resultFourier = np.array(map(f, dataFourier.flatten())) + resultFourier = resultFourier.reshape(dataFourier.shape) + return resultFourier + + ############################################################################### + + def filterCollapsed4PtFunc(self, ell, lMean): + # lsigma = 50. + # lsigma = 100. + # lsigma = 200. + # lsigma = 300. + lsigma = 500.0 + result = np.exp(-((ell - lMean) ** 2) / (4.0 * lsigma**2)) + result *= (2.0 * np.pi / lsigma**2) ** 0.25 / np.sqrt(lMean) + # result *= (l>=lMean - 4.*lsigma) * (l<=lMean + 4.*lsigma) + return result + + def collapsed4PtFunc( + self, + lMean=1000.0, + dataFourier=None, + theory=None, + name="test", + save=False, + nBins=51, + lRange=None, + ): + if dataFourier is None: + dataFourier = self.dataFourier.copy() + + # # define filter function + # lMean = 1.e3 + # sl = 2.e2 + # # defined so that int d^2l W^2 = 1? + # # !!! normalization is wrong here, but doesn't matter for ratios + # fW2 = lambda l: np.exp(-0.5*(l-lMean)**2/sl**2) / (2.*np.pi*sl**2) + # fW = lambda l: np.sqrt(fW2(l)) + + # filter in Fourier space, to keep only l,l' \simeq lmean + f = lambda ell: self.filterCollapsed4PtFunc(ell, lMean) + dataFourier = self.filterFourierIsotropic(f, dataFourier) + # update real space map + data = self.inverseFourier(dataFourier=dataFourier) + # square map in real space + data = data**2 + # update Fourier map + dataFourier = self.fourier(data=data) + + # compute power spectrum of squared map + lCen, Cl, sCl = self.powerSpectrum( + dataFourier=dataFourier, + theory=theory, + name=name, + save=save, + nBins=nBins, + lRange=lRange, + ) + return lCen, Cl, sCl + + def saveTrispectrum( + self, + dataFourier=None, + gaussDataFourier=None, + path="./output/flat_map/test_", + nBins=51, + lRange=None, + ): + # print "computing trispectrum" + + # array of the high map ells to evaluate + lMean = np.logspace(np.log10(100.0), np.log10(3.0e3), 9, 10.0) # 9 + + # choose the potentially non-Gaussian map to analyze + if dataFourier is None: + dataFourier = self.dataFourier.copy() + + # analyze the potentially non-Gaussian map + collapsed4PtFunc_NG = {} + sCollapsed4PtFunc_NG = {} + for ilMean in range(len(lMean)): + # here L is the low ell, while lMean is the high ell + ( + L, + collapsed4PtFunc_NG[ilMean], + sCollapsed4PtFunc_NG[ilMean], + ) = self.collapsed4PtFunc( + lMean=lMean[ilMean], dataFourier=dataFourier, nBins=nBins, lRange=lRange + ) + # print "done "+str(ilMean+1)+" of "+str(len(lMean)) + + # analyze the Gaussian mock with the same power spectrum + collapsed4PtFunc_G = {} + sCollapsed4PtFunc_G = {} + for ilMean in range(len(lMean)): + # here L is the low ell, while lMean is the high ell + ( + L, + collapsed4PtFunc_G[ilMean], + sCollapsed4PtFunc_G[ilMean], + ) = self.collapsed4PtFunc( + lMean=lMean[ilMean], + dataFourier=gaussDataFourier, + nBins=nBins, + lRange=lRange, + ) + # print "done "+str(ilMean+1)+" of "+str(len(lMean)) + + # the 4pt func is 2 * C^2/Nmodes + trispectrum. + # Infer C^2/Nmodes, C^2 and trispectrum + C2Nmodes = np.zeros((len(L), len(lMean))) + C2 = np.zeros((len(L), len(lMean))) + Trispec = np.zeros((len(L), len(lMean))) + # + sC2Nmodes = np.zeros((len(L), len(lMean))) + sC2 = np.zeros((len(L), len(lMean))) + sTrispec = np.zeros((len(L), len(lMean))) + for ilMean in range(len(lMean)): + # + # Number of modes included in filter around lMean + f = ( + lambda lnl: np.exp(lnl) ** 2 + / (2.0 * np.pi) + * self.filterCollapsed4PtFunc(np.exp(lnl), lMean[ilMean]) ** 2 + ) + # test = integrate.quad(f, np.log(1.), np.log(1.e5), epsabs=0., epsrel=1.e-3) + ##print + ##print test[0], test[1] / test[0] + Nmodes = integrate.quad( + f, np.log(1.0), np.log(1.0e5), epsabs=0.0, epsrel=1.0e-3 + )[0] + Nmodes **= 2 + f = ( + lambda lnl: np.exp(lnl) ** 2 + / (2.0 * np.pi) + * self.filterCollapsed4PtFunc(np.exp(lnl), lMean[ilMean]) ** 4 + ) + # test = integrate.quad(f, np.log(1.), np.log(1.e5), epsabs=0., epsrel=1.e-3) + ##print test[0], test[1] / test[0] + Nmodes /= integrate.quad( + f, np.log(1.0), np.log(1.0e5), epsabs=0.0, epsrel=1.0e-3 + )[0] + # + C2Nmodes[:, ilMean] = collapsed4PtFunc_G[ilMean].copy() + sC2Nmodes[:, ilMean] = sCollapsed4PtFunc_G[ilMean].copy() + Trispec[:, ilMean] = ( + collapsed4PtFunc_NG[ilMean] - collapsed4PtFunc_G[ilMean] + ) + sTrispec[:, ilMean] = np.sqrt(2.0) * sC2Nmodes[:, ilMean] + # + C2[:, ilMean] = C2Nmodes[:, ilMean] * Nmodes + sC2[:, ilMean] = sC2Nmodes[:, ilMean] * Nmodes + + # save everything + # print "saving trispectrum to "+path + np.savetxt(path + "_Llow.txt", L) + np.savetxt(path + "_lhigh.txt", lMean) + np.savetxt(path + "_C2Nmodes.txt", C2Nmodes) + np.savetxt(path + "_sC2Nmodes.txt", sC2Nmodes) + np.savetxt(path + "_C2.txt", C2) + np.savetxt(path + "_sC2.txt", sC2) + np.savetxt(path + "_Trispec.txt", Trispec) + np.savetxt(path + "_sTrispec.txt", sTrispec) + + def loadTrispectrum(self, path="./output/flat_map/test_"): + # read everything + # print "loading trispectrum from "+path + L = np.genfromtxt(path + "_Llow.txt") + ell = np.genfromtxt(path + "_lhigh.txt") + C2Nmodes = np.genfromtxt(path + "_C2Nmodes.txt") + sC2Nmodes = np.genfromtxt(path + "_sC2Nmodes.txt") + C2 = np.genfromtxt(path + "_C2.txt") + sC2 = np.genfromtxt(path + "_sC2.txt") + Trispec = np.genfromtxt(path + "_Trispec.txt") + sTrispec = np.genfromtxt(path + "_sTrispec.txt") + + return L, ell, C2Nmodes, sC2Nmodes, C2, sC2, Trispec, sTrispec + + ############################################################################### + + def computeGradient(self, dataFourier=None): + """returns the Fourier maps of the 2 components of the gradient""" + if dataFourier is None: + dataFourier = self.dataFourier.copy() + dx = dataFourier * 1.0j * self.lx + dy = dataFourier * 1.0j * self.ly + return dx, dy + + def computeDivergence(self, dxFourier, dyFourier): + divFourier = dxFourier * 1.0j * self.lx + dyFourier * 1.0j * self.ly + return divFourier + + def computeLaplacian(self, dataFourier=None): + if dataFourier is None: + dataFourier = self.dataFourier.copy() + laplacianFourier = -self.ell**2 * dataFourier + return laplacianFourier + + def phiFromKappa(self, kappaFourier): + phiFourier = -2.0 * kappaFourier / self.ell**2 + phiFourier[np.where(np.isfinite(phiFourier) == False)] = 0.0 + return phiFourier + + def kappaFromPhi(self, phiFourier): + """get kappa from phi, with: + Delta phi = 2 kappa + """ + kappaFourier = -0.5 * self.ell**2 * phiFourier + return kappaFourier + + def deflectionFromKappa(self, kappaFourier, test=False): + """gives the two components of the lensing deflection + should be called from a kappa map + Delta phi = 2 kappa + d = grad phi + """ + dxFourier = -2.0j * self.lx / self.ell**2 * kappaFourier + dyFourier = -2.0j * self.ly / self.ell**2 * kappaFourier + # + dxFourier[np.where(np.isfinite(dxFourier) == False)] = 0.0 + dyFourier[np.where(np.isfinite(dyFourier) == False)] = 0.0 + # + dx = self.inverseFourier(dxFourier) + dy = self.inverseFourier(dyFourier) + # + if test: + self.plot(dx) + self.plot(dy) + return dx, dy + + # def deflectionFromPhi(self, phiFourier): + # """gives the two components of the lensing deflection + # should be called from a phi map + # d = grad phi + # """ + # dxFourier = 1j * phiMap.lx * phiFourier + # dyFourier = 1j * phiMap.ly * phiFourier + # # + # dxFourier = np.nan_to_num(dxFourier) + # dyFourier = np.nan_to_num(dyFourier) + # # + # dx = self.inverseFourier(dxFourier) + # dy = self.inverseFourier(dyFourier) + # return dx, dy + + ############################################################################### + + + ############################################################################### + ############################################################################### + # Non-normalized quadratic estimator + + + def quadEstPhiNonNorm( + self, + fC0, + fCtot, + lMin=1.0, + lMax=1.0e5, + dataFourier=None, + dataFourier2=None, + test=False, + ): + """non-normalized quadratic estimator + fC0: ulensed power spectrum + fCtot: lensed power spectrum + noise + """ + if dataFourier is None: + dataFourier = self.dataFourier.copy() + if dataFourier2 is None: + dataFourier2 = dataFourier.copy() + + # inverse-var weighted map + def f(ell): + # cut off the high ells from input map + if (ell < lMin) or (ell > lMax): + return 0.0 + result = 1.0 / fCtot(ell) + if not np.isfinite(result): + result = 0.0 + return result + + iVarDataFourier = self.filterFourierIsotropic( + f, dataFourier=dataFourier, test=test + ) + iVarData = self.inverseFourier(iVarDataFourier) + if test: + # print "showing the inverse var. weighted map" + self.plot(data=iVarData) + # print "checking the power spectrum of this map" + self.powerSpectrum(theory=f, dataFourier=iVarDataFourier, plot=True) + + # Wiener-filter the map + def f(ell): + # cut off the high ells from input map + if (ell < lMin) or (ell > lMax): + return 0.0 + result = fC0(ell) / fCtot(ell) + if not np.isfinite(result): + result = 0.0 + return result + + WFDataFourier = self.filterFourierIsotropic( + f, dataFourier=dataFourier2, test=test + ) + if test: + # print "showing the WF map" + WFData = self.inverseFourier(WFDataFourier) + self.plot(data=WFData) + # print "checking the power spectrum of this map" + theory = lambda ell: f(ell) * fC0(ell) + self.powerSpectrum(theory=theory, dataFourier=WFDataFourier, plot=True) + + # get Wiener-filtered gradient map + WFDataXFourier, WFDataYFourier = self.computeGradient(dataFourier=WFDataFourier) + WFDataX = self.inverseFourier(dataFourier=WFDataXFourier) + WFDataY = self.inverseFourier(dataFourier=WFDataYFourier) + if test: + # print "showing x gradient of WF map" + self.plot(data=WFDataX) + # print "checking power spectrum of this map" + theory = ( + lambda ell: 0.5 * ell**2 * f(ell) * fC0(ell) + ) # 0.5 is from average of cos^2 + self.powerSpectrum(theory=theory, dataFourier=WFDataXFourier, plot=True) + # print "showing y gradient of WF map" + self.plot(data=WFDataY) + # print "checking power spectrum of this map" + theory = ( + lambda ell: 0.5 * ell**2 * f(ell) * fC0(ell) + ) # 0.5 is from average of sin^2 + self.powerSpectrum(theory=theory, dataFourier=WFDataYFourier, plot=True) + + # product in real space + productDataX = iVarData * WFDataX + productDataY = iVarData * WFDataY + productDataXFourier = self.fourier(data=productDataX) + productDataYFourier = self.fourier(data=productDataY) + + # take divergence + divergenceDataFourier = self.computeDivergence( + productDataXFourier, productDataYFourier + ) + if test: + # print "showing divergence map" + divergenceData = self.inverseFourier(dataFourier=divergenceDataFourier) + self.plot(data=divergenceData) + # print "checking the power spectrum of divergence map" + self.powerSpectrum(dataFourier=divergenceDataFourier, plot=True) + + # cut off the high ells from phi map + f = lambda ell: (ell <= 2.0 * lMax) + divergenceDataFourier = self.filterFourierIsotropic( + f, dataFourier=divergenceDataFourier, test=test + ) + + return divergenceDataFourier + + ############################################################################### + # Normalization for quadratic estimator, i.e. N_l^{0 phiphi} + + def computeQuadEstPhiNormalizationAna(self, fN_phi_TT, test=False): + """the normalization is N_l^phiphi, + obtained by evaluating the analytical calculation for the reconstruction noise + """ + W = np.array(map(fN_phi_TT, self.ell.flatten())) + W = np.nan_to_num(W) + W = W.reshape(self.ell.shape) + if test: + self.plotFourier(dataFourier=W) + # + plt.plot(self.ell.flatten(), W.flatten(), "b.") + plt.show() + return W + + def computeQuadEstPhiNormalizationRand( + self, + fC0, + fCtot, + lMin=1.0, + lMax=1.0e5, + nRand=1, + path=None, + dataFourier=None, + test=False, + ): + """the normalization is N_l^phiphi, + obtained as the inverse of the average power of the quad est on randoms, + nRand: number of random realizations to average + """ + if dataFourier is None: + dataFourier = self.dataFourier.copy() + + result = np.zeros_like(self.dataFourier) + + for iRand in range(nRand): + # # randomize the phases + # randomFourier = self.randomizePhases(dataFourier=dataFourier) + if path is not None: + # print "read mock", iRand, "from", path+str(iRand)+".fits" + randomFourier = self.loadDataFourier(path + str(iRand) + ".fits") + else: + # print "generate mock", iRand + randomFourier = self.genGRF(fCtot, test=test) + + # get the non-normalized quad. est. for this map + randomQuadFourier = self.quadEstPhiNonNorm( + fC0, fCtot, lMin=lMin, lMax=lMax, dataFourier=randomFourier, test=test + ) + # measure its power spectrum, + # with enough bins... + lCen, Cl, sCl = self.powerSpectrum( + dataFourier=randomQuadFourier, nBins=0.1 * self.nX, plot=test + ) + # interpolate the inverse power spectrum + f = UnivariateSpline(lCen, Cl, k=1, s=0) + + if test: + # see if we used too many/few bins for power spectrum + L = np.linspace(0.0, 8.0e3, 10001) + F = np.array(map(f, L)) + plt.loglog(L, F, "b") + plt.loglog(L, -F, "r") + plt.show() + + # make it the normalization map + W = np.array(map(f, self.ell.flatten())) + W = np.nan_to_num(1.0 / W) + W = W.reshape(self.ell.shape) + result += W + result /= nRand + + if test: + # print "showing normalization map" + self.plotFourier(dataFourier=result) + # print "checking power spectrum" + self.powerSpectrum(dataFourier=result) + + return result + + def computeQuadEstPhiNormalizationFFT( + self, fC0, fCtot, lMin=1.0, lMax=1.0e5, test=False, cache=None + ): + """the normalization is N_l^phiphi, + Works great, and super fast. + """ + + # Actual calculation + def doCalculation(): + # print "Doing full calculation: computeQuadEstPhiNormalizationFFT" + # inverse-var weighted map + def f(ell): + if (ell < lMin) or (ell > lMax): + return 0.0 + result = 1.0 / fCtot(ell) + if not np.isfinite(result): + result = 0.0 + return result + + iVarFourier = np.array(map(f, self.ell.flatten())) + iVarFourier = iVarFourier.reshape(self.ell.shape) + iVar = self.inverseFourier(dataFourier=iVarFourier) + + # C map + def f(ell): + if (ell < lMin) or (ell > lMax): + return 0.0 + result = fC0(ell) ** 2 / fCtot(ell) + if not np.isfinite(result): + result = 0.0 + return result + + CFourier = np.array(map(f, self.ell.flatten())) + CFourier = CFourier.reshape(self.ell.shape) + + # term 1x + term1x = self.inverseFourier(dataFourier=self.lx**2 * CFourier) + term1x *= iVar + term1xFourier = self.fourier(data=term1x) + term1xFourier *= self.lx**2 + # + # term 1y + term1y = self.inverseFourier(dataFourier=self.ly**2 * CFourier) + term1y *= iVar + term1yFourier = self.fourier(data=term1y) + term1yFourier *= self.ly**2 + # + # term 1xy + term1xy = self.inverseFourier( + dataFourier=2.0 * self.lx * self.ly * CFourier + ) + term1xy *= iVar + term1xyFourier = self.fourier(data=term1xy) + term1xyFourier *= self.lx * self.ly + + if test: + self.plotFourier(term1xFourier) + self.plotFourier(term1yFourier) + self.plotFourier(term1xyFourier) + self.plotFourier(term1xFourier + term1yFourier + term1xyFourier) + + # WF map + def f(ell): + if (ell < lMin) or (ell > lMax): + return 0.0 + result = fC0(ell) / fCtot(ell) + # artificial factor of i such that f(-l) = f(l)*, + # such that f(x) is real + result *= 1.0j + if not np.isfinite(result): + result = 0.0 + return result + + WFFourier = np.array(map(f, self.ell.flatten())) + WFFourier = WFFourier.reshape(self.ell.shape) + + # term 2 + term2_x = self.inverseFourier(dataFourier=self.lx * WFFourier) + term2_y = self.inverseFourier(dataFourier=self.ly * WFFourier) + # + # term 2x + term2x = term2_x**2 + term2xFourier = self.fourier(data=term2x) + term2xFourier *= self.lx**2 + # minus sign to correct the artificial i**2 + term2xFourier *= -1.0 + # + # term 2y + term2y = term2_y**2 + term2yFourier = self.fourier(data=term2y) + term2yFourier *= self.ly**2 + # minus sign to correct the artificial i**2 + term2yFourier *= -1.0 + # + # term 2xy + term2xy = 2.0 * term2_x * term2_y + term2xyFourier = self.fourier(data=term2xy) + term2xyFourier *= self.lx * self.ly + # minus sign to correct the artificial i**2 + term2xyFourier *= -1.0 + + if test: + self.plotFourier(term2xFourier) + self.plotFourier(term2yFourier) + self.plotFourier(term2xyFourier) + self.plotFourier(term2xFourier + term2yFourier + term2xyFourier) + + # add all terms + resultFourier = term1xFourier + term1yFourier + term1xyFourier + resultFourier += term2xFourier + term2yFourier + term2xyFourier + + # cut off the high ells from phi normalization map + f = lambda ell: (ell <= 2.0 * lMax) + resultFourier = self.filterFourierIsotropic( + f, dataFourier=resultFourier, test=False + ) + + # invert + resultFourier = 1.0 / resultFourier + resultFourier[np.where(np.isfinite(resultFourier) == False)] = 0.0 + + if test: + plt.loglog(self.ell.flatten(), resultFourier.flatten(), "b.") + plt.show() + + return resultFourier + + # Caching boiler plate + # if no caching is desired, just compute + if cache is None: + resultFourier = doCalculation() + # if caching is desired + else: + # if first call with caching, set up the cache dictionary + if not hasattr(self.computeQuadEstPhiNormalizationFFT.__func__, "cache"): + self.computeQuadEstPhiNormalizationFFT.__func__.cache = {} + # if the calculation has been done before + if self.computeQuadEstPhiNormalizationFFT.cache.has_key(cache): + resultFourier = self.computeQuadEstPhiNormalizationFFT.cache[ + cache + ].copy() + # if this calculation was not done before + else: + resultFourier = doCalculation() + self.computeQuadEstPhiNormalizationFFT.cache[ + cache + ] = resultFourier.copy() + + return resultFourier + + def forecastN0Kappa( + self, fC0, fCtot, fCfg=None, lMin=1.0, lMax=1.0e5, test=False, cache=None + ): + """Interpolates the result for N_L^kappa = f(L), + to be used for forecasts on lensing reconstruction + """ + # print "computing the reconstruction noise" + # Standard reconstruction noise + if fCfg is None: + n0Phi = self.computeQuadEstPhiNormalizationFFT( + fC0, fCtot, lMin=lMin, lMax=lMax, test=test, cache=cache + ) + # Gaussian noise contribution from the foregrounds only + else: + n0Phi = self.computeQuadEstPhiNormalizationFgNoiseFFT( + fC0, fCtot, fCfg, lMin=lMin, lMax=lMax, test=test + ) + # keep only the real part (the imag. part should be zero, but is tiny in practice) + n0Phi = np.real(n0Phi) + # remove the nans + n0Phi = np.nan_to_num(n0Phi) + # make sure every value is positive + n0Phi = np.abs(n0Phi) + + # convert from phi to kappa + n0Kappa = 0.25 * self.ell**4 * n0Phi + + # fix the issue of the wrong ell=0 value + # replace it by the value lx=0, ly=fundamental + n0Kappa[0, 0] = n0Kappa[0, 1] + + # interpolate + where = (self.ell.flatten() > 0.0) * (self.ell.flatten() < 2.0 * lMax) + L = self.ell.flatten()[where] + N = n0Kappa.flatten()[where] + lnfln = interp1d( + np.log(L), np.log(N), kind="linear", bounds_error=False, fill_value=np.inf + ) + f = lambda ell: np.exp(lnfln(np.log(ell))) + return f + + def computeQuadEstKappaNorm( + self, + fC0, + fCtot, + lMin=1.0, + lMax=1.0e5, + dataFourier=None, + dataFourier2=None, + path=None, + test=False, + cache=None, + ): + """Returns the normalized quadratic estimator for kappa in Fourier space, + and saves it to file if needed. + """ + # non-normalized QE for phi + resultFourier = self.quadEstPhiNonNorm( + fC0, + fCtot, + lMin=lMin, + lMax=lMax, + dataFourier=dataFourier, + dataFourier2=dataFourier2, + test=test, + ) + # convert from phi to kappa + resultFourier = self.kappaFromPhi(resultFourier) + # compute normalization + normalizationFourier = self.computeQuadEstPhiNormalizationFFT( + fC0, fCtot, lMin=lMin, lMax=lMax, test=test, cache=cache + ) + # normalized (not mean field-subtracted) QE for kappa + resultFourier *= normalizationFourier + # save to file if needed + if path is not None: + self.saveDataFourier(resultFourier, path) + return resultFourier + + ############################################################################### + # Gaussian N0 term from foregrounds in the CMB map + + def computeQuadEstPhiNormalizationFgNoiseFFT( + self, fC0, fCtot, fCfg, lMin=1.0, lMax=1.0e5, test=False, cache=None + ): + """Computes the N^{0 fg phiphi}_L due to foreground power in the temperature map + N^{0 fg phiphi}_L = N^{0 phiphi}_L^2 * integral. + """ + # inverse-var weighted map + def f(ell): + if (ell < lMin) or (ell > lMax): + return 0.0 + result = 1.0 / fCtot(ell) + result *= fCfg(ell) / fCtot(ell) + if not np.isfinite(result): + result = 0.0 + return result + + iVarFourier = np.array(map(f, self.ell.flatten())) + iVarFourier = iVarFourier.reshape(self.ell.shape) + iVar = self.inverseFourier(dataFourier=iVarFourier) + + # C map + def f(ell): + if (ell < lMin) or (ell > lMax): + return 0.0 + result = fC0(ell) ** 2 / fCtot(ell) + result *= fCfg(ell) / fCtot(ell) + if not np.isfinite(result): + result = 0.0 + return result + + CFourier = np.array(map(f, self.ell.flatten())) + CFourier = CFourier.reshape(self.ell.shape) + + # term 1x + term1x = self.inverseFourier(dataFourier=self.lx**2 * CFourier) + term1x *= iVar + term1xFourier = self.fourier(data=term1x) + term1xFourier *= self.lx**2 + # + # term 1y + term1y = self.inverseFourier(dataFourier=self.ly**2 * CFourier) + term1y *= iVar + term1yFourier = self.fourier(data=term1y) + term1yFourier *= self.ly**2 + # + # term 1xy + term1xy = self.inverseFourier(dataFourier=2.0 * self.lx * self.ly * CFourier) + term1xy *= iVar + term1xyFourier = self.fourier(data=term1xy) + term1xyFourier *= self.lx * self.ly + + if test: + self.plotFourier(term1xFourier) + self.plotFourier(term1yFourier) + self.plotFourier(term1xyFourier) + self.plotFourier(term1xFourier + term1yFourier + term1xyFourier) + + # WF map + def f(ell): + if (ell < lMin) or (ell > lMax): + return 0.0 + result = fC0(ell) / fCtot(ell) + result *= fCfg(ell) / fCtot(ell) + # artificial factor of i such that f(-l) = f(l)*, + # such that f(x) is real + result *= 1.0j + if not np.isfinite(result): + result = 0.0 + return result + + WFFourier = np.array(map(f, self.ell.flatten())) + WFFourier = WFFourier.reshape(self.ell.shape) + + # term 2 + term2_x = self.inverseFourier(dataFourier=self.lx * WFFourier) + term2_y = self.inverseFourier(dataFourier=self.ly * WFFourier) + # + # term 2x + term2x = term2_x**2 + term2xFourier = self.fourier(data=term2x) + term2xFourier *= self.lx**2 + # minus sign to correct the artificial i**2 + term2xFourier *= -1.0 + # + # term 2y + term2y = term2_y**2 + term2yFourier = self.fourier(data=term2y) + term2yFourier *= self.ly**2 + # minus sign to correct the artificial i**2 + term2yFourier *= -1.0 + # + # term 2xy + term2xy = 2.0 * term2_x * term2_y + term2xyFourier = self.fourier(data=term2xy) + term2xyFourier *= self.lx * self.ly + # minus sign to correct the artificial i**2 + term2xyFourier *= -1.0 + + if test: + self.plotFourier(term2xFourier) + self.plotFourier(term2yFourier) + self.plotFourier(term2xyFourier) + self.plotFourier(term2xFourier + term2yFourier + term2xyFourier) + + # add all terms + result = term1xFourier + term1yFourier + term1xyFourier + result += term2xFourier + term2yFourier + term2xyFourier + + # cut off the high ells from phi normalization map + f = lambda ell: (ell <= 2.0 * lMax) + result = self.filterFourierIsotropic(f, dataFourier=result, test=False) + + # normalize by the standard N0 squared + result *= ( + self.computeQuadEstPhiNormalizationFFT( + fC0, fCtot, lMin=lMin, lMax=lMax, test=False, cache=cache + ) + ** 2 + ) + + # clean up + result[np.where(np.isfinite(result) == False)] = 0.0 + + if test: + plt.loglog(self.ell.flatten(), result.flatten(), "b.") + plt.show() + + return result + + ############################################################################### + # Lensing multiplicative bias from lensed foregrounds + + def computeMultBiasLensedForegrounds( + self, fC0, fCtot, fCfgBias, lMin=1.0, lMax=1.0e5, test=False, cache=None + ): + """Multiplicative bias to CMB lensing + from lensed foregrounds in the CMB map. + fCfgBias: unlensed foreground power spectrum. + Such that: + = kappa_CMB + (multiplicative bias) * kappa_foreground + noise. + """ + # inverse-var weighted map + def f(ell): + if (ell < lMin) or (ell > lMax): + return 0.0 + result = 1.0 / fCtot(ell) + if not np.isfinite(result): + result = 0.0 + return result + + iVarFourier = np.array(map(f, self.ell.flatten())) + iVarFourier = iVarFourier.reshape(self.ell.shape) + iVar = self.inverseFourier(dataFourier=iVarFourier) + + # C map + def f(ell): + if (ell < lMin) or (ell > lMax): + return 0.0 + result = fC0(ell) * fCfgBias(ell) / fCtot(ell) + if not np.isfinite(result): + result = 0.0 + return result + + CFourier = np.array(map(f, self.ell.flatten())) + CFourier = CFourier.reshape(self.ell.shape) + + # term 1x + term1x = self.inverseFourier(dataFourier=self.lx**2 * CFourier) + term1x *= iVar + term1xFourier = self.fourier(data=term1x) + term1xFourier *= self.lx**2 + # + # term 1y + term1y = self.inverseFourier(dataFourier=self.ly**2 * CFourier) + term1y *= iVar + term1yFourier = self.fourier(data=term1y) + term1yFourier *= self.ly**2 + # + # term 1xy + term1xy = self.inverseFourier(dataFourier=2.0 * self.lx * self.ly * CFourier) + term1xy *= iVar + term1xyFourier = self.fourier(data=term1xy) + term1xyFourier *= self.lx * self.ly + + if test: + self.plotFourier(term1xFourier) + self.plotFourier(term1yFourier) + self.plotFourier(term1xyFourier) + self.plotFourier(term1xFourier + term1yFourier + term1xyFourier) + + # WF map + def f(ell): + if (ell < lMin) or (ell > lMax): + return 0.0 + result = fC0(ell) / fCtot(ell) + # artificial factor of i such that f(-l) = f(l)*, + # such that f(x) is real + result *= 1.0j + if not np.isfinite(result): + result = 0.0 + return result + + WFFourier = np.array(map(f, self.ell.flatten())) + WFFourier = WFFourier.reshape(self.ell.shape) + # term 2 + term2_x = self.inverseFourier(dataFourier=self.lx * WFFourier) + term2_y = self.inverseFourier(dataFourier=self.ly * WFFourier) + + # WF map, for foreground + def f(ell): + if (ell < lMin) or (ell > lMax): + return 0.0 + result = fCfgBias(ell) / fCtot(ell) + # artificial factor of i such that f(-l) = f(l)*, + # such that f(x) is real + result *= 1.0j + if not np.isfinite(result): + result = 0.0 + return result + + fgWFFourier = np.array(map(f, self.ell.flatten())) + fgWFFourier = fgWFFourier.reshape(self.ell.shape) + # term 2 + term2_x_fg = self.inverseFourier(dataFourier=self.lx * fgWFFourier) + term2_y_fg = self.inverseFourier(dataFourier=self.ly * fgWFFourier) + + # term 2x + term2x = term2_x * term2_x_fg + term2xFourier = self.fourier(data=term2x) + term2xFourier *= self.lx**2 + # minus sign to correct the artificial i**2 + term2xFourier *= -1.0 + # + # term 2y + term2y = term2_y * term2_y_fg + term2yFourier = self.fourier(data=term2y) + term2yFourier *= self.ly**2 + # minus sign to correct the artificial i**2 + term2yFourier *= -1.0 + # + # term 2xy + term2xy = term2_x * term2_y_fg + term2xy += term2_x_fg * term2_y + term2xyFourier = self.fourier(data=term2xy) + term2xyFourier *= self.lx * self.ly + # minus sign to correct the artificial i**2 + term2xyFourier *= -1.0 + + if test: + self.plotFourier(term2xFourier) + self.plotFourier(term2yFourier) + self.plotFourier(term2xyFourier) + self.plotFourier(term2xFourier + term2yFourier + term2xyFourier) + + # add all terms + result = term1xFourier + term1yFourier + term1xyFourier + result += term2xFourier + term2yFourier + term2xyFourier + + # cut off the high ells from phi normalization map + f = lambda ell: (ell <= 2.0 * lMax) + result = self.filterFourierIsotropic(f, dataFourier=result, test=False) + + # normalize by the standard N0 + result *= self.computeQuadEstPhiNormalizationFFT( + fC0, fCtot, lMin=lMin, lMax=lMax, test=False, cache=cache + ) + result[np.where(np.isfinite(result) == False)] = 0.0 + + if test: + plt.loglog(self.ell.flatten(), np.abs(result.flatten()), "k.") + plt.loglog(self.ell.flatten(), np.abs(np.real(result.flatten())), "b.") + plt.loglog(self.ell.flatten(), np.imag(result.flatten()), "g.") + plt.show() + + return result + + def forecastMultBiasLensedForegrounds( + self, fC0, fCtot, fCfgBias, lMin=1.0, lMax=1.0e5, test=False + ): + """Interpolates the multiplicative bias to CMB lensing + due to lensed foregrounds in the map + """ + # print "computing the multiplicative bias" + result = self.computeMultBiasLensedForegrounds( + fC0, fCtot, fCfgBias, lMin=lMin, lMax=lMax, test=test + ) + # keep only the real part (the imag. part should be zero, but is tiny in practice) + result = np.real(result) + # remove the nans + result = np.nan_to_num(result) + + # fix the issue of the wrong ell=0 value + # replace it by the value lx=0, ly=fundamental + result[0, 0] = result[0, 1] + + # interpolate, preserving the sign + lnfln = interp1d( + np.log(self.ell.flatten()), + np.log(np.abs(result).flatten()), + kind="linear", + bounds_error=False, + fill_value=0.0, + ) + signln = interp1d( + np.log(self.ell.flatten()), + np.sign(result.flatten()), + kind="linear", + bounds_error=False, + fill_value=0.0, + ) + f = lambda ell: np.exp(lnfln(np.log(ell))) * signln(np.log(ell)) + return f + + ############################################################################### + # Converting a map trispectrum into a non-Gaussian N_L^kappa + + def computeConversionTrispecToNoisePhi( + self, fC0, fCtot, lMin=1.0, lMax=1.0e5, test=False, cache=None + ): + """Computes conversion factor, such that: + N_L^phi = Gaussian + (conversion factor) * (white trispectrum of temperature map) + The conversion factor has units of 1/C_l^2 (omitting radians) + The integral to compute is very similar to that for the + Gaussian N_L^phi. + """ + # inverse-var weighted map + def f(ell): + if (ell < lMin) or (ell > lMax): + return 0.0 + result = 1.0 / fCtot(ell) + if not np.isfinite(result): + result = 0.0 + return result + + iVarFourier = np.array(map(f, self.ell.flatten())) + iVarFourier = iVarFourier.reshape(self.ell.shape) + iVar = self.inverseFourier(dataFourier=iVarFourier) + + # C map + def f(ell): + if (ell < lMin) or (ell > lMax): + return 0.0 + result = fC0(ell) / fCtot(ell) + # artificial factor of i such that f(-l) = f(l)*, + # such that f(x) is real + result *= 1.0j + if not np.isfinite(result): + result = 0.0 + return result + + CFourier = np.array(map(f, self.ell.flatten())) + CFourier = CFourier.reshape(self.ell.shape) + + # term x + termx = self.inverseFourier(dataFourier=self.lx * CFourier) + termx *= iVar + termxFourier = self.fourier(data=termx) + termxFourier *= self.lx + # correct for the artificial factor of i + termxFourier *= -1.0j + # + # term y + termy = self.inverseFourier(dataFourier=self.ly * CFourier) + termy *= iVar + termyFourier = self.fourier(data=termy) + termyFourier *= self.ly + # correct for the artificial factor of i + termyFourier *= -1.0j + + if test: + self.plotFourier(termxFourier) + self.plotFourier(termyFourier) + self.plotFourier(termxFourier + termyFourier) + + # add all terms + resultFourier = termxFourier + termyFourier + + # cut off the high ells from phi normalization map + f = lambda ell: (ell <= 2.0 * lMax) + resultFourier = self.filterFourierIsotropic( + f, dataFourier=resultFourier, test=False + ) + + if test: + plt.loglog(self.ell.flatten(), resultFourier.flatten(), "b.") + plt.show() + + # take modulus squared + resultFourier = np.abs(resultFourier) ** 2 + + # normalize by the squared Gaussian N0 + n0GFourier = self.computeQuadEstPhiNormalizationFFT( + fC0, fCtot, lMin=lMin, lMax=lMax, test=False, cache=cache + ) + n0GFourier = np.abs(n0GFourier) ** 2 + resultFourier *= n0GFourier + + if test: + plt.loglog(self.ell.flatten(), resultFourier.flatten(), "b.") + plt.show() + + return resultFourier + + def computeConversionTrispecToNoiseKappa( + self, fC0, fCtot, lMin=1.0, lMax=1.0e5, test=False + ): + # print "computing the conversion factor" + convPhi = self.computeConversionTrispecToNoisePhi( + fC0, fCtot, lMin=lMin, lMax=lMax, test=test + ) + # keep only the real part (the imag. part should be zero, but is tiny in practice) + convPhi = np.real(convPhi) + # remove the nans + convPhi = np.nan_to_num(convPhi) + # make sure every value is positive + convPhi = np.abs(convPhi) + + # convert from phi to kappa + convKappa = 0.25 * self.ell**4 * convPhi + + # fix the issue of the wrong ell=0 value + # replace it by the value lx=0, ly=fundamental + convKappa[0, 0] = convKappa[0, 1] + + # interpolate + f = interp1d( + self.ell.flatten(), + convKappa.flatten(), + kind="linear", + bounds_error=False, + fill_value=0.0, + ) + return f + + ############################################################################### + # Non-normalized mean field for quadratic estimator + + def computeQuadEstPhiMeanFieldNonNormRand( + self, + fC0, + fCtot, + lMin=1.0, + lMax=1.0e5, + nRand=1, + path=None, + dataFourier=None, + test=False, + ): + """Mean field for quadratic estimator from randoms""" + if dataFourier is None: + dataFourier = self.dataFourier.copy() + + result = np.zeros_like(self.dataFourier) + + for iRand in range(nRand): + # # randomize the phases + # randomFourier = self.randomizePhases(dataFourier=dataFourier) + if path is not None: + # print "read mock", iRand, "from", path+str(iRand)+".fits" + randomFourier = self.loadDataFourier(path + str(iRand) + ".fits") + else: + # print "generate mock", iRand + randomFourier = self.genGRF(fCtot, test=test) + + # get the non-normalized quad. est. for this map + randomQuadFourier = self.quadEstPhiNonNorm( + fC0, fCtot, lMin=lMin, lMax=lMax, dataFourier=randomFourier, test=test + ) + result += randomQuadFourier + result /= nRand + + if test: + # print "showing normalization map" + self.plotFourier(dataFourier=result) + # print "checking power spectrum" + self.powerSpectrum(dataFourier=result) + + return result + + ############################################################################### + # Smarter estimator for C_L^kk, smarter than the auto of the kappa QE: + # avoids the N0 bias, and the secondary foreground biases. + + def computeQuadEstPhiInverseDataNormalizationFFT( + self, fC0, fCtot, lMin=1.0, lMax=1.0e5, dataFourier=None, test=False + ): + """Analogous to the standard QE normalization 1/N_L, + but from the data. + Used to estimate C_L^phiphi auto from cross only, + and to reduce secondary foreground bias. + """ + if dataFourier is None: + dataFourier = self.dataFourier.copy() + + # inverse-var weighted map + def f(ell): + if (ell < lMin) or (ell > lMax): + return 0.0 + result = 1.0 / fCtot(ell) ** 2 + if not np.isfinite(result): + result = 0.0 + return result + + iVarFourier = np.array(map(f, self.ell.flatten())) + iVarFourier = iVarFourier.reshape(self.ell.shape) + # multiply by data squared modulus + iVarFourier *= np.abs(dataFourier) ** 2 + iVar = self.inverseFourier(dataFourier=iVarFourier) + + # C map + def f(ell): + if (ell < lMin) or (ell > lMax): + return 0.0 + result = fC0(ell) ** 2 / fCtot(ell) ** 2 + if not np.isfinite(result): + result = 0.0 + return result + + CFourier = np.array(map(f, self.ell.flatten())) + CFourier = CFourier.reshape(self.ell.shape) + # multiply by data squared modulus + CFourier *= np.abs(dataFourier) ** 2 + + # term 1x + term1x = self.inverseFourier(dataFourier=self.lx**2 * CFourier) + term1x *= iVar + term1xFourier = self.fourier(data=term1x) + term1xFourier *= self.lx**2 + # + # term 1y + term1y = self.inverseFourier(dataFourier=self.ly**2 * CFourier) + term1y *= iVar + term1yFourier = self.fourier(data=term1y) + term1yFourier *= self.ly**2 + # + # term 1xy + term1xy = self.inverseFourier(dataFourier=2.0 * self.lx * self.ly * CFourier) + term1xy *= iVar + term1xyFourier = self.fourier(data=term1xy) + term1xyFourier *= self.lx * self.ly + + if test: + self.plotFourier(term1xFourier) + self.plotFourier(term1yFourier) + self.plotFourier(term1xyFourier) + self.plotFourier(term1xFourier + term1yFourier + term1xyFourier) + + # WF map + def f(ell): + if (ell < lMin) or (ell > lMax): + return 0.0 + result = fC0(ell) / fCtot(ell) ** 2 + # artificial factor of i such that f(-l) = f(l)*, + # such that f(x) is real + result *= 1.0j + if not np.isfinite(result): + result = 0.0 + return result + + WFFourier = np.array(map(f, self.ell.flatten())) + WFFourier = WFFourier.reshape(self.ell.shape) + # multiply by data squared modulus + WFFourier *= np.abs(dataFourier) ** 2 + + # term 2 + term2_x = self.inverseFourier(dataFourier=self.lx * WFFourier) + term2_y = self.inverseFourier(dataFourier=self.ly * WFFourier) + # + # term 2x + term2x = term2_x**2 + term2xFourier = self.fourier(data=term2x) + term2xFourier *= self.lx**2 + # minus sign to correct the artificial i**2 + term2xFourier *= -1.0 + # + # term 2y + term2y = term2_y**2 + term2yFourier = self.fourier(data=term2y) + term2yFourier *= self.ly**2 + # minus sign to correct the artificial i**2 + term2yFourier *= -1.0 + # + # term 2xy + term2xy = 2.0 * term2_x * term2_y + term2xyFourier = self.fourier(data=term2xy) + term2xyFourier *= self.lx * self.ly + # minus sign to correct the artificial i**2 + term2xyFourier *= -1.0 + + if test: + self.plotFourier(term2xFourier) + self.plotFourier(term2yFourier) + self.plotFourier(term2xyFourier) + self.plotFourier(term2xFourier + term2yFourier + term2xyFourier) + + # add all terms + resultFourier = term1xFourier + term1yFourier + term1xyFourier + resultFourier += term2xFourier + term2yFourier + term2xyFourier + + # cut off the high ells from phi normalization map + f = lambda ell: (ell <= 2.0 * lMax) + resultFourier = self.filterFourierIsotropic( + f, dataFourier=resultFourier, test=False + ) + + # invert + # resultFourier = 1./resultFourier + resultFourier[np.where(np.isfinite(resultFourier) == False)] = 0.0 + + if test: + plt.loglog(self.ell.flatten(), resultFourier.flatten(), "b.") + plt.show() + + return resultFourier + + def computeQuadEstKappaAutoCorrectionMap( + self, + fC0, + fCtot, + lMin=1.0, + lMax=1.0e5, + dataFourier=None, + path=None, + test=False, + cache=None, + ): + # non-normalized phi inverse normalization map + resultFourier = self.computeQuadEstPhiInverseDataNormalizationFFT( + fC0, fCtot, lMin=lMin, lMax=lMax, dataFourier=dataFourier, test=test + ) + # convert from phi to kappa + resultFourier = self.kappaFromPhi(resultFourier) + # do it again, since this is effectively a "power spectrum map" + resultFourier = self.kappaFromPhi(resultFourier) + # compute normalization + normalizationFourier = self.computeQuadEstPhiNormalizationFFT( + fC0, fCtot, lMin=lMin, lMax=lMax, test=test, cache=cache + ) + # normalized correction for QE kappa auto-spectrum correction map + resultFourier *= normalizationFourier**2 + #!!!!!!!!! weird factor needed. I haven't figured out why + resultFourier /= self.sizeX * self.sizeY + # take square root, so that all you have to do is to take the power spectrum + resultFourier = np.sqrt(np.real(resultFourier)) + # save to file if needed + if path is not None: + self.saveDataFourier(resultFourier, path) + return resultFourier + + ############################################################################### + ############################################################################### + # Dilation-only estimator + + def quadEstPhiDilationNonNorm( + self, + fC0, + fCtot, + lMin=5.0e2, + lMax=3000.0, + dataFourier=None, + dataFourier2=None, + test=False, + ): + """Non-normalized quadratic estimator for phi from dilation only + fC0: ulensed power spectrum + fCtot: lensed power spectrum + noise + ell cuts are performed to remain in the regime L_phi < l_T + """ + if dataFourier is None: + dataFourier = self.dataFourier.copy() + if dataFourier2 is None: + dataFourier2 = dataFourier.copy() + + # cut off high ells + f = lambda ell: 1.0 * (ell >= lMin) * (ell <= lMax) + FDataFourier = self.filterFourierIsotropic( + f, dataFourier=dataFourier, test=test + ) + FData = self.inverseFourier(FDataFourier) + + # weight function for dilation + def fdLnl2C0dLnl(ell): + e = 0.01 + lup = ell * (1.0 + e) + ldown = ell * (1.0 - e) + result = lup**2 * fC0(lup) + result /= ldown**2 * fC0(ldown) + result = np.log(result) / (2.0 * e) + return result + + if test: + # print "testing derivative" + F = np.array(map(fdLnl2C0dLnl, self.ell.flatten())) + plt.semilogx(self.ell.flatten(), F, "b.") + plt.show() + + # sort of Dilation Wiener-filter + def f(ell): + # cut off the high ells from input map + if (ell < lMin) or (ell > lMax): + return 0.0 + result = fC0(ell) / fCtot(ell) ** 2 + result *= fdLnl2C0dLnl(ell) # for isotropic dilation + result /= 0.5 + if not np.isfinite(result): + result = 0.0 + return result + + WFDataFourier = self.filterFourierIsotropic( + f, dataFourier=dataFourier2, test=test + ) + WFData = self.inverseFourier(WFDataFourier) + if test: + # print "showing the WF map" + self.plot(data=WFData) + # print "checking the power spectrum of this map" + theory = lambda ell: f(ell) ** 2 * fC0(ell) + self.powerSpectrum(dataFourier=WFDataFourier, theory=[theory], plot=True) + + # product in real space + product = FData * WFData + productFourier = self.fourier(product) + + # get phi from kappa + productFourier *= -2.0 / self.ell**2 + productFourier[np.where(np.isfinite(productFourier) == False)] = 0.0 + if test: + # print "checking the power spectrum of phi map" + self.powerSpectrum(dataFourier=productFourier, plot=True) + + if test: + "Show real space map" + product = self.inverseFourier(productFourier) + self.plot(product) + + return productFourier + + def computeQuadEstPhiDilationNormalizationFFT( + self, fC0, fCtot, fC0wg=None, lMin=5.0e2, lMax=3.0e3, test=False + ): + """Multiplicative normalization for phi estimator from dilation only, + computed with FFT. + This normalization does not correct for the estimator's multiplicative bias. + C0wg: a hypothetical wiggle-only or no-wiggle unlensed CMB power spectrum. + """ + if fC0wg is None: + fC0wg = fC0 + + def fdLnl2C0dLnl(ell): + e = 0.01 + lup = ell * (1.0 + e) + ldown = ell * (1.0 - e) + result = lup**2 * fC0wg(lup) + result /= ldown**2 * fC0wg(ldown) + result = np.log(result) / (2.0 * e) + return result + + def f(ell): + # cut off the high ells from input map + if (ell < lMin) or (ell > lMax): + return 0.0 + result = fC0wg(ell) / fCtot(ell) + result *= fdLnl2C0dLnl(ell) # for isotropic dilation + result = result**2 + result /= 0.5 + if not np.isfinite(result): + result = 0.0 + return result + + dataFourier = self.filterFourierIsotropic( + f, dataFourier=np.ones_like(self.ell), test=test + ) + + # compute the integral + integral = np.sum(dataFourier) + integral /= self.sizeX * self.sizeY + if test: + integralFFT = self.inverseFourier(dataFourier) + integralFFT = ( + integralFFT[0, 0] / 2.0 + ) # I think the factor 2 is about the half/full Fourier plane + # print "Integral from sum=", integral + # print "Same integral from FFT=", integralFFT + + # fill a Fourier map with this value + resultFourier = np.ones_like(dataFourier) + resultFourier /= integral + + return resultFourier + + def computeQuadEstPhiDilationNormalizationCorrectedFFT( + self, fC0, fCtot, fC0wg=None, lMin=5.0e2, lMax=3.0e3, test=False, cache=None + ): + """Multiplicative normalization for phi estimator from dilation only, + computed with FFT. + This normalization corrects for the multiplicative bias in the estimator. + C0wg: a hypothetical wiggle-only or no-wiggle unlensed CMB power spectrum. + """ + + if fC0wg is None: + fC0wg = fC0 + + def doCalculation(): + def fdLnl2C0dLnl(ell): + e = 0.01 + lup = ell * (1.0 + e) + ldown = ell * (1.0 - e) + result = lup**2 * fC0wg(lup) + result /= ldown**2 * fC0wg(ldown) + result = np.log(result) / (2.0 * e) + return result + + def fDilation(ell): + # cut off the high ells from input map + if (ell < lMin) or (ell > lMax): + return 0.0 + result = fC0wg(ell) / fCtot(ell) ** 2 + result *= fdLnl2C0dLnl(ell) # for isotropic dilation + result /= 0.5 + if not np.isfinite(result): + result = 0.0 + return result + + # generate dilation map + dilationFourier = self.filterFourierIsotropic( + fDilation, dataFourier=np.ones_like(self.ell), test=test + ) + dilation = self.inverseFourier(dilationFourier) + + # generate gradient C0 map + f = lambda ell: fC0(ell) * (ell >= lMin) * (ell <= lMax) + c0Fourier = self.filterFourierIsotropic( + f, dataFourier=np.ones_like(self.ell), test=test + ) + # the factor i in the gradient makes the Fourier function Hermitian + gradXFourier, gradYFourier = self.computeGradient(dataFourier=c0Fourier) + gradX = self.inverseFourier( + gradXFourier + ) # extra factor of i will be cancelled later + gradY = self.inverseFourier( + gradYFourier + ) # extra factor of i will be cancelled later + + # generate ell limit map + f = lambda ell: (ell >= lMin) * (ell <= lMax) + ellLimitsFourier = self.filterFourierIsotropic( + f, dataFourier=np.ones_like(self.ell), test=test + ) + ellLimits = self.inverseFourier(ellLimitsFourier) + + # First, the asymmetric term + # term1x + term1XFourier = self.fourier(gradX * dilation) + term1XFourier *= ( + 2.0 * self.lx / self.ell**2 / 1.0j + ) # factor of i to cancel the one in the gradient + term1XFourier[np.where(np.isfinite(term1XFourier) == False)] = 0.0 + # term1y + term1YFourier = self.fourier(gradY * dilation) + term1YFourier *= ( + 2.0 * self.ly / self.ell**2 / 1.0j + ) # factor of i to cancel the one in the gradient + term1YFourier[np.where(np.isfinite(term1YFourier) == False)] = 0.0 + # sum + term1Fourier = term1XFourier + term1YFourier + if test: + # print "showing term1XFourier" + self.plotFourier(term1XFourier) + # print "showing term1YFourier" + self.plotFourier(term1YFourier) + # print "showing term1Fourier" + self.plotFourier(term1Fourier) + + # Second, the symmetric term + # term2x + term2X = self.inverseFourier(gradXFourier * dilationFourier) + term2XFourier = self.fourier(term2X * ellLimits) + term2XFourier *= ( + 2.0 * self.lx / self.ell**2 / 1.0j + ) # factor of i to cancel the one in the gradient + term2XFourier[np.where(np.isfinite(term2XFourier) == False)] = 0.0 + # term2y + term2Y = self.inverseFourier(gradYFourier * dilationFourier) + term2YFourier = self.fourier(term2Y * ellLimits) + term2YFourier *= ( + 2.0 * self.ly / self.ell**2 / 1.0j + ) # factor of i to cancel the one in the gradient + term2YFourier[np.where(np.isfinite(term2XFourier) == False)] = 0.0 + # sum + term2Fourier = term2XFourier + term2YFourier + if test: + # print "showing term2XFourier" + self.plotFourier(term2XFourier) + # print "showing term2YFourier" + self.plotFourier(term2YFourier) + # print "showing term2Fourier" + self.plotFourier(term2Fourier) + + # sum and invert + resultFourier = 1.0 / (term1Fourier + term2Fourier) + resultFourier[np.where(np.isfinite(resultFourier) == False)] = 0.0 + # remove L > 2 lMax + f = lambda ell: (ell <= 2.0 * lMax) + resultFourier = self.filterFourierIsotropic( + f, dataFourier=resultFourier, test=test + ) + if test: + # print "showing sum" + self.plotFourier(resultFourier) + + return resultFourier + + # Caching boiler plate + # if no caching is desired, just compute + if cache is None: + resultFourier = doCalculation() + # if caching is desired + else: + # if first call with caching, set up the cache dictionary + if not hasattr( + self.computeQuadEstPhiDilationNormalizationCorrectedFFT.__func__, + "cache", + ): + self.computeQuadEstPhiDilationNormalizationCorrectedFFT.__func__.cache = ( + {} + ) + # if the calculation has been done before + if self.computeQuadEstPhiDilationNormalizationCorrectedFFT.cache.has_key( + cache + ): + resultFourier = ( + self.computeQuadEstPhiDilationNormalizationCorrectedFFT.cache[ + cache + ].copy() + ) + # if this calculation was not done before + else: + resultFourier = doCalculation() + self.computeQuadEstPhiDilationNormalizationCorrectedFFT.cache[ + cache + ] = resultFourier.copy() + + return resultFourier + + def computeQuadEstKappaDilationNormCorr( + self, + fC0, + fCtot, + fC0wg=None, + lMin=1.0, + lMax=1.0e5, + dataFourier=None, + dataFourier2=None, + path=None, + corr=True, + test=False, + cache=None, + ): + # non-normalized QE for phi + resultFourier = self.quadEstPhiDilationNonNorm( + fC0, + fCtot, + lMin=lMin, + lMax=lMax, + dataFourier=dataFourier, + dataFourier2=dataFourier2, + test=test, + ) + # convert from phi to kappa + resultFourier = self.kappaFromPhi(resultFourier) + # compute normalization (no mean field subtraction) + if corr: + resultFourier *= self.computeQuadEstPhiDilationNormalizationCorrectedFFT( + fC0, fCtot, fC0wg=fC0wg, lMin=lMin, lMax=lMax, test=test, cache=cache + ) + else: + resultFourier *= self.computeQuadEstPhiDilationNormalizationFFT( + fC0, fCtot, fC0wg=fC0wg, lMin=lMin, lMax=lMax, test=test + ) + # save to file if needed + if path is not None: + self.saveDataFourier(resultFourier, path) + return resultFourier + + def computeQuadEstKappaDilationNoiseFFT( + self, + fC0, + fCtot, + fCfg=None, + fC0wg=None, + lMin=1.0, + lMax=1.0e5, + corr=True, + test=False, + cache=None, + ): + """Computes the lensing noise power spectrum N_L^kappa + for the dilation estimator. + fCfg: gives the N0 due only to the foreground component in the map. + fC0wg: to replace consistently C0 by a wiggles-only power spectrum. + """ + if fCfg is None: + fCfg = fCtot + if fC0wg is None: + fC0wg = fC0 + + def fdLnl2C0dLnl(ell): + e = 0.01 + lup = ell * (1.0 + e) + ldown = ell * (1.0 - e) + result = lup**2 * fC0wg(lup) + result /= ldown**2 * fC0wg(ldown) + result = np.log(result) / (2.0 * e) + return result + + def g(ell): + # cut off the high ells from input map + if (ell < lMin) or (ell > lMax): + return 0.0 + result = fC0wg(ell) / fCtot(ell) ** 2 + result *= fdLnl2C0dLnl(ell) # for isotropic dilation + result /= 0.5 + if not np.isfinite(result): + result = 0.0 + return result + + # First, the symmetric term + f = lambda ell: g(ell) * fCfg(ell) + term1Fourier = self.filterFourierIsotropic( + f, dataFourier=np.ones_like(self.ell), test=test + ) + term1 = self.inverseFourier(term1Fourier) + term1Fourier = self.fourier(term1**2) + if test: + self.plotFourier(term1Fourier) + + # Second, the asymmetric term + f = lambda ell: g(ell) ** 2 * fCfg(ell) + term2aFourier = self.filterFourierIsotropic( + f, dataFourier=np.ones_like(self.ell), test=test + ) + term2a = self.inverseFourier(term2aFourier) + # + f = lambda ell: fCfg(ell) * (ell >= lMin) * (ell <= lMax) + term2bFourier = self.filterFourierIsotropic( + f, dataFourier=np.ones_like(self.ell), test=test + ) + term2b = self.inverseFourier(term2bFourier) + # + term2Fourier = self.fourier(term2a * term2b) + if test: + self.plotFourier(term2Fourier) + + # add term1 and term2 + resultFourier = term1Fourier + term2Fourier + if test: + self.plotFourier(term1Fourier + term2Fourier) + + # compute normalization + if corr: + normalizationFourier = ( + self.computeQuadEstPhiDilationNormalizationCorrectedFFT( + fC0, + fCtot, + fC0wg=fC0wg, + lMin=lMin, + lMax=lMax, + test=test, + cache=cache, + ) + ) + else: + normalizationFourier = self.computeQuadEstPhiDilationNormalizationFFT( + fC0, fCtot, fC0wg=fC0wg, lMin=lMin, lMax=lMax, test=test + ) + + # multiply by squared normalization + resultFourier *= normalizationFourier**2 + + # remove L > 2 lMax + f = lambda ell: (ell <= 2.0 * lMax) + resultFourier = self.filterFourierIsotropic( + f, dataFourier=resultFourier, test=test + ) + + return resultFourier + + def forecastN0KappaDilation( + self, + fC0, + fCtot, + fCfg=None, + fC0wg=None, + lMin=1.0, + lMax=1.0e5, + corr=True, + test=False, + cache=None, + ): + """Interpolates the result for N_L^{kappa_dilation} = f(l), + to be used for forecasts on lensing reconstruction. + fCfg: gives the N0 due only to the foreground component in the map. + fC0wg: to replace consistently C0 by a wiggles-only power spectrum. + """ + # print "computing the reconstruction noise" + n0Kappa = self.computeQuadEstKappaDilationNoiseFFT( + fC0, + fCtot, + fCfg=fCfg, + fC0wg=fC0wg, + lMin=lMin, + lMax=lMax, + corr=corr, + test=test, + cache=cache, + ) + # keep only the real part (the imag. part should be zero, but is tiny in practice) + n0Kappa = np.real(n0Kappa) + # remove the nans + n0Kappa = np.nan_to_num(n0Kappa) + # make sure every value is positive + n0Kappa = np.abs(n0Kappa) + + # fix the issue of the wrong ell=0 value + # replace it by the value lx=0, ly=fundamental + n0Kappa[0, 0] = n0Kappa[0, 1] + + # interpolate + where = (self.ell.flatten() > 0.0) * (self.ell.flatten() < 2.0 * lMax) + L = self.ell.flatten()[where] + N = n0Kappa.flatten()[where] + lnfln = interp1d( + np.log(L), np.log(N), kind="linear", bounds_error=False, fill_value=np.inf + ) + f = lambda ell: np.exp(lnfln(np.log(ell))) + return f + + ############################################################################### + # Lensing multiplicative bias from lensed foregrounds + + def computeMultBiasLensedForegroundsDilation( + self, fC0, fCtot, fCfgBias, lMin=5.0e2, lMax=3.0e3, test=False, cache=None + ): + """Multiplicative bias to the dilation estimator + from lensed foregrounds in the CMB map. + fCfgBias: unlensed foreground power spectrum. + Such that: + = kappa_CMB + (multiplicative bias) * kappa_foreground + noise. + """ + + def fdLnl2C0dLnl(ell): + e = 0.01 + lup = ell * (1.0 + e) + ldown = ell * (1.0 - e) + result = lup**2 * fC0(lup) + result /= ldown**2 * fC0(ldown) + result = np.log(result) / (2.0 * e) + return result + + def fDilation(ell): + # cut off the high ells from input map + if (ell < lMin) or (ell > lMax): + return 0.0 + result = fC0(ell) / fCtot(ell) ** 2 + result *= fdLnl2C0dLnl(ell) # for isotropic dilation + result /= 0.5 + if not np.isfinite(result): + result = 0.0 + return result + + # generate dilation map + dilationFourier = self.filterFourierIsotropic( + fDilation, dataFourier=np.ones_like(self.ell), test=test + ) + dilation = self.inverseFourier(dilationFourier) + + # generate gradient C0 map + f = lambda ell: fCfgBias(ell) * (ell >= lMin) * (ell <= lMax) + c0Fourier = self.filterFourierIsotropic( + f, dataFourier=np.ones_like(self.ell), test=test + ) + # the factor i in the gradient makes the Fourier function Hermitian + gradXFourier, gradYFourier = self.computeGradient(dataFourier=c0Fourier) + gradX = self.inverseFourier( + gradXFourier + ) # extra factor of i will be cancelled later + gradY = self.inverseFourier( + gradYFourier + ) # extra factor of i will be cancelled later + + # generate ell limit map + f = lambda ell: (ell >= lMin) * (ell <= lMax) + ellLimitsFourier = self.filterFourierIsotropic( + f, dataFourier=np.ones_like(self.ell), test=test + ) + ellLimits = self.inverseFourier(ellLimitsFourier) + + # First, the asymmetric term + # term1x + term1XFourier = self.fourier(gradX * dilation) + term1XFourier *= ( + 2.0 * self.lx / self.ell**2 / 1.0j + ) # factor of i to cancel the one in the gradient + term1XFourier[np.where(np.isfinite(term1XFourier) == False)] = 0.0 + # term1y + term1YFourier = self.fourier(gradY * dilation) + term1YFourier *= ( + 2.0 * self.ly / self.ell**2 / 1.0j + ) # factor of i to cancel the one in the gradient + term1YFourier[np.where(np.isfinite(term1YFourier) == False)] = 0.0 + # sum + term1Fourier = term1XFourier + term1YFourier + if test: + # print "showing term1XFourier" + self.plotFourier(term1XFourier) + # print "showing term1YFourier" + self.plotFourier(term1YFourier) + # print "showing term1Fourier" + self.plotFourier(term1Fourier) + + # Second, the symmetric term + # term2x + term2X = self.inverseFourier(gradXFourier * dilationFourier) + term2XFourier = self.fourier(term2X * ellLimits) + term2XFourier *= ( + 2.0 * self.lx / self.ell**2 / 1.0j + ) # factor of i to cancel the one in the gradient + term2XFourier[np.where(np.isfinite(term2XFourier) == False)] = 0.0 + # term2y + term2Y = self.inverseFourier(gradYFourier * dilationFourier) + term2YFourier = self.fourier(term2Y * ellLimits) + term2YFourier *= ( + 2.0 * self.ly / self.ell**2 / 1.0j + ) # factor of i to cancel the one in the gradient + term2YFourier[np.where(np.isfinite(term2XFourier) == False)] = 0.0 + # sum + term2Fourier = term2XFourier + term2YFourier + if test: + # print "showing term2XFourier" + self.plotFourier(term2XFourier) + # print "showing term2YFourier" + self.plotFourier(term2YFourier) + # print "showing term2Fourier" + self.plotFourier(term2Fourier) + + # sum + resultFourier = term1Fourier + term2Fourier + resultFourier[np.where(np.isfinite(resultFourier) == False)] = 0.0 + # remove L > 2 lMax + f = lambda ell: (ell <= 2.0 * lMax) + resultFourier = self.filterFourierIsotropic( + f, dataFourier=resultFourier, test=test + ) + + # normalize with the standard normalization + resultFourier *= self.computeQuadEstPhiDilationNormalizationCorrectedFFT( + fC0, fCtot, fC0wg=None, lMin=lMin, lMax=lMax, test=False, cache=cache + ) + + resultFourier[np.where(np.isfinite(resultFourier) == False)] = 0.0 + if test: + # print "showing sum" + self.plotFourier(resultFourier) + + return resultFourier + + def forecastMultBiasLensedForegroundsDilation( + self, fC0, fCtot, fCfgBias, lMin=1.0, lMax=1.0e5, test=False + ): + """Interpolates the multiplicative bias to the dilation estimator + due to lensed foregrounds in the map + """ + # print "computing the multiplicative bias" + result = self.computeMultBiasLensedForegroundsDilation( + fC0, fCtot, fCfgBias, lMin=lMin, lMax=lMax, test=test + ) + # keep only the real part (the imag. part should be zero, but is tiny in practice) + result = np.real(result) + # remove the nans + result = np.nan_to_num(result) + + # fix the issue of the wrong ell=0 value + # replace it by the value lx=0, ly=fundamental + result[0, 0] = result[0, 1] + + # interpolate, preserving the sign + lnfln = interp1d( + np.log(self.ell.flatten()), + np.log(np.abs(result).flatten()), + kind="linear", + bounds_error=False, + fill_value=0.0, + ) + signln = interp1d( + np.log(self.ell.flatten()), + np.sign(result.flatten()), + kind="linear", + bounds_error=False, + fill_value=0.0, + ) + f = lambda ell: np.exp(lnfln(np.log(ell))) * signln(np.log(ell)) + return f + + ############################################################################### + ############################################################################### + # Shear-only estimator + + def quadEstPhiShearNonNorm( + self, + fC0, + fCtot, + lMin=5.0e2, + lMax=3000.0, + dataFourier=None, + dataFourier2=None, + test=False, + ): + """Non-normalized quadratic estimator for phi from shear only + fC0: ulensed power spectrum + fCtot: lensed power spectrum + noise + """ + if dataFourier is None: + dataFourier = self.dataFourier.copy() + if dataFourier2 is None: + dataFourier2 = dataFourier.copy() + + # cut off high ells + f = lambda ell: 1.0 * (ell >= lMin) * (ell <= lMax) + FDataFourier = self.filterFourierIsotropic( + f, dataFourier=dataFourier, test=test + ) + FData = self.inverseFourier(FDataFourier) + if test: + # print "show Fourier data" + self.plotFourier(dataFourier=FDataFourier) + + # weight function for shear + def fdLnC0dLnl(ell): + e = 0.01 + lup = ell * (1.0 + e) + ldown = ell * (1.0 - e) + result = fC0(lup) / fC0(ldown) + result = np.log(result) / (2.0 * e) + return result + + if test: + F = np.array(map(fdLnC0dLnl, self.ell.flatten())) + plt.semilogx(self.ell.flatten(), F, "b.") + plt.show() + + # sort of shear Wiener-filter + def f(ell): + # cut off the high ells from input map + if (ell < lMin) or (ell > lMax): + return 0.0 + result = fC0(ell) / fCtot(ell) ** 2 + result *= fdLnC0dLnl(ell) # for shear + result /= 0.5 + if not np.isfinite(result): + result = 0.0 + return result + + WFDataFourier = self.filterFourierIsotropic( + f, dataFourier=dataFourier2, test=test + ) + if test: + # print "showing the WF map" + self.plotFourier(dataFourier=WFDataFourier) + # print "checking the power spectrum of this map" + theory = lambda ell: f(ell) ** 2 * fC0(ell) + self.powerSpectrum(theory=[theory], dataFourier=WFDataFourier, plot=True) + + # multiplication by cos 2 theta_{L,l} + # + # term 1 + def f(lx, ly): + l2 = lx**2 + ly**2 + if l2 == 0: + return 0.0 + else: + return (lx**2 - ly**2) / l2 + + term1Fourier = self.filterFourier(f, dataFourier=WFDataFourier, test=test) + term1 = self.inverseFourier(dataFourier=term1Fourier) + term1 *= FData + term1Fourier = self.fourier(data=term1) + term1Fourier = self.filterFourier(f, dataFourier=term1Fourier, test=test) + # + + # term 2 + def f(lx, ly): + l2 = lx**2 + ly**2 + if l2 == 0: + return 0.0 + else: + return lx * ly / l2 + + term2Fourier = self.filterFourier(f, dataFourier=WFDataFourier, test=test) + term2 = self.inverseFourier(dataFourier=term2Fourier) + term2 *= FData + term2Fourier = self.fourier(data=term2) + term2Fourier = self.filterFourier(f, dataFourier=term2Fourier, test=test) + term2Fourier *= 4.0 + # + # add term1 and term2 to get kappa + resultFourier = term1Fourier + term2Fourier + # + if test: + # print "show term 1" + self.plotFourier(term1Fourier) + # print "show term 2" + self.plotFourier(term2Fourier) + # print "show term 1 + term 2" + self.plotFourier(term1Fourier + term2Fourier) + + # get phi from kappa + resultFourier *= -2.0 / self.ell**2 + resultFourier[np.where(np.isfinite(resultFourier) == False)] = 0.0 + if test: + # print "checking the power spectrum of phi map" + self.powerSpectrum(dataFourier=resultFourier, plot=True) + + if test: + # print "Show real-space phi map" + result = self.inverseFourier(resultFourier) + self.plot(result) + + return resultFourier + + def computeQuadEstPhiShearNormalizationFFT( + self, fC0, fCtot, lMin=5.0e2, lMax=3.0e3, test=False + ): + """Multiplicative normalization for phi estimator from shear only, + computed with FFT + ell cuts are performed to remain in the regime L_phi < l_T + """ + + # weight function for shear + def fdLnC0dLnl(ell): + e = 0.01 + lup = ell * (1.0 + e) + ldown = ell * (1.0 - e) + result = fC0(lup) / fC0(ldown) + result = np.log(result) / (2.0 * e) + return result + + if test: + F = np.array(map(fdLnC0dLnl, self.ell.flatten())) + plt.semilogx(self.ell.flatten(), F, "b.") + plt.show() + + # sort of shear Wiener-filter + def f(ell): + # cut off the high ells from input map + if ell < lMin or ell > lMax: + return 0.0 + result = fC0(ell) / fCtot(ell) + result *= fdLnC0dLnl(ell) # for shear + result **= 2.0 + result /= 0.5 + if not np.isfinite(result): + result = 0.0 + return result + + WFDataFourier = self.filterFourierIsotropic( + f, dataFourier=np.ones_like(self.ell), test=test + ) + if test: + # print "showing the WF map" + self.plotFourier(dataFourier=WFDataFourier) + + # mean value of cos^2 is 1/2 + integral = np.sum(0.5 * WFDataFourier) + integral /= self.sizeX * self.sizeY + if test: + integralFFT = self.inverseFourier(WFDataFourier) + integralFFT = ( + integralFFT[0, 0] / 2.0 + ) # I think the factor 2 is about the half/full Fourier plane + # print "Integral from sum=", integral + # print "Same integral from FFT=", integralFFT + + # fill a Fourier map with this value + resultFourier = np.ones_like(WFDataFourier) + resultFourier /= integral + + return resultFourier + + def computeQuadEstPhiShearNormalizationCorrectedFFT( + self, fC0, fCtot, lMin=5.0e2, lMax=3.0e3, test=False, cache=None + ): + """Multiplicative normalization for phi estimator from shear only, + computed with FFT. + This normalization corrects for the multiplicative bias in the estimator. + ell cuts are performed to remain in the regime L_phi < l_T + """ + + def doCalculation(): + # weight function for shear + def fdLnC0dLnl(ell): + e = 0.01 + lup = ell * (1.0 + e) + ldown = ell * (1.0 - e) + result = fC0(lup) / fC0(ldown) + result = np.log(result) / (2.0 * e) + return result + + def f(ell): + # cut off the high ells from input map + if (ell < lMin) or (ell > lMax): + return 0.0 + result = fC0(ell) / fCtot(ell) ** 2 + result *= fdLnC0dLnl(ell) # for shear + result /= 0.5 + if not np.isfinite(result): + result = 0.0 + return result + + # generate shear maps + shearFourier = self.filterFourierIsotropic( + f, dataFourier=np.ones_like(self.ell), test=test + ) + # + cosXFourier = shearFourier * (self.lx**2 - self.ly**2) / self.ell**2 + cosXFourier[np.where(np.isfinite(cosXFourier) == False)] = 0.0 + cosX = self.inverseFourier(cosXFourier) + # + cosYFourier = shearFourier * self.lx * self.ly / self.ell**2 + cosYFourier[np.where(np.isfinite(cosYFourier) == False)] = 0.0 + cosY = self.inverseFourier(cosYFourier) + + # generate gradient C0 map + f = lambda ell: fC0(ell) * (ell >= lMin) * (ell <= lMax) + c0Fourier = self.filterFourierIsotropic( + f, dataFourier=np.ones_like(self.ell), test=test + ) + # the factor i in the gradient makes the Fourier function Hermitian + gradXFourier, gradYFourier = self.computeGradient(dataFourier=c0Fourier) + gradX = self.inverseFourier( + gradXFourier + ) # extra factor of i will be cancelled later + gradY = self.inverseFourier( + gradYFourier + ) # extra factor of i will be cancelled later + + # generate ell limit map + f = lambda ell: (ell >= lMin) * (ell <= lMax) + ellLimitsFourier = self.filterFourierIsotropic( + f, dataFourier=np.ones_like(self.ell), test=test + ) + ellLimits = self.inverseFourier(ellLimitsFourier) + + # Term 1 + grad1XcosXFourier = self.fourier(gradX * cosX) + grad1XcosXFourier *= self.lx**2 - self.ly**2 # for cos + grad1XcosXFourier *= 2.0 * self.lx / 1.0j / self.ell**4 # for grad + grad1XcosXFourier[np.where(np.isfinite(grad1XcosXFourier) == False)] = 0.0 + # + grad1YcosXFourier = self.fourier(gradY * cosX) + grad1YcosXFourier *= self.lx**2 - self.ly**2 # for cos + grad1YcosXFourier *= 2.0 * self.ly / 1.0j / self.ell**4 # for grad + grad1YcosXFourier[np.where(np.isfinite(grad1YcosXFourier) == False)] = 0.0 + # + grad1XcosYFourier = self.fourier(gradX * cosY) + grad1XcosYFourier *= 4.0 * self.lx * self.ly # for cos + grad1XcosYFourier *= 2.0 * self.lx / 1.0j / self.ell**4 # for grad + grad1XcosYFourier[np.where(np.isfinite(grad1XcosYFourier) == False)] = 0.0 + # + grad1YcosYFourier = self.fourier(gradY * cosY) + grad1YcosYFourier *= 4.0 * self.lx * self.ly # for cos + grad1YcosYFourier *= 2.0 * self.ly / 1.0j / self.ell**4 # for grad + grad1YcosYFourier[np.where(np.isfinite(grad1YcosYFourier) == False)] = 0.0 + # + # sum + term1Fourier = ( + grad1XcosXFourier + + grad1YcosXFourier + + grad1XcosYFourier + + grad1YcosYFourier + ) + if test: + # print "showing various terms" + self.plotFourier(grad1XcosXFourier) + self.plotFourier(grad1YcosXFourier) + self.plotFourier(grad1XcosYFourier) + self.plotFourier(grad1YcosYFourier) + self.plotFourier(term1Fourier) + + # Term 2 + grad2XcosX = self.inverseFourier(gradXFourier * cosXFourier) + grad2XcosXFourier = self.fourier(grad2XcosX * ellLimits) + grad2XcosXFourier *= self.lx**2 - self.ly**2 # for cos + grad2XcosXFourier *= 2.0 * self.lx / 1.0j / self.ell**4 # for grad + grad2XcosXFourier[np.where(np.isfinite(grad2XcosXFourier) == False)] = 0.0 + # + grad2YcosX = self.inverseFourier(gradYFourier * cosXFourier) + grad2YcosXFourier = self.fourier(grad2YcosX * ellLimits) + grad2YcosXFourier *= self.lx**2 - self.ly**2 # for cos + grad2YcosXFourier *= 2.0 * self.ly / 1.0j / self.ell**4 # for grad + grad2YcosXFourier[np.where(np.isfinite(grad2YcosXFourier) == False)] = 0.0 + # + grad2XcosY = self.inverseFourier(gradXFourier * cosYFourier) + grad2XcosYFourier = self.fourier(grad2XcosY * ellLimits) + grad2XcosYFourier *= 4.0 * self.lx * self.ly # for cos + grad2XcosYFourier *= 2.0 * self.lx / 1.0j / self.ell**4 # for grad + grad2XcosYFourier[np.where(np.isfinite(grad2XcosYFourier) == False)] = 0.0 + # + grad2YcosY = self.inverseFourier(gradYFourier * cosYFourier) + grad2YcosYFourier = self.fourier(grad2YcosY * ellLimits) + grad2YcosYFourier *= 4.0 * self.lx * self.ly # for cos + grad2YcosYFourier *= 2.0 * self.ly / 1.0j / self.ell**4 # for grad + grad2YcosYFourier[np.where(np.isfinite(grad2YcosXFourier) == False)] = 0.0 + # + # sum + term2Fourier = ( + grad2XcosXFourier + + grad2YcosXFourier + + grad2XcosYFourier + + grad2YcosYFourier + ) + if test: + # print "showing various terms" + self.plotFourier(grad2XcosXFourier) + self.plotFourier(grad2YcosXFourier) + self.plotFourier(grad2XcosYFourier) + self.plotFourier(grad2YcosYFourier) + self.plotFourier(term2Fourier) + + # sum and invert + resultFourier = 1.0 / (term1Fourier + term2Fourier) + resultFourier[np.where(np.isfinite(resultFourier) == False)] = 0.0 + + # remove L > 2 lMax + f = lambda ell: (ell <= 2.0 * lMax) + resultFourier = self.filterFourierIsotropic( + f, dataFourier=resultFourier, test=test + ) + if test: + # print "showing result" + self.plotFourier(resultFourier) + return resultFourier + + # Caching boiler plate + # if no caching is desired, just compute + if cache is None: + resultFourier = doCalculation() + # if caching is desired + else: + # if first call with caching, set up the cache dictionary + if not hasattr( + self.computeQuadEstPhiShearNormalizationCorrectedFFT.__func__, "cache" + ): + self.computeQuadEstPhiShearNormalizationCorrectedFFT.__func__.cache = {} + # if the calculation has been done before + if self.computeQuadEstPhiShearNormalizationCorrectedFFT.cache.has_key( + cache + ): + resultFourier = ( + self.computeQuadEstPhiShearNormalizationCorrectedFFT.cache[ + cache + ].copy() + ) + # if this calculation was not done before + else: + resultFourier = doCalculation() + self.computeQuadEstPhiShearNormalizationCorrectedFFT.cache[ + cache + ] = resultFourier.copy() + + return resultFourier + + def computeQuadEstKappaShearNormCorr( + self, + fC0, + fCtot, + lMin=1.0, + lMax=1.0e5, + dataFourier=None, + dataFourier2=None, + path=None, + corr=True, + test=False, + cache=None, + ): + # non-normalized QE for phi + resultFourier = self.quadEstPhiShearNonNorm( + fC0, + fCtot, + lMin=lMin, + lMax=lMax, + dataFourier=dataFourier, + dataFourier2=dataFourier2, + test=test, + ) + # convert from phi to kappa + resultFourier = self.kappaFromPhi(resultFourier) + # compute normalization + if corr: + normalizationFourier = self.computeQuadEstPhiShearNormalizationCorrectedFFT( + fC0, fCtot, lMin=lMin, lMax=lMax, test=test, cache=cache + ) + else: + normalizationFourier = self.computeQuadEstPhiShearNormalizationFFT( + fC0, fCtot, lMin=lMin, lMax=lMax, test=test + ) + # normalized (not mean field-subtracted) QE for kappa + resultFourier *= normalizationFourier + # save to file if needed + if path is not None: + self.saveDataFourier(resultFourier, path) + return resultFourier + + def computeQuadEstKappaShearNoiseFFT( + self, + fC0, + fCtot, + fCfg=None, + lMin=1.0, + lMax=1.0e5, + corr=True, + test=False, + cache=None, + ): + """Computes the lensing noise power spectrum N_L^kappa + for the shear estimator. + """ + if fCfg is None: + fCfg = fCtot + + def fdLnC0dLnl(ell): + e = 0.01 + lup = ell * (1.0 + e) + ldown = ell * (1.0 - e) + result = fC0(lup) / fC0(ldown) + result = np.log(result) / (2.0 * e) + return result + + def g(ell): + # cut off the high ells from input map + if (ell < lMin) or (ell > lMax): + return 0.0 + result = fC0(ell) / fCtot(ell) ** 2 + result *= fdLnC0dLnl(ell) # for isotropic dilation + result /= 0.5 + if not np.isfinite(result): + result = 0.0 + return result + + # useful terms + fx = (self.lx**2 - self.ly**2) / self.ell**2 + fx[np.where(np.isfinite(fx) == False)] = 0.0 + # + fy = self.lx * self.ly / self.ell**2 + fy[np.where(np.isfinite(fy) == False)] = 0.0 + + # First, the symmetric term + f = lambda ell: g(ell) * fCfg(ell) + term1Fourier = self.filterFourierIsotropic( + f, dataFourier=np.ones_like(self.ell), test=test + ) + # + term1xFourier = term1Fourier * fx + term1x = self.inverseFourier(term1xFourier) + # + term1yFourier = term1Fourier * fy + term1y = self.inverseFourier(term1yFourier) + + term1xxFourier = self.fourier(term1x**2) + term1xxFourier *= fx**2 + # + term1yyFourier = self.fourier(term1y**2) + term1yyFourier *= (4.0 * fy) ** 2 + # + term1xyFourier = self.fourier(term1x * term1y) + term1xyFourier *= fx * 4.0 * fy + term1xyFourier *= 2.0 + # + term1Fourier = term1xxFourier + term1yyFourier + term1xyFourier + if test: + self.plotFourier(term1xxFourier) + self.plotFourier(term1yyFourier) + self.plotFourier(term1xyFourier) + self.plotFourier(term1Fourier) + + # Second, the asymmetric term + f = lambda ell: fCfg(ell) * (ell >= lMin) * (ell <= lMax) + term2bFourier = self.filterFourierIsotropic( + f, dataFourier=np.ones_like(self.ell), test=test + ) + term2b = self.inverseFourier(term2bFourier) + # + f = lambda ell: g(ell) ** 2 * fCfg(ell) + term2aFourier = self.filterFourierIsotropic( + f, dataFourier=np.ones_like(self.ell), test=test + ) + + term2xxFourier = term2aFourier * fx**2 + term2xx = self.inverseFourier(term2xxFourier) + term2xxFourier = self.fourier(term2xx * term2b) + term2xxFourier *= fx**2 + # + term2yyFourier = term2aFourier * fy**2 + term2yy = self.inverseFourier(term2yyFourier) + term2yyFourier = self.fourier(term2yy * term2b) + term2yyFourier *= (4.0 * fy) ** 2 + # + term2xyFourier = term2aFourier * fx * fy + term2xy = self.inverseFourier(term2xyFourier) + term2xyFourier = self.fourier(term2xy * term2b) + term2xyFourier *= fx * 4.0 * fy + term2xyFourier *= 2.0 + # + term2Fourier = term2xxFourier + term2yyFourier + term2xyFourier + if test: + self.plotFourier(term2xxFourier) + self.plotFourier(term2yyFourier) + self.plotFourier(term2xyFourier) + self.plotFourier(term2Fourier) + + # add terms + resultFourier = term1Fourier + term2Fourier + + # compute normalization + if corr: + normalizationFourier = self.computeQuadEstPhiShearNormalizationCorrectedFFT( + fC0, fCtot, lMin=lMin, lMax=lMax, test=test, cache=cache + ) + else: + normalizationFourier = self.computeQuadEstPhiShearNormalizationFFT( + fC0, fCtot, lMin=lMin, lMax=lMax, test=test + ) + + # divide by squared normalization + resultFourier *= normalizationFourier**2 + + # remove L > 2 lMax + f = lambda ell: (ell <= 2.0 * lMax) + resultFourier = self.filterFourierIsotropic( + f, dataFourier=resultFourier, test=test + ) + + return resultFourier + + def forecastN0KappaShear( + self, + fC0, + fCtot, + fCfg=None, + lMin=1.0, + lMax=1.0e5, + corr=True, + test=False, + cache=None, + ): + """Interpolates the result for N_L^{kappa_shear} = f(l), + to be used for forecasts on lensing reconstruction. + """ + # print "computing the reconstruction noise" + n0Kappa = self.computeQuadEstKappaShearNoiseFFT( + fC0, + fCtot, + fCfg=fCfg, + lMin=lMin, + lMax=lMax, + corr=corr, + test=test, + cache=cache, + ) + # keep only the real part (the imag. part should be zero, but is tiny in practice) + n0Kappa = np.real(n0Kappa) + # remove the nans + n0Kappa = np.nan_to_num(n0Kappa) + # make sure every value is positive + n0Kappa = np.abs(n0Kappa) + + # fix the issue of the wrong ell=0 value + # replace it by the value lx=0, ly=fundamental + n0Kappa[0, 0] = n0Kappa[0, 1] + + # interpolate + where = (self.ell.flatten() > 0.0) * (self.ell.flatten() < 2.0 * lMax) + L = self.ell.flatten()[where] + N = n0Kappa.flatten()[where] + lnfln = interp1d( + np.log(L), np.log(N), kind="linear", bounds_error=False, fill_value=np.inf + ) + f = lambda ell: np.exp(lnfln(np.log(ell))) + return f + + ############################################################################### + # Lensing multiplicative bias from lensed foregrounds + + def computeMultBiasLensedForegroundsShear( + self, fC0, fCtot, fCfgBias, lMin=5.0e2, lMax=3.0e3, test=False, cache=None + ): + """Multiplicative bias to the shear estimator + from lensed foregrounds in the CMB map. + fCfgBias: unlensed foreground power spectrum. + Such that: + = kappa_CMB + (multiplicative bias) * kappa_foreground + noise. + """ + + # weight function for shear + def fdLnC0dLnl(ell): + e = 0.01 + lup = ell * (1.0 + e) + ldown = ell * (1.0 - e) + result = fC0(lup) / fC0(ldown) + result = np.log(result) / (2.0 * e) + return result + + def f(ell): + # cut off the high ells from input map + if (ell < lMin) or (ell > lMax): + return 0.0 + result = fC0(ell) / fCtot(ell) ** 2 + result *= fdLnC0dLnl(ell) # for shear + result /= 0.5 + if not np.isfinite(result): + result = 0.0 + return result + + # generate shear maps + shearFourier = self.filterFourierIsotropic( + f, dataFourier=np.ones_like(self.ell), test=test + ) + # + cosXFourier = shearFourier * (self.lx**2 - self.ly**2) / self.ell**2 + cosXFourier[np.where(np.isfinite(cosXFourier) == False)] = 0.0 + cosX = self.inverseFourier(cosXFourier) + # + cosYFourier = shearFourier * self.lx * self.ly / self.ell**2 + cosYFourier[np.where(np.isfinite(cosYFourier) == False)] = 0.0 + cosY = self.inverseFourier(cosYFourier) + + # generate gradient C0 map + f = lambda ell: fCfgBias(ell) * (ell >= lMin) * (ell <= lMax) + c0Fourier = self.filterFourierIsotropic( + f, dataFourier=np.ones_like(self.ell), test=test + ) + # the factor i in the gradient makes the Fourier function Hermitian + gradXFourier, gradYFourier = self.computeGradient(dataFourier=c0Fourier) + gradX = self.inverseFourier( + gradXFourier + ) # extra factor of i will be cancelled later + gradY = self.inverseFourier( + gradYFourier + ) # extra factor of i will be cancelled later + + # generate ell limit map + f = lambda ell: (ell >= lMin) * (ell <= lMax) + ellLimitsFourier = self.filterFourierIsotropic( + f, dataFourier=np.ones_like(self.ell), test=test + ) + ellLimits = self.inverseFourier(ellLimitsFourier) + + # Term 1 + grad1XcosXFourier = self.fourier(gradX * cosX) + grad1XcosXFourier *= self.lx**2 - self.ly**2 # for cos + grad1XcosXFourier *= 2.0 * self.lx / 1.0j / self.ell**4 # for grad + grad1XcosXFourier[np.where(np.isfinite(grad1XcosXFourier) == False)] = 0.0 + # + grad1YcosXFourier = self.fourier(gradY * cosX) + grad1YcosXFourier *= self.lx**2 - self.ly**2 # for cos + grad1YcosXFourier *= 2.0 * self.ly / 1.0j / self.ell**4 # for grad + grad1YcosXFourier[np.where(np.isfinite(grad1YcosXFourier) == False)] = 0.0 + # + grad1XcosYFourier = self.fourier(gradX * cosY) + grad1XcosYFourier *= 4.0 * self.lx * self.ly # for cos + grad1XcosYFourier *= 2.0 * self.lx / 1.0j / self.ell**4 # for grad + grad1XcosYFourier[np.where(np.isfinite(grad1XcosYFourier) == False)] = 0.0 + # + grad1YcosYFourier = self.fourier(gradY * cosY) + grad1YcosYFourier *= 4.0 * self.lx * self.ly # for cos + grad1YcosYFourier *= 2.0 * self.ly / 1.0j / self.ell**4 # for grad + grad1YcosYFourier[np.where(np.isfinite(grad1YcosYFourier) == False)] = 0.0 + # + # sum + term1Fourier = ( + grad1XcosXFourier + + grad1YcosXFourier + + grad1XcosYFourier + + grad1YcosYFourier + ) + if test: + # print "showing various terms" + self.plotFourier(grad1XcosXFourier) + self.plotFourier(grad1YcosXFourier) + self.plotFourier(grad1XcosYFourier) + self.plotFourier(grad1YcosYFourier) + self.plotFourier(term1Fourier) + + # Term 2 + grad2XcosX = self.inverseFourier(gradXFourier * cosXFourier) + grad2XcosXFourier = self.fourier(grad2XcosX * ellLimits) + grad2XcosXFourier *= self.lx**2 - self.ly**2 # for cos + grad2XcosXFourier *= 2.0 * self.lx / 1.0j / self.ell**4 # for grad + grad2XcosXFourier[np.where(np.isfinite(grad2XcosXFourier) == False)] = 0.0 + # + grad2YcosX = self.inverseFourier(gradYFourier * cosXFourier) + grad2YcosXFourier = self.fourier(grad2YcosX * ellLimits) + grad2YcosXFourier *= self.lx**2 - self.ly**2 # for cos + grad2YcosXFourier *= 2.0 * self.ly / 1.0j / self.ell**4 # for grad + grad2YcosXFourier[np.where(np.isfinite(grad2YcosXFourier) == False)] = 0.0 + # + grad2XcosY = self.inverseFourier(gradXFourier * cosYFourier) + grad2XcosYFourier = self.fourier(grad2XcosY * ellLimits) + grad2XcosYFourier *= 4.0 * self.lx * self.ly # for cos + grad2XcosYFourier *= 2.0 * self.lx / 1.0j / self.ell**4 # for grad + grad2XcosYFourier[np.where(np.isfinite(grad2XcosYFourier) == False)] = 0.0 + # + grad2YcosY = self.inverseFourier(gradYFourier * cosYFourier) + grad2YcosYFourier = self.fourier(grad2YcosY * ellLimits) + grad2YcosYFourier *= 4.0 * self.lx * self.ly # for cos + grad2YcosYFourier *= 2.0 * self.ly / 1.0j / self.ell**4 # for grad + grad2YcosYFourier[np.where(np.isfinite(grad2YcosXFourier) == False)] = 0.0 + # + # sum + term2Fourier = ( + grad2XcosXFourier + + grad2YcosXFourier + + grad2XcosYFourier + + grad2YcosYFourier + ) + if test: + # print "showing various terms" + self.plotFourier(grad2XcosXFourier) + self.plotFourier(grad2YcosXFourier) + self.plotFourier(grad2XcosYFourier) + self.plotFourier(grad2YcosYFourier) + self.plotFourier(term2Fourier) + + # sum + resultFourier = term1Fourier + term2Fourier + resultFourier[np.where(np.isfinite(resultFourier) == False)] = 0.0 + + # remove L > 2 lMax + f = lambda ell: (ell <= 2.0 * lMax) + resultFourier = self.filterFourierIsotropic( + f, dataFourier=resultFourier, test=test + ) + + # normalize with the standard normalization + resultFourier *= self.computeQuadEstPhiShearNormalizationCorrectedFFT( + fC0, fCtot, lMin=lMin, lMax=lMax, test=False, cache=cache + ) + if test: + # print "showing result" + self.plotFourier(resultFourier) + + return resultFourier + + def forecastMultBiasLensedForegroundsShear( + self, fC0, fCtot, fCfgBias, lMin=1.0, lMax=1.0e5, test=False + ): + """Interpolates the multiplicative bias to the shear estimator + due to lensed foregrounds in the map + """ + # print "computing the multiplicative bias" + result = self.computeMultBiasLensedForegroundsShear( + fC0, fCtot, fCfgBias, lMin=lMin, lMax=lMax, test=test + ) + # keep only the real part (the imag. part should be zero, but is tiny in practice) + result = np.real(result) + # remove the nans + result = np.nan_to_num(result) + + # fix the issue of the wrong ell=0 value + # replace it by the value lx=0, ly=fundamental + result[0, 0] = result[0, 1] + + # interpolate, preserving the sign + lnfln = interp1d( + np.log(self.ell.flatten()), + np.log(np.abs(result).flatten()), + kind="linear", + bounds_error=False, + fill_value=0.0, + ) + signln = interp1d( + np.log(self.ell.flatten()), + np.sign(result.flatten()), + kind="linear", + bounds_error=False, + fill_value=0.0, + ) + f = lambda ell: np.exp(lnfln(np.log(ell))) * signln(np.log(ell)) + return f + + ############################################################################### + ############################################################################### + # Noise cross-power between shear and dilation + + def computeQuadEstKappaShearDilationNoiseFFT( + self, + fC0, + fCtot, + fCfg=None, + lMin=1.0, + lMaxS=1.0e5, + lMaxD=1.0e5, + corr=True, + test=False, + cache=None, + ): + """Computes the lensing noise cross-spectrum N_L^{kappa_shear * kappa_dilation} + between the shear and the dilation estimators. + lMaxS: maximum temperature multipole used in the shear estimator + lMaxD: maximum temperature multipole used in the dilation estimator + """ + if fCfg is None: + fCfg = fCtot + + # Numerator: keep the minimum lMax + lMax = min(lMaxS, lMaxD) + + def fdLnC0dLnl(ell): + e = 0.01 + lup = ell * (1.0 + e) + ldown = ell * (1.0 - e) + result = fC0(lup) / fC0(ldown) + result = np.log(result) / (2.0 * e) + return result + + def fdLnl2C0dLnl(ell): + e = 0.01 + lup = ell * (1.0 + e) + ldown = ell * (1.0 - e) + result = lup**2 * fC0(lup) + result /= ldown**2 * fC0(ldown) + result = np.log(result) / (2.0 * e) + return result + + def gShear(ell): + # cut off the high ells from input map + if (ell < lMin) or (ell > lMax): + return 0.0 + result = fC0(ell) / fCtot(ell) ** 2 + result *= fdLnC0dLnl(ell) # for isotropic dilation + result /= 0.5 + if not np.isfinite(result): + result = 0.0 + return result + + def gDilation(ell): + # cut off the high ells from input map + if (ell < lMin) or (ell > lMax): + return 0.0 + result = fC0(ell) / fCtot(ell) ** 2 + result *= fdLnl2C0dLnl(ell) # for isotropic dilation + result /= 0.5 + if not np.isfinite(result): + result = 0.0 + return result + + # useful terms + fx = (self.lx**2 - self.ly**2) / self.ell**2 + fx[np.where(np.isfinite(fx) == False)] = 0.0 + # + fy = self.lx * self.ly / self.ell**2 + fy[np.where(np.isfinite(fy) == False)] = 0.0 + + # First, the symmetric term + f = lambda ell: gDilation(ell) * fCfg(ell) + term1DFourier = self.filterFourierIsotropic( + f, dataFourier=np.ones_like(self.ell), test=test + ) + term1D = self.inverseFourier(term1DFourier) + # + f = lambda ell: gShear(ell) * fCfg(ell) + term1SFourier = self.filterFourierIsotropic( + f, dataFourier=np.ones_like(self.ell), test=test + ) + + term1xFourier = term1SFourier * fx + term1x = self.inverseFourier(term1xFourier) + term1xFourier = self.fourier(term1x * term1D) + term1xFourier *= fx + # + term1yFourier = term1SFourier * fy + term1y = self.inverseFourier(term1yFourier) + term1yFourier = self.fourier(term1y * term1D) + term1yFourier *= 4.0 * fy + # + term1Fourier = term1xFourier + term1yFourier + if test: + self.plotFourier(term1xFourier) + self.plotFourier(term1yFourier) + self.plotFourier(term1Fourier) + + # Second, the asymmetric term + f = lambda ell: fCfg(ell) * (ell >= lMin) * (ell <= lMax) + term2CFourier = self.filterFourierIsotropic( + f, dataFourier=np.ones_like(self.ell), test=test + ) + term2C = self.inverseFourier(term2CFourier) + # + f = lambda ell: gShear(ell) * gDilation(ell) * fCfg(ell) + term2SDFourier = self.filterFourierIsotropic( + f, dataFourier=np.ones_like(self.ell), test=test + ) + + term2xFourier = term2SDFourier * fx + term2x = self.inverseFourier(term2xFourier) + term2xFourier = self.fourier(term2x * term2C) + term2xFourier *= fx + # + term2yFourier = term2SDFourier * fy + term2y = self.inverseFourier(term2yFourier) + term2yFourier = self.fourier(term2y * term2C) + term2yFourier *= 4.0 * fy + # + term2Fourier = term2xFourier + term2yFourier + if test: + self.plotFourier(term2xFourier) + self.plotFourier(term2yFourier) + self.plotFourier(term2Fourier) + self.plotFourier(term1Fourier + term2Fourier) + + # add terms + resultFourier = term1Fourier + term2Fourier + + # compute normalization + if corr: + normalizationFourier = self.computeQuadEstPhiShearNormalizationCorrectedFFT( + fC0, fCtot, lMin=lMin, lMax=lMaxS, test=test, cache=cache + ) + normalizationFourier *= ( + self.computeQuadEstPhiDilationNormalizationCorrectedFFT( + fC0, fCtot, lMin=lMin, lMax=lMaxD, test=test, cache=cache + ) + ) + else: + normalizationFourier = self.computeQuadEstPhiShearNormalizationFFT( + fC0, fCtot, lMin=lMin, lMax=lMaxS, test=test + ) + normalizationFourier = self.computeQuadEstPhiDilationNormalizationFFT( + fC0, fCtot, lMin=lMin, lMax=lMaxD, test=test + ) + + # divide by squared normalization + resultFourier *= normalizationFourier + + # remove L > 2 lMax + f = lambda ell: (ell <= 2.0 * lMax) + resultFourier = self.filterFourierIsotropic( + f, dataFourier=resultFourier, test=test + ) + + return resultFourier + + def forecastN0KappaShearDilation( + self, + fC0, + fCtot, + fCfg=None, + lMin=1.0, + lMaxS=1.0e5, + lMaxD=1.0e5, + corr=True, + test=False, + cache=None, + ): + """Interpolates the result for N_L^{kappa_shear * kappa_dilation} = f(l), + to be used for forecasts on lensing reconstruction. + lMaxS: maximum temperature multipole used in the shear estimator + lMaxD: maximum temperature multipole used in the dilation estimator + """ + # One can only form the cross-correlation for lensing modes < 2*lMax, + # where lMax = min(lMaxS, lMaxD). + lMax = min(lMaxS, lMaxD) + + # print "computing the reconstruction noise" + n0Kappa = self.computeQuadEstKappaShearDilationNoiseFFT( + fC0, + fCtot, + fCfg=fCfg, + lMin=lMin, + lMaxS=lMaxS, + lMaxD=lMaxD, + corr=corr, + test=test, + cache=cache, + ) + # keep only the real part (the imag. part should be zero, but is tiny in practice) + n0Kappa = np.real(n0Kappa) + # remove the nans + n0Kappa = np.nan_to_num(n0Kappa) + + # fix the issue of the wrong ell=0 value + # replace it by the value lx=0, ly=fundamental + n0Kappa[0, 0] = n0Kappa[0, 1] + + # interpolate + where = (self.ell.flatten() > 0.0) * (self.ell.flatten() < 2.0 * lMax) + L = self.ell.flatten()[where] + N = n0Kappa.flatten()[where] + lnfln = interp1d( + np.log(L), + np.log(np.abs(N)), + kind="linear", + bounds_error=False, + fill_value=np.inf, + ) + sgnfln = interp1d( + np.log(L), np.sign(N), kind="linear", bounds_error=False, fill_value=np.inf + ) + f = lambda ell: np.exp(lnfln(np.log(ell))) * sgnfln(np.log(ell)) + return f + + ############################################################################### + ############################################################################### + # Shear B-mode estimator + + def quadEstPhiShearBNonNorm( + self, + fC0, + fCtot, + lMin=5.0e2, + lMax=3000.0, + dataFourier=None, + dataFourier2=None, + test=False, + ): + """Non-normalized quadratic estimator for phi from shear B-mode only + fC0: ulensed power spectrum + fCtot: lensed power spectrum + noise + """ + if dataFourier is None: + dataFourier = self.dataFourier.copy() + if dataFourier2 is None: + dataFourier2 = dataFourier.copy() + + # cut off high ells + f = lambda ell: 1.0 * (ell >= lMin) * (ell <= lMax) + FDataFourier = self.filterFourierIsotropic( + f, dataFourier=dataFourier, test=test + ) + FData = self.inverseFourier(FDataFourier) + if test: + # print "show Fourier data" + self.plotFourier(dataFourier=FDataFourier) + + # weight function for shear + def fdLnC0dLnl(ell): + e = 0.01 + lup = ell * (1.0 + e) + ldown = ell * (1.0 - e) + result = fC0(lup) / fC0(ldown) + result = np.log(result) / (2.0 * e) + return result + + if test: + F = np.array(map(fdLnC0dLnl, self.ell.flatten())) + plt.semilogx(self.ell.flatten(), F, "b.") + plt.show() + + # sort of shear Wiener-filter + def f(ell): + # cut off the high ells from input map + if (ell < lMin) or (ell > lMax): + return 0.0 + result = fC0(ell) / fCtot(ell) ** 2 + result *= fdLnC0dLnl(ell) # for shear + result /= 0.5 + if not np.isfinite(result): + result = 0.0 + return result + + WFDataFourier = self.filterFourierIsotropic( + f, dataFourier=dataFourier2, test=test + ) + if test: + # print "showing the WF map" + self.plotFourier(dataFourier=WFDataFourier) + # print "checking the power spectrum of this map" + theory = lambda ell: f(ell) ** 2 * fC0(ell) + self.powerSpectrum(theory=[theory], dataFourier=WFDataFourier, plot=True) + + # multiplication by sin 2 theta_{L,l} + + def fDiff(lx, ly): + l2 = lx**2 + ly**2 + if l2 == 0: + return 0.0 + else: + return (lx**2 - ly**2) / l2 + + def fProd(lx, ly): + l2 = lx**2 + ly**2 + if l2 == 0: + return 0.0 + else: + return lx * ly / l2 + + # term 1 + term1Fourier = self.filterFourier(fProd, dataFourier=WFDataFourier, test=test) + term1 = self.inverseFourier(dataFourier=term1Fourier) + term1 *= FData + term1Fourier = self.fourier(data=term1) + term1Fourier = self.filterFourier(fDiff, dataFourier=term1Fourier, test=test) + term1Fourier *= 2.0 + + # term 2 + term2Fourier = self.filterFourier(fDiff, dataFourier=WFDataFourier, test=test) + term2Fourier *= -1.0 # because we want ly^2-lx^2 here + term2 = self.inverseFourier(dataFourier=term2Fourier) + term2 *= FData + term2Fourier = self.fourier(data=term2) + term2Fourier = self.filterFourier(fProd, dataFourier=term2Fourier, test=test) + term2Fourier *= 2.0 + + # add term1 and term2 to get kappa + resultFourier = term1Fourier + term2Fourier + # + if test: + # print "show term 1" + self.plotFourier(term1Fourier) + # print "show term 2" + self.plotFourier(term2Fourier) + # print "show term 1 + term 2" + self.plotFourier(term1Fourier + term2Fourier) + + # get phi from kappa + resultFourier *= -2.0 / self.ell**2 + resultFourier[np.where(np.isfinite(resultFourier) == False)] = 0.0 + if test: + # print "checking the power spectrum of phi map" + self.powerSpectrum(dataFourier=resultFourier, plot=True) + + if test: + # print "Show real-space phi map" + result = self.inverseFourier(resultFourier) + self.plot(result) + + return resultFourier + + def computeQuadEstPhiShearBNormalizationFFT( + self, fC0, fCtot, lMin=5.0e2, lMax=3.0e3, test=False + ): + """Multiplicative normalization for phi estimator from shear B-mode only, + computed with FFT + ell cuts are performed to remain in the regime L_phi < l_T. + Same as for shear E-mode, since = = 1/2. + """ + resultFourier = self.computeQuadEstPhiShearNormalizationFFT( + fC0, fCtot, lMin=lMin, lMax=lMax, test=test + ) + return resultFourier + + def computeQuadEstPhiShearBNormalizationCorrectedFFT( + self, fC0, fCtot, lMin=5.0e2, lMax=3.0e3, test=False, cache=None + ): + """The shear B-mode estimator has zero response to lensing, + so the corrected normalization would be to divide by zero, + i.e. the unbiased normalization should be +infinity. + Instead, here we decide to use the same normalization as the shear E-mode, + so that shear E-mode and B-modes can be compared. + """ + resultFourier = self.computeQuadEstPhiShearNormalizationCorrectedFFT( + fC0, fCtot, lMin=lMin, lMax=lMax, test=test, cache=cache + ) + return resultFourier + + def computeQuadEstKappaShearBNormCorr( + self, + fC0, + fCtot, + lMin=1.0, + lMax=1.0e5, + dataFourier=None, + dataFourier2=None, + path=None, + corr=True, + test=False, + cache=None, + ): + # non-normalized QE for phi + resultFourier = self.quadEstPhiShearBNonNorm( + fC0, + fCtot, + lMin=lMin, + lMax=lMax, + dataFourier=dataFourier, + dataFourier2=dataFourier2, + test=test, + ) + # convert from phi to kappa + resultFourier = self.kappaFromPhi(resultFourier) + # compute normalization + if corr: + normalizationFourier = ( + self.computeQuadEstPhiShearBNormalizationCorrectedFFT( + fC0, fCtot, lMin=lMin, lMax=lMax, test=test, cache=cache + ) + ) + else: + normalizationFourier = self.computeQuadEstPhiShearBNormalizationFFT( + fC0, fCtot, lMin=lMin, lMax=lMax, test=test + ) + # normalized (not mean field-subtracted) QE for kappa + resultFourier *= normalizationFourier + # save to file if needed + if path is not None: + self.saveDataFourier(resultFourier, path) + return resultFourier + + def computeQuadEstKappaShearBNoiseFFT( + self, + fC0, + fCtot, + fCfg=None, + lMin=1.0, + lMax=1.0e5, + corr=True, + test=False, + cache=None, + ): + """Computes the lensing noise power spectrum N_L^kappa + for the shear B-mode estimator. + """ + if fCfg is None: + fCfg = fCtot + + def fdLnC0dLnl(ell): + e = 0.01 + lup = ell * (1.0 + e) + ldown = ell * (1.0 - e) + result = fC0(lup) / fC0(ldown) + result = np.log(result) / (2.0 * e) + return result + + def g(ell): + # cut off the high ells from input map + if (ell < lMin) or (ell > lMax): + return 0.0 + result = fC0(ell) / fCtot(ell) ** 2 + result *= fdLnC0dLnl(ell) # for isotropic dilation + result /= 0.5 + if not np.isfinite(result): + result = 0.0 + return result + + # useful terms + fDiff = (self.lx**2 - self.ly**2) / self.ell**2 + fDiff[np.where(np.isfinite(fDiff) == False)] = 0.0 + # + fProd = self.lx * self.ly / self.ell**2 + fProd[np.where(np.isfinite(fProd) == False)] = 0.0 + + # First, the symmetric term + f = lambda ell: g(ell) * fCfg(ell) + term1Fourier = self.filterFourierIsotropic( + f, dataFourier=np.ones_like(self.ell), test=test + ) + # + term1xFourier = term1Fourier * fDiff + term1xFourier *= -1.0 + term1x = self.inverseFourier(term1xFourier) + # + term1yFourier = term1Fourier * fProd + term1y = self.inverseFourier(term1yFourier) + + term1xxFourier = self.fourier(term1x**2) + term1xxFourier *= fProd**2 + term1xxFourier *= 4.0 + # + term1yyFourier = self.fourier(term1y**2) + term1yyFourier *= fDiff**2 + term1yyFourier *= 4.0 + # + term1xyFourier = self.fourier(term1x * term1y) + term1xyFourier *= fDiff * fProd + term1xyFourier *= 4.0 * 2.0 + # + term1Fourier = term1xxFourier + term1yyFourier + term1xyFourier + if test: + self.plotFourier(term1xxFourier) + self.plotFourier(term1yyFourier) + self.plotFourier(term1xyFourier) + self.plotFourier(term1Fourier) + + # Second, the asymmetric term + f = lambda ell: fCfg(ell) * (ell >= lMin) * (ell <= lMax) + term2bFourier = self.filterFourierIsotropic( + f, dataFourier=np.ones_like(self.ell), test=test + ) + term2b = self.inverseFourier(term2bFourier) + # + f = lambda ell: g(ell) ** 2 * fCfg(ell) + term2aFourier = self.filterFourierIsotropic( + f, dataFourier=np.ones_like(self.ell), test=test + ) + + term2xxFourier = term2aFourier * fProd**2 + term2xx = self.inverseFourier(term2xxFourier) + term2xxFourier = self.fourier(term2xx * term2b) + term2xxFourier *= fDiff**2 + term2xxFourier *= 4.0 + # + term2yyFourier = term2aFourier * fDiff**2 + term2yy = self.inverseFourier(term2yyFourier) + term2yyFourier = self.fourier(term2yy * term2b) + term2yyFourier *= fProd**2 + term2yyFourier *= 4.0 + # + term2xyFourier = term2aFourier * fDiff * fProd + term2xyFourier *= -1.0 + term2xy = self.inverseFourier(term2xyFourier) + term2xyFourier = self.fourier(term2xy * term2b) + term2xyFourier *= fDiff * fProd + term2xyFourier *= 4.0 + term2xyFourier *= 2.0 + # + term2Fourier = term2xxFourier + term2yyFourier + term2xyFourier + if test: + self.plotFourier(term2xxFourier) + self.plotFourier(term2yyFourier) + self.plotFourier(term2xyFourier) + self.plotFourier(term2Fourier) + self.plotFourier(term1Fourier + term2Fourier) + + # add terms + resultFourier = term1Fourier + term2Fourier + + # compute normalization + if corr: + normalizationFourier = ( + self.computeQuadEstPhiShearBNormalizationCorrectedFFT( + fC0, fCtot, lMin=lMin, lMax=lMax, test=test, cache=cache + ) + ) + else: + normalizationFourier = self.computeQuadEstPhiShearBNormalizationFFT( + fC0, fCtot, lMin=lMin, lMax=lMax, test=test + ) + + # divide by squared normalization + resultFourier *= normalizationFourier**2 + + # remove L > 2 lMax + f = lambda ell: (ell <= 2.0 * lMax) + resultFourier = self.filterFourierIsotropic( + f, dataFourier=resultFourier, test=test + ) + + return resultFourier + + def forecastN0KappaShearB( + self, + fC0, + fCtot, + fCfg=None, + lMin=1.0, + lMax=1.0e5, + corr=True, + test=False, + cache=None, + ): + """Interpolates the result for N_L^{kappa_shearB} = f(l), + to be used for forecasts on lensing reconstruction. + """ + # print "computing the reconstruction noise" + n0Kappa = self.computeQuadEstKappaShearBNoiseFFT( + fC0, + fCtot, + fCfg=fCfg, + lMin=lMin, + lMax=lMax, + corr=corr, + test=test, + cache=cache, + ) + # keep only the real part (the imag. part should be zero, but is tiny in practice) + n0Kappa = np.real(n0Kappa) + # remove the nans + n0Kappa = np.nan_to_num(n0Kappa) + # make sure every value is positive + n0Kappa = np.abs(n0Kappa) + + # fix the issue of the wrong ell=0 value + # replace it by the value lx=0, ly=fundamental + n0Kappa[0, 0] = n0Kappa[0, 1] + + # interpolate + where = (self.ell.flatten() > 0.0) * (self.ell.flatten() < 2.0 * lMax) + L = self.ell.flatten()[where] + N = n0Kappa.flatten()[where] + lnfln = interp1d( + np.log(L), np.log(N), kind="linear", bounds_error=False, fill_value=np.inf + ) + f = lambda ell: np.exp(lnfln(np.log(ell))) + return f + + ############################################################################### + ############################################################################### + # Noise cross-spectrum for Shear E and Shear B estimators + #!!! Not working: gives weird results, that are not independent of the direction of L. + + def computeQuadEstKappaShearShearBNoiseFFT( + self, + fC0, + fCtot, + fCfg=None, + lMin=1.0, + lMaxS=1.0e5, + lMaxSB=1.0e5, + corr=True, + test=False, + cache=None, + ): + """Computes the lensing noise power spectrum N_L^kappa""" + if fCfg is None: + fCfg = fCtot + + # Numerator: keep the minimum lMax + lMax = min(lMaxS, lMaxSB) + + def fdLnC0dLnl(ell): + e = 0.01 + lup = ell * (1.0 + e) + ldown = ell * (1.0 - e) + result = fC0(lup) / fC0(ldown) + result = np.log(result) / (2.0 * e) + return result + + def g(ell): + # cut off the high ells from input map + if (ell < lMin) or (ell > lMax): + return 0.0 + result = fC0(ell) / fCtot(ell) ** 2 + result *= fdLnC0dLnl(ell) # for isotropic dilation + result /= 0.5 + if not np.isfinite(result): + result = 0.0 + return result + + # useful terms + fDiff = (self.lx**2 - self.ly**2) / self.ell**2 + fDiff[np.where(np.isfinite(fDiff) == False)] = 0.0 + # + fProd = self.lx * self.ly / self.ell**2 + fProd[np.where(np.isfinite(fProd) == False)] = 0.0 + + # First, the symmetric term + f = lambda ell: g(ell) * fCfg(ell) + termSFourier = self.filterFourierIsotropic( + f, dataFourier=np.ones_like(self.ell), test=test + ) + + termEx = self.inverseFourier(termSFourier * fDiff) + termEy = self.inverseFourier(termSFourier * 4.0 * fProd) + # + termBx = self.inverseFourier(termSFourier * 2.0 * fProd) + termBy = self.inverseFourier(termSFourier * (-2.0 * fDiff)) + + term1xxFourier = self.fourier(termEx * termBx) + term1xxFourier *= fDiff**2 + # + term1yyFourier = self.fourier(termEy * termBy) + term1yyFourier *= fProd**2 + # + term1xyFourier = self.fourier(termEx * termBy) + term1xyFourier *= fDiff * fProd + # + term1yxFourier = self.fourier(termEy * termBx) + term1yxFourier *= fProd * fDiff + # + term1Fourier = term1xxFourier + term1yyFourier + term1xyFourier + term1yxFourier + if test: + self.plotFourier(term1xxFourier) + self.plotFourier(term1yyFourier) + self.plotFourier(term1xyFourier) + self.plotFourier(term1yxFourier) + self.plotFourier(term1Fourier) + + # Second, the asymmetric term + f = lambda ell: g(ell) ** 2 * fCfg(ell) + term2aFourier = self.filterFourierIsotropic( + f, dataFourier=np.ones_like(self.ell), test=test + ) + # + f = lambda ell: fCfg(ell) * (ell >= lMin) * (ell <= lMax) + term2bFourier = self.filterFourierIsotropic( + f, dataFourier=np.ones_like(self.ell), test=test + ) + term2b = self.inverseFourier(term2bFourier) + + term2xx = self.inverseFourier(term2aFourier * fDiff * 2.0 * fProd) + term2xxFourier = self.fourier(term2xx * term2b) + term2xxFourier *= fDiff**2 + # + term2yy = self.inverseFourier(term2aFourier * 4.0 * fProd * (-2.0 * fDiff)) + term2yyFourier = self.fourier(term2yy * term2b) + term2yyFourier *= fProd**2 + # + term2xy = self.inverseFourier(term2aFourier * fDiff * (-2.0 * fDiff)) + term2xyFourier = self.fourier(term2xy * term2b) + term2xyFourier *= fDiff * fProd + # + term2yx = self.inverseFourier(term2aFourier * 4.0 * fProd * 2.0 * fProd) + term2yxFourier = self.fourier(term2yx * term2b) + term2yxFourier *= fProd * fDiff + # + term2Fourier = term2xxFourier + term2yyFourier + term2xyFourier + term2yxFourier + if test: + self.plotFourier(term2xxFourier) + self.plotFourier(term2yyFourier) + self.plotFourier(term2xyFourier) + self.plotFourier(term2yxFourier) + self.plotFourier(term2Fourier) + + # add terms + resultFourier = term1Fourier + term2Fourier + + if test: + self.plotFourier(resultFourier) + plt.loglog(self.ell.flatten(), resultFourier.flatten(), ".") + plt.show() + + # compute normalization: same for shear and shear B + if corr: + normalizationFourier = self.computeQuadEstPhiShearNormalizationCorrectedFFT( + fC0, fCtot, lMin=lMin, lMax=lMaxS, test=test, cache=cache + ) + if lMaxS == lMaxSB: + normalizationFourier *= normalizationFourier + else: + normalizationFourier *= ( + self.computeQuadEstPhiShearNormalizationCorrectedFFT( + fC0, fCtot, lMin=lMin, lMax=lMaxSB, test=test, cache=cache + ) + ) + + else: + normalizationFourier = self.computeQuadEstPhiShearNormalizationFFT( + fC0, fCtot, lMin=lMin, lMax=lMaxS, test=test + ) + if lMaxS == lMaxSB: + normalizationFourier *= normalizationFourier + else: + normalizationFourier *= self.computeQuadEstPhiShearNormalizationFFT( + fC0, fCtot, lMin=lMin, lMax=lMaxSB, test=test + ) + + # divide by squared normalization + resultFourier *= normalizationFourier + + # remove L > 2 lMax + f = lambda ell: (ell <= 2.0 * lMax) + resultFourier = self.filterFourierIsotropic( + f, dataFourier=resultFourier, test=test + ) + + return resultFourier + + def forecastN0KappaShearShearB( + self, + fC0, + fCtot, + fCfg=None, + lMin=1.0, + lMaxS=1.0e5, + lMaxSB=1.0e5, + corr=True, + test=False, + cache=None, + ): + """Interpolates the result for N_L^{kappa_shear * kappa_shearB} = f(l), + to be used for forecasts on lensing reconstruction. + """ + # One can only form the cross-correlation for lensing modes < 2*lMax, + # where lMax = min(lMaxS, lMaxD). + lMax = min(lMaxS, lMaxSB) + # print "computing the reconstruction noise" + n0Kappa = self.computeQuadEstKappaShearShearBNoiseFFT( + fC0, + fCtot, + fCfg=fCfg, + lMin=lMin, + lMaxS=lMaxS, + lMaxSB=lMaxSB, + corr=corr, + test=test, + cache=cache, + ) + # keep only the real part (the imag. part should be zero, but is tiny in practice) + n0Kappa = np.real(n0Kappa) + # remove the nans + n0Kappa = np.nan_to_num(n0Kappa) + # make sure every value is positive + n0Kappa = np.abs(n0Kappa) + + # fix the issue of the wrong ell=0 value + # replace it by the value lx=0, ly=fundamental + n0Kappa[0, 0] = n0Kappa[0, 1] + + # interpolate + where = (self.ell.flatten() > 0.0) * (self.ell.flatten() < 2.0 * lMax) + L = self.ell.flatten()[where] + N = n0Kappa.flatten()[where] + lnfln = interp1d( + np.log(L), np.log(N), kind="linear", bounds_error=False, fill_value=np.inf + ) + f = lambda ell: np.exp(lnfln(np.log(ell))) + return f diff --git a/soliket/szlike/gnfw.py b/soliket/szlike/gnfw.py new file mode 100644 index 00000000..7f40539c --- /dev/null +++ b/soliket/szlike/gnfw.py @@ -0,0 +1,139 @@ +import numpy as np +from .cosmo import rho_cz + +from ..constants import MSUN_CGS, G_CGS, MPC2CM + +kpc_cgs = MPC2CM * 1.0e-3 + + +def r200(M, z, provider): + """radius of a sphere with density 200 times the critical density of the universe. + Input mass in solar masses. Output radius in cm. + """ + M_cgs = M * MSUN_CGS + om = provider.get_param("Omega_m") + ol = provider.get_param("Omega_L") + Ez2 = om * (1 + z) ** 3 + ol + rhocrit = 1.87847e-29 * provider.get_param("hh") ** 2 + ans = (3 * M_cgs / (4 * np.pi * 200.0 * rhocrit * Ez2)) ** (1.0 / 3.0) + return ans + + +def rho_gnfw1h(x, M, z, theta, provider): + # Bins and weights specific to CMASS + b_cen = np.array([[12.27689266, 12.67884686, 13.16053855, 13.69871423]]).T + p = np.array([4.13431979e-03, 1.31666601e-01, 3.36540698e-01, 8.13760167e-02]) + rho = [] + + fb = provider.get_param("Omega_b") / provider.get_param("Omega_m") + for i in range(0, len(b_cen)): + m = 10 ** b_cen[i] + r200c = r200(m, z, provider) + rvir = r200c / kpc_cgs / 1e3 # Mpc + xc = 0.5 + al = 0.88 * (float(M) / 1e14) ** (-0.03) * (1 + z) ** 0.19 # note change of m to M + gm = -0.2 + rho0, bt = theta + rho.append( + 10**rho0 + * (x / rvir / xc) ** gm + / ((1 + (x / rvir / xc) ** al) ** ((bt - gm) / al)) + * rho_cz(z, provider) + * fb + ) + rho = np.array(rho) + rho_av = np.average(rho, weights=p, axis=0) + return rho_av + + +def rho_gnfw2h(xx, theta2h, twohalo_term): + # Average 2h specific to CMASS, make option for importing file + # or using two_halo script to calculate + rho_file = np.genfromtxt(twohalo_term) + x1 = rho_file[:, 0] + rho2h = rho_file[:, 1] + ans = np.interp(xx, x1, rho2h) + return theta2h * ans + + +def rho_gnfw(xx, M, z, theta, twohalo_term, provider): + theta1h = theta[0], theta[1] + theta2h = theta[2] + ans = rho_gnfw1h(xx, M, z, theta1h, provider) + rho_gnfw2h( + xx, theta2h, twohalo_term + ) + return ans + + +def Pth_gnfw1h(x, M, z, theta, provider): + # Bins and weights specific to CMASS + """ + b_cen = np.array( + [ + [ + 11.31932504, + 11.43785913, + 11.57526319, + 11.74539764, + 11.97016907, + 12.27689266, + 12.67884686, + 13.16053855, + 13.69871423, + ] + ] + ).T + p = np.array( + [ + 2.94467222e-06, + 2.94467222e-06, + 2.94467222e-06, + 1.47233611e-05, + 3.38637305e-05, + 4.13431979e-03, + 1.31666601e-01, + 3.36540698e-01, + 8.13760167e-02, + ] + ) + """ + b_cen = np.array([[12.27689266, 12.67884686, 13.16053855, 13.69871423]]).T + p = np.array([4.13431979e-03, 1.31666601e-01, 3.36540698e-01, 8.13760167e-02]) + + pth = [] + fb = provider.get_param("Omega_b") / provider.get_param("Omega_m") + for i in range(0, len(b_cen)): + m = 10 ** b_cen[i] + r200c = r200(m, z, provider) + rvir = r200c / kpc_cgs / 1e3 # Mpc + M_cgs = m * MSUN_CGS + P200c = G_CGS * M_cgs * 200.0 * rho_cz(z, provider) * fb / (2.0 * r200c) + P0, bt = theta + al = 1.0 + xc = 0.497 * (float(M) / 1e14) ** (-0.00865) * (1 + z) ** 0.731 # note change m-M + gm = -0.3 + pth.append( + P0 * (x / rvir / xc) ** gm * (1 + (x / rvir / xc) ** al) ** (-bt) * P200c + ) + pth = np.array(pth) + pth_av = np.average(pth, weights=p, axis=0) + return pth_av + + +def Pth_gnfw2h(xx, theta2h, twohalo_term): + # 2h specific to CMASS, make option for importing file + # or using two_halo script to calculate + pth_file = np.genfromtxt(twohalo_term) + x1 = pth_file[:, 0] + pth2h = pth_file[:, 2] + ans = np.interp(xx, x1, pth2h) + return theta2h * ans + + +def Pth_gnfw(xx, M, z, theta, twohalo_term, provider): + theta1h = theta[0], theta[1] + theta2h = theta[2] + ans = Pth_gnfw1h(xx, M, z, theta1h, provider) + Pth_gnfw2h( + xx, theta2h, twohalo_term + ) + return ans diff --git a/soliket/szlike/gnfw_test_projections.txt b/soliket/szlike/gnfw_test_projections.txt new file mode 100644 index 00000000..6076c0d4 --- /dev/null +++ b/soliket/szlike/gnfw_test_projections.txt @@ -0,0 +1,16 @@ +# theta (arcmin), kSZ_FFT (muK*sqarcmin), tSZ_FFT (muK*sqarcmin), kSZ_Hankel (muK*sqarcmin), tSZ_Hankel (muK*sqarcmin) +# These are GNFW projections with Battaglia parameters for M=1e13, z=0.55 +# NOTE alpha for kSZ and xc for tSZ are calculated at M and z +# theta_rho=[np.log10(rho0),al,bt,a2h]= [3.1864410705359125, 1.02481552231, 3.45493977635, 1.0] +# theta_pth=[P0,xc,bt,a2h]=[9.10766709141, 0.698455690636, 4.76625102519, 1.0] +# Projections calculated with updated r_ext from basemaps r (with size-res 30-1000), rad = (-3,1) +# Cosmology: omegam = 0.25. Resolution factor (2 for kSZ, 3.5 for tSZ) +1.000000000000000000e+00 6.097523925458924188e-01 -7.832120001071506410e-01 6.154810533013465035e-01 -7.884739859002836138e-01 +1.625000000000000000e+00 1.547954878380804811e+00 -1.899817506621524821e+00 1.558948741006732197e+00 -1.909655669779672493e+00 +2.250000000000000000e+00 2.414207223394933877e+00 -2.798071288943773993e+00 2.419637515397978333e+00 -2.802610721400195626e+00 +2.875000000000000000e+00 3.267267889487348143e+00 -3.602125403857157959e+00 3.275817048117358965e+00 -3.609044610546465304e+00 +3.500000000000000000e+00 4.155779582440343844e+00 -4.390248935673390918e+00 4.164751949877880222e+00 -4.397349380014884623e+00 +4.125000000000000000e+00 5.123799459481552354e+00 -5.229688012569605249e+00 5.134906123788639754e+00 -5.238143105088985507e+00 +4.750000000000000000e+00 6.194313806653685361e+00 -6.151398101568434917e+00 6.206637701471446888e+00 -6.160307472963987507e+00 +5.375000000000000000e+00 7.379992997585159031e+00 -7.172451803705417461e+00 7.397142084626286085e+00 -7.184687401138124763e+00 +6.000000000000000000e+00 8.685564559320285483e+00 -8.299995977161024996e+00 8.705530376792859570e+00 -8.313191742312117327e+00 diff --git a/soliket/szlike/hankel_transform_class.py b/soliket/szlike/hankel_transform_class.py new file mode 100644 index 00000000..af1a3c39 --- /dev/null +++ b/soliket/szlike/hankel_transform_class.py @@ -0,0 +1,101 @@ +# This class is taken from Pixell, written by Sigurd Naess. We don't need all of pixell +# for this, so for now we just take the Hankel transform class. +import numpy as np + + +class RadialFourierTransform: + def __init__(self, lrange=None, rrange=None, n=512, pad=256): + """Construct an object for transforming between radially + symmetric profiles in real-space and fourier space using a + fast Hankel transform. Aside from being fast, this is also + good for representing both cuspy and very extended profiles + due to the logarithmically spaced sample points the fast + Hankel transform uses. A cost of this is that the user can't + freely choose the sample points. Instead one passes the + multipole range or radial range of interest as well as the + number of points to use. + + The function currently assumes two dimensions with flat geometry. + That means the function is only approximate for spherical + geometries, and will only be accurate up to a few degrees + in these cases. + + Arguments: + * lrange = [lmin, lmax]: The multipole range to use. Defaults + to [0.01, 1e6] if no rrange is given. + * rrange = [rmin, rmax]: The radius range to use if lrange is + not specified, in radians. Example values: [1e-7,10]. + Since we don't use spherical geometry r is not limited to 2 pi. + * n: The number of logarithmically equi-spaced points to use + in the given range. Default: 512. The Hankel transform usually + doesn't need many points for good accuracy, and can suffer if + too many points are used. + * pad: How many extra points to pad by on each side of the range. + Padding is useful to get good accuracy in a Hankel transform. + The transforms this function does will return padded output, + which can be unpadded using the unpad method. Default: 256 + """ + if lrange is None and rrange is None: + lrange = [0.1, 1e7] + if lrange is None: + lrange = [1 / rrange[1], 1 / rrange[0]] + logl1, logl2 = np.log(lrange) + logl0 = (logl2 + logl1) / 2 + self.dlog = (logl2 - logl1) / n + i0 = (n + 1) / 2 + pad + self.ell = np.exp(logl0 + (np.arange(1, n + 2 * pad + 1) - i0) * self.dlog) + self.r = 1 / self.ell[::-1] + self.pad = pad + + def real2harm(self, rprof): + """Perform a forward (real -> harmonic) transform, taking us from the + provided real-space radial profile rprof(r) to a harmonic-space profile + lprof(l). rprof can take two forms: + 1. A function rprof(r) that can be called to evalute the profile at + arbitrary points. + 2. An array rprof[self.r] that provides the profile evaluated at the + points given by this object's .r member. + The transform is done along the last axis of the profile. + Returns lprof[self.ell]. This includes padding, which can be removed + using self.unpad""" + import scipy.fft + + try: + rprof = rprof(self.r) + except TypeError: + pass + lprof = 2 * np.pi * scipy.fft.fht(rprof * self.r, self.dlog, 0) / self.ell + return lprof + + def harm2real(self, lprof): + """Perform a backward (harmonic -> real) transform, taking us from the + provided harmonic-space radial profile lprof(l) to a real-space profile + rprof(r). lprof can take two forms: + 1. A function lprof(l) that can be called to evalute the profile at + arbitrary points. + 2. An array lprof[self.ell] that provides the profile evaluated at the + points given by this object's .l member. + The transform is done along the last axis of the profile. + Returns rprof[self.r]. This includes padding, which can be removed + using self.unpad""" + import scipy.fft + + try: + lprof = lprof(self.ell) + except TypeError: + pass + rprof = scipy.fft.ifht(lprof / (2 * np.pi) * self.ell, self.dlog, 0) / self.r + return rprof + + def unpad(self, *arrs): + """Remove the padding from arrays used by this object. The + values in the padded areas of the output of the transform have + unreliable values, but they're not cropped automatically to + allow for round-trip transforms. Example: + r = unpad(r_padded) + r, l, vals = unpad(r_padded, l_padded, vals_padded)""" + if self.pad == 0: + res = arrs + else: + res = tuple([arr[..., self.pad: -self.pad] for arr in arrs]) + return res[0] if len(arrs) == 1 else res diff --git a/soliket/szlike/obb.py b/soliket/szlike/obb.py new file mode 100644 index 00000000..9e16e8cd --- /dev/null +++ b/soliket/szlike/obb.py @@ -0,0 +1,558 @@ +import numpy as np +from scipy.special import spence +from scipy.optimize import fmin +from scipy.optimize import newton + +from ..constants import MSUN_CGS, G_CGS, C_M_S, MPC2CM + +delx = 0.01 +C_CGS = C_M_S * 1.0e2 +kpc_cgs = MPC2CM * 1.0e-3 + + +def nfw(x): + """shape of a NFW profile (NFW 1997, ApJ,490, 493)""" + ans = 1.0 / (x * (1 + x) ** 2) + return ans + + +def gx(x): + ans = np.log(1.0 + x) - x / (1.0 + x) + return ans + + +def gc(c): + ans = 1.0 / (np.log(1.0 + c) - c / (1.0 + c)) + return ans + + +def Hc(c): + ans = ( + -1.0 * np.log(1 + c) / (1.0 + c) + c * (1.0 + 0.5 * c) / ((1.0 + c) ** 2) + ) / gx(c) + return ans + + +def Sc(c): + ans = ( + 0.5 * np.pi**2 + - np.log(c) / 2.0 + - 0.5 / c + - 0.5 / (1 + c) ** 2 + - 3 / (1 + c) + + np.log(1 + c) * (0.5 + 0.5 / c**2 - 2 / c - 1 / (1 + c)) + + 1.5 * (np.log(1 + c)) ** 2 + + 3.0 * spence(c + 1) + ) + + return ans + + +def del_s(c): + ans = Sc(c) / (Sc(c) + (1.0 / c**3) * Hc(c) * gx(c)) + return ans + + +def K_c(c): + ans = 1.0 / 3.0 * Hc(c) / (1.0 - del_s(c)) + return ans + + +def sig_dm2(x, c): + """EQ 14 Lokas & Mamon 2001""" + ans = ( + 0.5 + * x + * c + * gc(c) + * (1 + x) ** 2 + * ( + np.pi**2 + - np.log(x) + - (1.0 / x) + - (1.0 / (1.0 + x) ** 2) + - (6.0 / (1.0 + x)) + + np.log(1.0 + x) * (1.0 + (1.0 / x**2) - 4.0 / x - 2 / (1 + x)) + + 3.0 * (np.log(1.0 + x)) ** 2 + + 6.0 * spence(x + 1) + ) + ) + return ans + + +def r200(M, z, provider): + """radius of a sphere with density 200 times the critical density of the universe. + Input mass in solar masses. Output radius in cm. + """ + M_cgs = M * MSUN_CGS + om = provider.get_param("Omega_m") + ol = provider.get_param("Omega_L") + Ez2 = om * (1 + z) ** 3 + ol + rhocrit = 1.87847e-29 * provider.get_param("hh") ** 2 + ans = (3 * M_cgs / (4 * np.pi * 200.0 * rhocrit * Ez2)) ** (1.0 / 3.0) + return ans + + +def con(M, z): + """ + concentration parameter from Duffy et al. (2008) + input mass in solar masses + """ + ans = 5.71 / (1 + z) ** (0.47) * (M / 2e12) ** (-0.084) + return ans + + +def jx(x, c): + ans = 1.0 - np.log(1.0 + x) / x + ind = np.where(x > c) # [0] + if len(ind) > 0: + ans[ind] = 1.0 - 1.0 / (1.0 + c) - (np.log(1.0 + c) - c / (1.0 + c)) / x[ind] + return ans + + +def jx_f(x, c): + if x <= c: + ans = 1.0 - np.log(1.0 + x) / x + else: + ans = 1.0 - 1.0 / (1.0 + c) - (np.log(1.0 + c) - c / (1.0 + c)) / x + return ans + + +def fx(x, c): + ans = np.log(1.0 + x) / x - 1.0 / (1.0 + c) + ind = np.where(x > c)[0] + if len(ind) > 0: + ans = (np.log(1.0 + c) / c - 1.0 / (1.0 + c)) * c / x + return ans + + +def fstar_func(M): + """Giodini 2009, modified by 0.5""" + ans = 2.5e-2 * (M / (7e13)) ** (-0.37) + return ans + + +def xs_min_func(x, M, z): + c = con(M, z) + fstar = fstar_func(M) + ans = gx(c) * fstar / (1.0 + fstar) - gx(x) + return ans + + +def xs_func(M, z): + x0 = 1.0 + xs = newton( + xs_min_func, + x0, + args=( + M, + z, + ), + ) + return xs + + +def Ks(x_s, M, z): + c = con(M, z) + xx = np.arange(delx / 2.0, x_s, delx) + ans = ( + 1.0 + / gx(c) + * ( + np.sum(Sc(xx) * xx**2) + - 2.0 / 3.0 * np.sum(fx(xx, c) * xx / (1.0 + xx) ** 2) + ) + * delx + ) + return ans + + +def n_exp(gamma): + """exponent of the polytopic e.o.s.""" + ans = 1.0 / (gamma - 1) + return ans + + +def theta_func(x, M, z, theta, theta2): + """polytropic variable""" + gamma, alpha, Ef = theta + beta, x_f = theta2 + c = con(M, z) + nn = n_exp(gamma) + # print(jx(x,c).shape) + ans = 1.0 - beta * jx(x, c) / (1.0 + nn) + return ans + + +def theta_func_f(x, M, z, theta, theta2): + gamma, alpha, Ef = theta + beta, x_f = theta2 + c = con(M, z) + nn = n_exp(gamma) + ans = 1.0 - beta * jx_f(x, c) / (1.0 + nn) + return ans + + +def rho_use(x, M, z, theta, theta2): + gamma, alpha, Ef = theta + beta, x_f = theta2 + nn = n_exp(gamma) + ans = (theta_func(x, M, z, theta, theta2)) ** nn + return ans + + +def rho1h_one_mass(x, M, z, theta, provider): + gamma, alpha, Ef = theta + theta2 = find_params_M(M, z, theta) + P_0, rho_0, x_f = theta2 + nn = n_exp(gamma) + c = con(M, z) + rvir = r200(M, z, provider) + M_cgs = M * MSUN_CGS + beta = rho_0 / P_0 * G_CGS * M_cgs / rvir * c / gx(c) + theta2_use = beta, x_f + ans = ( + rho_0 + * (theta_func(x / (rvir / kpc_cgs / 1e3) * c, M, z, theta, theta2_use)) ** nn + ) + return ans + + +def rho1h(x, M, z, theta, provider): + gamma, alpha, Ef = theta + nn = n_exp(gamma) + # CMASS specific bins and weights + b_cen = np.array([[12.27689266, 12.67884686, 13.16053855, 13.69871423]]).T + p = np.array([4.13431979e-03, 1.31666601e-01, 3.36540698e-01, 8.13760167e-02]) + rh = [] + for i in range(0, len(b_cen)): + m = 10 ** b_cen[i] + thta2 = find_params_M(m, z, theta) + P_0, rho_0, x_f = thta2 + c = con(m, z) + rvir = r200(m, z, provider) + M_cgs = m * MSUN_CGS + beta = rho_0 / P_0 * G_CGS * M_cgs / rvir * c / gx(c) + theta2_use = beta, x_f + rh.append( + rho_0 + * (theta_func(x / (rvir / kpc_cgs / 1e3) * c, M, z, theta, theta2_use)) + ** nn + ) + rh = np.array(rh) + rh_av = np.average(rh, weights=p, axis=0) + return rh_av + + +def rho2h(xx, theta2h): + # CMASS specific 2h + rho_file = np.genfromtxt( + "/home/cemoser/Repositories/SOLikeT/soliket/szlike/twohalo_cmass_average.txt" + ) + x1 = rho_file[:, 0] + rho2h = rho_file[:, 1] + ans = np.interp(xx, x1, rho2h) + return theta2h * ans + + +def rho(x, M, z, theta, provider): + theta1h = theta[0], theta[1], theta[2] + theta2h = theta[3] + ans = rho1h_one_mass(x, M, z, theta1h, provider) + rho2h(x, theta2h) + return ans + + +def rho_outtest(x, M, z, theta, theta2, provider): + gamma, alpha, Ef = theta + P_0, rho_0, x_f = theta2 + nn = n_exp(gamma) + c = con(M, z) + rvir = r200(M, z, provider) + M_cgs = M * MSUN_CGS + beta = rho_0 / P_0 * G_CGS * M_cgs / rvir * c / gx(c) + theta2_use = beta, x_f + ans = rho_0 * (theta_func(x, M, z, theta, theta2_use)) ** nn + return ans + + +def Pnth_th(x, M, z, theta): + gamma, alpha, Ef = theta + c = con(M, z) + ans = 1.0 - alpha * (x / c) ** 0.8 + return ans + + +def Pth1h_one_mass(x, M, z, theta, provider): + gamma, alpha, Ef = theta + theta2 = find_params_M(M, z, theta) + P_0, rho_0, x_f = theta2 + c = con(M, z) + rvir = r200(M, z, provider) + nn = n_exp(gamma) + M_cgs = M * MSUN_CGS + beta = rho_0 / P_0 * G_CGS * M_cgs / rvir * c / gx(c) + theta2_use = beta, x_f + ans = ( + P_0 + * (theta_func(x, M, z, theta, theta2_use)) ** (nn + 1.0) + * Pnth_th(x, M, z, theta) + ) + return ans + + +def Pth1h(x, M, z, theta, provider): + gamma, alpha, Ef = theta + nn = n_exp(gamma) + # CMASS specific bins and weights + b_cen = np.array( + [ + [ + 11.31932504, + 11.43785913, + 11.57526319, + 11.74539764, + 11.97016907, + 12.27689266, + 12.67884686, + 13.16053855, + 13.69871423, + ] + ] + ).T + p = np.array( + [ + 2.94467222e-06, + 2.94467222e-06, + 2.94467222e-06, + 1.47233611e-05, + 3.38637305e-05, + 4.13431979e-03, + 1.31666601e-01, + 3.36540698e-01, + 8.13760167e-02, + ] + ) + pth = [] + for i in range(0, len(b_cen)): + m = 10 ** b_cen[i] + thta2 = find_params_M(m, z, theta) + P_0, rho_0, x_f = thta2 + c = con(m, z) + rvir = r200(m, z, provider) + M_cgs = m * MSUN_CGS + beta = rho_0 / P_0 * G_CGS * M_cgs / rvir * c / gx(c) + theta2_use = beta, x_f + pth.append( + P_0 + * (theta_func(x / (rvir / kpc_cgs / 1e3) * c, m, z, theta, theta2_use)) + ** (nn + 1.0) + * Pnth_th(x / (rvir / kpc_cgs / 1e3) * c, m, z, theta) + ) + pth = np.array(pth) + pth_av = np.average(pth, weights=p, axis=0) + return pth_av + + +def Pth2h(xx, theta2h): + # CMASS specific 2h file + pth_file = np.genfromtxt( + "/home/cemoser/Repositories/SOLikeT/soliket/szlike/twohalo_cmass_average.txt" + ) + x1 = pth_file[:, 0] + pth2h = pth_file[:, 2] + ans = np.interp(xx, x1, pth2h) + return theta2h * ans + + +def Pth(x, M, z, theta, provider): + theta1h = theta[0], theta[1], theta[2] + theta2h = theta[3] + ans = Pth1h(x, M, z, theta1h, provider) + Pth2h(x, theta2h) + return ans + + +def Pth_use(x, M, z, theta, theta2): + gamma, alpha, Ef = theta + beta, x_f = theta2 + nn = n_exp(gamma) + ans = (theta_func(x, M, z, theta, theta2)) ** (nn + 1.0) * Pnth_th(x, M, z, theta) + return ans + + +def Ptot(M, z, theta, theta2, provider): + gamma, alpha, Ef = theta + P_0, rho_0, x_f = theta2 + nn = n_exp(gamma) + rvir = r200(M, z, provider) + c = con(M, z) + M_cgs = M * MSUN_CGS + beta = rho_0 / P_0 * G_CGS * M_cgs / rvir * c / gx(c) + theta2_use = beta, x_f + ans = P_0 * (theta_func_f(x_f, M, z, theta, theta2_use)) ** (nn + 1.0) + return ans + + +def Ptot_use(M, z, theta, theta2): + gamma, alpha, Ef = theta + beta, x_f = theta2 + nn = n_exp(gamma) + ans = (theta_func_f(x_f, M, z, theta, theta2)) ** (nn + 1.0) + return ans + + +def Pnth(x, M, z, theta, theta2, provider): + gamma, alpha, Ef = theta + P_0, rho_0, x_f = theta2 + c = con(M, z) + nn = n_exp(gamma) + rvir = r200(M, z, provider) + beta = ( + rho_0 / P_0 * G_CGS * M / rvir * c / gx(c) + ) # Mvir in original code, but undefined + theta2_use = beta, x_f + ans = ( + alpha + * (x / c) ** 0.8 + * P_0 + * (theta_func(x, M, z, theta, theta2_use)) ** (nn + 1.0) + ) + return ans + + +def Pnth_use(x, M, z, theta, theta2): + gamma, alpha, Ef = theta + c = con(M, z) + nn = n_exp(gamma) + ans = alpha * (x / c) ** 0.8 * (theta_func(x, M, z, theta, theta2)) ** (nn + 1.0) + return ans + + +def I2_int(M, z, theta, theta2): + gamma, alpha, Ef = theta + beta, x_f = theta2 + c = con(M, z) + xx = np.arange(delx / 2.0, x_f, delx) + ans = np.sum(fx(xx, c) * rho_use(xx, M, z, theta, theta2) * xx**2) * delx + return ans + + +def I3_int(M, z, theta, theta2): + gamma, alpha, Ef = theta + beta, x_f = theta2 + xx = np.arange(delx / 2.0, x_f, delx) + ans = np.sum(Pth_use(xx, M, z, theta, theta2) * xx**2) * delx + return ans + + +def I4_int(M, z, theta, theta2): + gamma, alpha, Ef = theta + beta, x_f = theta2 + xx = np.arange(delx / 2.0, x_f, delx) + ans = np.sum(Pnth_use(xx, M, z, theta, theta2) * xx**2) * delx + return ans + + +def L_int(M, z, theta, theta2): + gamma, alpha, Ef = theta + beta, x_f = theta2 + xx = np.arange(delx / 2.0, x_f, delx) + ans = np.sum(rho_use(xx, M, z, theta, theta2) * xx**2) * delx + return ans + + +def rho_0_func(theta0, theta2, provider): + M, z, gamma, alpha, Ef = theta0 + theta = [gamma, alpha, Ef] + c = con(M, z) + rvir = r200(M, z, provider) + fstar = fstar_func(M) + M_cgs = M * MSUN_CGS + fb = provider.get_param("Omega_b") / provider.get_param("Omega_m") + ans = ( + M_cgs + * (fb - fstar) + / (4.0 * np.pi * L_int(M, z, theta, theta2) * (rvir / c) ** 3) + ) + return ans + + +def P_0_func(theta0, theta2, rho_0, provider): + M, z, gamma, alpha, Ef = theta0 + beta, x_f = theta2 + c = con(M, z) + rvir = r200(M, z, provider) + M_cgs = M * MSUN_CGS + ans = rho_0 / beta * G_CGS * M_cgs / rvir * c / gx(c) + return ans + + +def findroots2(theta2, theta0, provider): + M, z, gamma, alpha, Ef = theta0 + theta = [gamma, alpha, Ef] + beta, x_f = theta2 + c = con(M, z) + rvir = r200(M, z, provider) + x_s = xs_func(M, z) + fstar = fstar_func(M) + M_cgs = M * MSUN_CGS + + E_inj = Ef * gx(c) * rvir * fstar / (G_CGS * M_cgs * c) * C_CGS**2 + + Eq1 = ( + 3.0 / 2.0 * (1.0 + fstar) * (K_c(c) * (3.0 - 4.0 * del_s(c)) + Ks(x_s, M, z)) + - E_inj + + 1.0 / 3.0 * (1.0 + fstar) * Sc(c) / gx(c) * (x_f**3 - c**3) + - I2_int(M, z, theta, theta2) / L_int(M, z, theta, theta2) + + 3.0 / 2.0 * I3_int(M, z, theta, theta2) / (beta * L_int(M, z, theta, theta2)) + + 3.0 * I4_int(M, z, theta, theta2) / (beta * L_int(M, z, theta, theta2)) + ) + + Eq2 = (1.0 + fstar) * Sc(c) / gx(c) * ( + beta * L_int(M, z, theta, theta2) + ) - Ptot_use(M, z, theta, theta2) + + ans = Eq1**2 + Eq2**2 + return ans + + +def return_prof_pars(theta2, theta0): + M, z, gamma, alpha, Ef = theta0 + beta, x_f = theta2 + ans = fmin(findroots2, theta2, args=(theta0,), disp=False) + beta_ans, x_f_ans = ans + rho_0 = rho_0_func(theta0, ans) + P_0 = P_0_func(theta0, ans, rho_0) + return P_0, rho_0, x_f_ans + + +def findroots(theta2, theta0, provider): + M, z, gamma, alpha, Ef = theta0 + theta = [gamma, alpha, Ef] + beta, x_f = theta2 + c = con(M, z) + rvir = r200(M, z, provider) + M_cgs = M * MSUN_CGS + + E_inj = Ef * gx(c) * rvir / (G_CGS * M_cgs * c) * C_CGS**2 + + Eq1 = ( + 3.0 / 2.0 * (K_c(c) * (3.0 - 4.0 * del_s(c))) + - E_inj + + 1.0 / 3.0 * Sc(c) / gx(c) * (x_f**3 - c**3) + - I2_int(M, z, theta, theta2) / L_int(M, z, theta, theta2) + + 3.0 / 2.0 * I3_int(M, z, theta, theta2) / (beta * L_int(M, z, theta, theta2)) + + 3.0 * I4_int(M, z, theta, theta2) / (beta * L_int(M, z, theta, theta2)) + ) + Eq2 = Sc(c) / gx(c) * (beta * L_int(M, z, theta, theta2)) - Ptot_use( + M, z, theta, theta2 + ) + return (Eq1, Eq2) + + +def find_params_M(M, z, theta_0): + theta0 = np.append([M, z], [theta_0]) + beta_0 = 1.1 + con_test = con(M, z) + theta2 = np.array([beta_0, con_test * 1.01]) + ans = return_prof_pars(theta2, theta0) + return ans diff --git a/soliket/szlike/projection_functions.py b/soliket/szlike/projection_functions.py new file mode 100644 index 00000000..ca2acc69 --- /dev/null +++ b/soliket/szlike/projection_functions.py @@ -0,0 +1,464 @@ +import numpy as np +from scipy import special +from scipy.interpolate import interp1d +from scipy.signal import convolve +from .cosmo import AngDist +from .gnfw import r200, rho_gnfw1h, Pth_gnfw1h, rho_gnfw, Pth_gnfw +from .obb import con, fstar_func, return_prof_pars, rho, Pth +from .beam import read_beam, f_beam, f_beam_fft, f_response +from .flat_map import FlatMap +from .hankel_transform_class import RadialFourierTransform +from ..constants import ( + MPC2CM, + C_M_S, + h_Planck, + k_Boltzmann, + electron_mass_kg, + proton_mass_kg, + hydrogen_fraction, + T_CMB, + ST_CGS, +) + +kpc_cgs = MPC2CM * 1.0e-3 +C_CGS = C_M_S * 1.0e2 +ME_CGS = electron_mass_kg * 1.0e3 +MP_CGS = proton_mass_kg * 1.0e3 +sr2sqarcmin = 3282.8 * 60.0**2 +XH = hydrogen_fraction + + +def radius_definition(transform_type): + # There is probably a cleaner way to do this... + if transform_type == "FFT": + # New default from 15-500 (Amodeo21) to 30-1000 (Moser23) + sizeArcmin = 30.0 # [degree] + sz = 1000 + baseMap = FlatMap( + nX=sz, + nY=sz, + sizeX=sizeArcmin * np.pi / 180.0 / 60.0, + sizeY=sizeArcmin * np.pi / 180.0 / 60.0, + ) + # map of the cluster, not yet convolved with the beam + xc = baseMap.sizeX / 2.0 + yc = baseMap.sizeY / 2.0 + r = np.sqrt( + (baseMap.x - xc) ** 2 + (baseMap.y - yc) ** 2 + ) # cluster centric radius in rad + r_max = np.max(r) + + elif transform_type == "Hankel": + baseMap = False + r = False + sizeArcmin = 30.0 + r_max = sizeArcmin * np.pi / 180.0 / 60.0 / np.sqrt(2) + return baseMap, r, r_max + + +def coth(x): + return 1 / np.tanh(x) + + +def fnu(nu): + """input frequency in GHz""" + nu *= 1e9 + x = h_Planck * nu / (k_Boltzmann * T_CMB) + ans = x * coth(x / 2.0) - 4.0 + return ans + + +def convolve_FFT( + r, thta_smooth, prof2D, beam_txt, baseMap, thta_use, beam_response=False +): + profMap = np.interp(r, thta_smooth, prof2D) + beamMapF = f_beam_fft(beam_txt, baseMap.ell) + # Fourier transform the profile + profMapF = baseMap.fourier(profMap) # see note in flat_map.py about precision + # multiply by the beam transfer function + convolvedProfMapF = profMapF * beamMapF + + if beam_response is not False: # for tSZ + respTF = f_response(beam_response, baseMap.ell) + # multiply by the response + convolvedProfMapF *= respTF + + # inverse fourier transform + convolvedProfMap = baseMap.inverseFourier(convolvedProfMapF) + prof2D_beam = interp1d( + r.flatten(), + convolvedProfMap.flatten(), + kind="linear", + bounds_error=False, + fill_value=0.0, + )(thta_use) + return prof2D_beam + + +def convolve_Hankel(thta_smooth, prof2D, beam_txt, thta_use, beam_response=False): + rht = RadialFourierTransform( + n=200, pad=100, lrange=[170.0, 1.4e6] + ) # note hard values, n here needs to be same size as rad in project + profMap = np.interp(rht.r, thta_smooth, prof2D) + lprofs = rht.real2harm(profMap) + lprofs *= f_beam_fft(beam_txt, rht.ell) + + if beam_response is not False: + respTF = f_response(beam_response, rht.ell) + # multiply by the response + lprofs *= respTF + + rprofs = rht.harm2real(lprofs) + # padding + r_unpad, rprofs = rht.unpad(rht.r, rprofs) + prof2D_beam = interp1d( + r_unpad.flatten(), + rprofs.flatten(), + kind="linear", + bounds_error=False, + fill_value=0.0, + )(thta_use) + return prof2D_beam + + +def project_ksz( + tht, M, z, beam_txt, transform_type, model_params, twohalo_term, provider +): # input_model + disc_fac = np.sqrt(2) + NNR = 100 + resolution_factor = 2.0 + NNR2 = resolution_factor * NNR + AngDis = AngDist(z, provider) + + baseMap, r, r_max = radius_definition(transform_type) + + r_use = AngDis * np.arctan(np.radians(tht / 60.0)) + r_use2 = AngDis * np.arctan(np.radians(tht * disc_fac / 60.0)) + r_ext = AngDis * np.arctan(r_max) # total profile + r_ext2 = r_ext + + rad = np.logspace(-3, 1, 200) # Mpc + rad2 = rad + + radlim = r_ext + radlim2 = r_ext2 + + dtht = np.arctan(radlim / AngDis) / NNR # rads + dtht2 = np.arctan(radlim2 / AngDis) / NNR # rads + dtht_use = np.arctan(r_use / AngDis) / NNR + dtht2_use = np.arctan(r_use2 / AngDis) / NNR + + thta_use = (np.arange(NNR) + 1.0) * dtht_use + thta2_use = (np.arange(NNR) + 1.0) * dtht2_use + + thta_smooth = (np.arange(NNR2) + 1.0) * dtht / resolution_factor + thta2_smooth = (np.arange(NNR2) + 1.0) * dtht2 / resolution_factor + thta_smooth = thta_smooth[:, None] + thta2_smooth = thta2_smooth[:, None] + + rint = np.sqrt(rad**2 + thta_smooth**2 * AngDis**2) + rint2 = np.sqrt(rad2**2 + thta2_smooth**2 * AngDis**2) + + # choose the model for projection + # models = { + # "gnfw":rho_gnfw, + # "obb":rho + # } + # chosen_model = models.get(input_model) + + rho2D = ( + 2 + * np.trapz( + rho_gnfw(rint, M, z, model_params, twohalo_term, provider), + x=rad * kpc_cgs, + axis=1, + ) + * 1e3 + ) + rho2D2 = ( + 2 + * np.trapz( + rho_gnfw(rint2, M, z, model_params, twohalo_term, provider), + x=rad2 * kpc_cgs, + axis=1, + ) + * 1e3 + ) + + thta_smooth = (np.arange(NNR2) + 1.0) * dtht / resolution_factor + thta2_smooth = (np.arange(NNR2) + 1.0) * dtht2 / resolution_factor + + if transform_type == "FFT": + rho2D_beam = convolve_FFT(r, thta_smooth, rho2D, beam_txt, baseMap, thta_use) + rho2D2_beam = convolve_FFT( + r, thta2_smooth, rho2D2, beam_txt, baseMap, thta2_use + ) + + elif transform_type == "Hankel": + rho2D_beam = convolve_Hankel(thta_smooth, rho2D, beam_txt, thta_use) + rho2D2_beam = convolve_Hankel(thta2_smooth, rho2D2, beam_txt, thta2_use) + + sig = 2.0 * np.pi * dtht_use * np.sum(thta_use * rho2D_beam) + sig2 = 2.0 * np.pi * dtht2_use * np.sum(thta2_use * rho2D2_beam) + + sig_all_beam = ( + (2 * sig - sig2) + * provider.get_param("v_rms") + * ST_CGS + * T_CMB + * 1e6 + * ((1.0 + XH) / 2) + / MP_CGS + ) # units in muK*sr + # sig_all_beam *= sr2sqarcmin #units in muK*sqarcmin + return sig_all_beam + + +def project_tsz( + tht, + M, + z, + nu, + beam_txt, + transform_type, + model_params, + beam_response, + twohalo_term, + provider, +): + disc_fac = np.sqrt(2) + NNR = 100 + resolution_factor = 3.5 + NNR2 = resolution_factor * NNR + AngDis = AngDist(z, provider) + + baseMap, r, r_max = radius_definition(transform_type) + + r_use = AngDis * np.arctan(np.radians(tht / 60.0)) + r_use2 = AngDis * np.arctan(np.radians(tht * disc_fac / 60.0)) + r_ext = AngDis * np.arctan(r_max) # total profile + r_ext2 = r_ext + + rad = np.logspace(-3, 1, 200) # Mpc + rad2 = rad + + radlim = r_ext + radlim2 = r_ext2 + + dtht = np.arctan(radlim / AngDis) / NNR # rads + dtht2 = np.arctan(radlim2 / AngDis) / NNR # rads + dtht_use = np.arctan(r_use / AngDis) / NNR + dtht2_use = np.arctan(r_use2 / AngDis) / NNR + + thta_use = (np.arange(NNR) + 1.0) * dtht_use + thta2_use = (np.arange(NNR) + 1.0) * dtht2_use + + thta_smooth = (np.arange(NNR2) + 1.0) * dtht / resolution_factor + thta2_smooth = (np.arange(NNR2) + 1.0) * dtht2 / resolution_factor + thta_smooth = thta_smooth[:, None] + thta2_smooth = thta2_smooth[:, None] + + rint = np.sqrt(rad**2 + thta_smooth**2 * AngDis**2) + rint2 = np.sqrt(rad2**2 + thta2_smooth**2 * AngDis**2) + + Pth2D = ( + 2 + * np.trapz( + Pth_gnfw(rint, M, z, model_params, twohalo_term, provider), + x=rad * kpc_cgs, + axis=1, + ) + * 1e3 + ) + Pth2D2 = ( + 2 + * np.trapz( + Pth_gnfw(rint2, M, z, model_params, twohalo_term, provider), + x=rad2 * kpc_cgs, + axis=1, + ) + * 1e3 + ) + + thta_smooth = (np.arange(NNR2) + 1.0) * dtht / resolution_factor + thta2_smooth = (np.arange(NNR2) + 1.0) * dtht2 / resolution_factor + + if transform_type == "FFT": + Pth2D_beam = convolve_FFT( + r, thta_smooth, Pth2D, beam_txt, baseMap, thta_use, beam_response + ) + Pth2D2_beam = convolve_FFT( + r, thta2_smooth, Pth2D2, beam_txt, baseMap, thta2_use, beam_response + ) + elif transform_type == "Hankel": + Pth2D_beam = convolve_Hankel( + thta_smooth, Pth2D, beam_txt, thta_use, beam_response + ) + Pth2D2_beam = convolve_Hankel( + thta2_smooth, Pth2D2, beam_txt, thta2_use, beam_response + ) + + sig_p = 2.0 * np.pi * dtht_use * np.sum(thta_use * Pth2D_beam) + sig2_p = 2.0 * np.pi * dtht2_use * np.sum(thta2_use * Pth2D2_beam) + sig_all_p_beam = ( + (2 * sig_p - sig2_p) + * ST_CGS + / (ME_CGS * C_CGS**2) + * ((2.0 + 2.0 * XH) / (3.0 + 5.0 * XH)) + * 1e6 + ) + + # sig_all_p_beam *= sr2sqarcmin #units in muK*sqarcmin + return sig_all_p_beam + + +def project_obb(tht, M, z, beam_txt, theta, nu, fbeam, provider): + # NOTE: this is convolving through analytical integral, not FFT or FHT + # NOTE: needs to be updated with new method for r from baseMap and Hankel transform + disc_fac = np.sqrt(2) + NNR = 100 + resolution_factor = 3.0 + NNR2 = resolution_factor * NNR + + theta_r = [theta[0], theta[1], theta[2], theta[3]] + theta_p = [theta[0], theta[1], theta[2], theta[4]] + + AngDis = AngDist(z) + + r_ext = AngDis * np.arctan(np.radians(tht / 60.0)) + r_ext2 = AngDis * np.arctan(np.radians(tht * disc_fac / 60.0)) + + rad = np.logspace(-3, 1, 200) + rad2 = np.logspace(-3, 1, 200) + + radlim = r_ext + radlim2 = r_ext2 + + dtht = np.arctan(radlim / AngDis) / NNR + dtht2 = np.arctan(radlim2 / AngDis) / NNR + + thta = (np.arange(NNR) + 1.0) * dtht + thta2 = (np.arange(NNR) + 1.0) * dtht2 + + thta_smooth = (np.arange(NNR2) + 1.0) * dtht / resolution_factor + thta2_smooth = (np.arange(NNR2) + 1.0) * dtht2 / resolution_factor + + thta_smooth = thta_smooth[:, None] + thta2_smooth = thta2_smooth[:, None] + + rint = np.sqrt(rad**2 + thta_smooth**2 * AngDis**2) + rint2 = np.sqrt(rad2**2 + thta2_smooth**2 * AngDis**2) + + rho2D = 2 * np.trapz(rho(rint, M, z, theta_r), x=rad * kpc_cgs, axis=1) * 1e3 + rho2D2 = 2 * np.trapz(rho(rint2, M, z, theta_r), x=rad2 * kpc_cgs, axis=1) * 1e3 + Pth2D = 2 * np.trapz(Pth(rint, M, z, theta_p), x=rad * kpc_cgs, axis=1) * 1e3 + Pth2D2 = 2 * np.trapz(Pth(rint2, M, z, theta_p), x=rad2 * kpc_cgs, axis=1) * 1e3 + + thta_smooth = (np.arange(NNR2) + 1.0) * dtht / resolution_factor + thta = thta[:, None, None] + thta2_smooth = (np.arange(NNR2) + 1.0) * dtht2 / resolution_factor + thta2 = thta2[:, None, None] + + phi = np.linspace(0.0, 2 * np.pi, 100) + phi = phi[None, None, :] + thta_smooth = thta_smooth[None, :, None] + thta2_smooth = thta2_smooth[None, :, None] + + rho2D = rho2D[None, :, None] + rho2D2 = rho2D2[None, :, None] + Pth2D = Pth2D[None, :, None] + Pth2D2 = Pth2D2[None, :, None] + + b = read_beam(beam_txt) + + rho2D_beam0 = np.trapz( + thta_smooth + * rho2D + * f_beam( + np.sqrt( + thta**2 + thta_smooth**2 - 2 * thta * thta_smooth * np.cos(phi) + ), + b, + ), + x=phi, + axis=2, + ) + + rho2D2_beam0 = np.trapz( + thta2_smooth + * rho2D2 + * f_beam( + np.sqrt( + thta2**2 + thta2_smooth**2 - 2 * thta2 * thta2_smooth * np.cos(phi) + ), + b, + ), + x=phi, + axis=2, + ) + + Pth2D_beam0 = np.trapz( + thta_smooth + * Pth2D + * f_beam( + np.sqrt( + thta**2 + thta_smooth**2 - 2 * thta * thta_smooth * np.cos(phi) + ), + b, + ), + x=phi, + axis=2, + ) + + Pth2D2_beam0 = np.trapz( + thta2_smooth + * Pth2D2 + * f_beam( + np.sqrt( + thta2**2 + thta2_smooth**2 - 2 * thta2 * thta2_smooth * np.cos(phi) + ), + b, + ), + x=phi, + axis=2, + ) + + thta_smooth = (np.arange(NNR2) + 1.0) * dtht / resolution_factor + thta2_smooth = (np.arange(NNR2) + 1.0) * dtht2 / resolution_factor + + rho2D_beam = np.trapz(rho2D_beam0, x=thta_smooth, axis=1) + rho2D2_beam = np.trapz(rho2D2_beam0, x=thta2_smooth, axis=1) + Pth2D_beam = np.trapz(Pth2D_beam0, x=thta_smooth, axis=1) + Pth2D2_beam = np.trapz(Pth2D2_beam0, x=thta2_smooth, axis=1) + + thta = (np.arange(NNR) + 1.0) * dtht + thta2 = (np.arange(NNR) + 1.0) * dtht2 + + sig = 2.0 * np.pi * dtht * np.sum(thta * rho2D_beam) + sig2 = 2.0 * np.pi * dtht2 * np.sum(thta2 * rho2D2_beam) + + sig_all_beam = ( + (2 * sig - sig2) + * provider.get_param("v_rms") + * ST_CGS + * T_CMB + * 1e6 + * ((1.0 + XH) / 2.0) + / MP_CGS + ) + + sig_p = 2.0 * np.pi * dtht * np.sum(thta * Pth2D_beam) + sig2_p = 2.0 * np.pi * dtht2 * np.sum(thta2 * Pth2D2_beam) + + sig_all_p_beam = ( + fnu(nu) + * (2 * sig_p - sig2_p) + * ST_CGS + / (ME_CGS * C_CGS**2) + * T_CMB + * 1e6 + * ((2.0 + 2.0 * XH) / (3.0 + 5.0 * XH)) + ) + # sig_all_beam *= sr2sqarcmin + # sig_all_p_beam *= sr2sqarcmin + return sig_all_beam, sig_all_p_beam diff --git a/soliket/szlike/szlike.py b/soliket/szlike/szlike.py new file mode 100644 index 00000000..3be4df2d --- /dev/null +++ b/soliket/szlike/szlike.py @@ -0,0 +1,139 @@ +""" +Likelihood for SZ model +""" +import numpy as np +from ..gaussian import GaussianData, GaussianLikelihood +from .projection_functions import project_ksz, project_tsz, project_obb + + +class SZLikelihood(GaussianLikelihood): + def initialize(self): + + self.beam_txt = self.beam_file + self.transform_type = self.transform_type + self.z = self.redshift + self.nu = self.frequency_GHz + self.M = self.mass_halo_mean_Msol + # self.input_model = self.input_model + self.beam_response = self.beam_response + self.twohalo_term = self.twohalo_term + + x, y, cov = self._get_data() + self.data = GaussianData("SZModel", x, y, cov) + + def logp(self, **params_values): + theory = self._get_theory(**params_values) + return self.data.loglike(theory) + + +class KSZLikelihood(SZLikelihood): # this is for GNFW model + def _get_data(self, **params_values): + # Modify these lines as necessary. Usecols 1 for FFT, 3 for Hankel + thta_arc, ksz_data = np.loadtxt(self.sz_data_file, usecols=(0, 3), unpack=True) + ksz_data /= 3282.8 * 60.0**2 # units muK*sr + cov_ksz = np.loadtxt(self.cov_ksz_file) # units muK*sr + + self.thta_arc = thta_arc + self.ksz_data = ksz_data + self.cov = cov_ksz + return self.thta_arc, self.ksz_data, self.cov + + def _get_theory(self, **params_values): + model_params = [ + params_values["gnfw_rho0"], + params_values["gnfw_bt_ksz"], + params_values["gnfw_A2h_ksz"], + ] + # model_params = model_params.get(self.input_model) + + rho = np.zeros(len(self.thta_arc)) + for ii in range(len(self.thta_arc)): + rho[ii] = project_ksz( + self.thta_arc[ii], + self.M, + self.z, + self.beam_txt, + self.transform_type, + model_params, + self.twohalo_term, + self.provider, + ) # self.input_model + return rho + + +class TSZLikelihood(SZLikelihood): # this is for GNFW model + def _get_data(self, **params_values): + # Modify these lines as necessary. Usecols 2 for FFT, 4 for Hankel + thta_arc, tsz_data = np.loadtxt(self.sz_data_file, usecols=(0, 4), unpack=True) + tsz_data /= 3282.8 * 60.0**2 # units muK*sr + cov_tsz = np.loadtxt(self.cov_tsz_file) # units muK*sr + + self.thta_arc = thta_arc + self.tsz_data = tsz_data + self.cov = cov_tsz + return self.thta_arc, self.tsz_data, self.cov + + def _get_theory(self, **params_values): + model_params = [ + params_values["gnfw_P0"], + params_values["gnfw_bt_tsz"], + params_values["gnfw_A2h_tsz"], + ] + + pth = np.zeros(len(self.thta_arc)) + for ii in range(len(self.thta_arc)): + pth[ii] = project_tsz( + self.thta_arc[ii], + self.M, + self.z, + self.nu, + self.beam_txt, + self.transform_type, + model_params, + self.beam_response, + self.twohalo_term, + self.provider, + ) + return pth + + +class OBBLikelihood(SZLikelihood): # OBB model, tSZ and kSZ together + def _get_data(self, **params_values): + thta_arc, ksz_data, tsz_data = np.loadtxt( + self.sz_data_file, usecols=(0, 1, 2), unpack=True + ) + ksz_data /= 3282.8 * 60.0**2 + tsz_data /= 3282.8 * 60.0**2 + cov_ksz = np.loadtxt(self.cov_ksz_file) # different cov file? combined? + cov_tsz = np.loadtxt(self.cov_tsz_file) + + self.thta_arc = thta_arc + self.ksz_data = ksz_data + self.tsz_data = tsz_data + self.cov_ksz = cov_ksz + self.cov_tsz = cov_tsz + return self.thta_arc, self.ksz_data, self.tsz_data, self.cov_ksz, self.cov_tsz + + def _get_theory(self, **params_values): + model_params = [ + params_values["obb_Gamma"], + params_values["obb_alpha_Nth"], + params_values["obb_logE"], + params_values["obb_Ak2h"], + params_values["obb_At2h"], + ] + + rho = np.zeros(len(self.thta_arc)) + pth = np.zeros(len(self.thta_arc)) + for ii in range(len(self.thta_arc)): + temp = project_obb( + self.thta_arc[ii], + self.M, + self.z, + self.beam_txt, + model_params, + self.provider, + ) + rho[ii] = temp[0] + pth[ii] = temp[1] + return rho, pth diff --git a/soliket/szlike/two_halo.py b/soliket/szlike/two_halo.py new file mode 100644 index 00000000..4060be4b --- /dev/null +++ b/soliket/szlike/two_halo.py @@ -0,0 +1,221 @@ +import numpy as np +from scipy.integrate import quad +from .cosmo import rho_cz +from .gnfw import r200 +from hmf import MassFunction, transfer +from colossus.lss import bias, peaks +from colossus.cosmology import cosmology + + +"""set cosmology""" +params = { + "flat": True, + "H0": 70.0, + "Om0": 0.25, + "Ob0": 0.044, + "sigma8": 0.8159, + "ns": 0.97, +} + +"""parameters used in Illustris TNG https://arxiv.org/pdf/1703.02970.pdf""" +# params = { +# 'flat': True, +# 'H0': 67.7, +# 'Om0': 0.31, +# 'Ob0': 0.0486, +# 'sigma8': 0.8159, +# 'ns': 0.97, +# } + +cosmology.setCosmology("myCosmo", params) + +hh = 0.7 +fb = params["Ob0"] / params["Om0"] +G_cgs = 6.67259e-8 # cm3/g/s2 +rhocrit = 1.87847e-29 * hh**2 +Msol_cgs = 1.989e33 +kpc_cgs = 3.086e21 + +################################################################################# +# Computing the 2-halo component of density and preassure profiles in cgs units # +################################################################################# + + +# From Battaglia 2016, Appendix A +def rho_gnfw(x, m, z): + rho0 = 4e3 * (m / 1e14) ** 0.29 * (1 + z) ** (-0.66) + al = 0.88 * (m / 1e14) ** (-0.03) * (1 + z) ** 0.19 + bt = 3.83 * (m / 1e14) ** 0.04 * (1 + z) ** (-0.025) + xc = 0.5 + gm = -0.2 + ans = rho0 * (x / xc) ** gm * (1 + (x / xc) ** al) ** (-(bt - gm) / al) + ans *= rho_cz(z) * fb + return ans + + +def rhoFourier(k, m, z): + ans = [] + for i in range(len(m)): + r200c = r200(m[i], z) / kpc_cgs / 1e3 + integrand = ( + lambda r: 4.0 + * np.pi + * r**2 + * rho_gnfw(r / r200c, m[i], z) + * np.sin(k * r) + / (k * r) + ) + res = quad(integrand, 0.0, 10 * r200c, epsabs=0.0, epsrel=1.0e-4, limit=10000)[ + 0 + ] + ans.append(res) + ans = np.array(ans) + return ans + + +def hmf(m, z): + """Shet, Mo & Tormen 2001""" + Mmin = np.log10(np.min(m)) + Mmax = np.log10(np.max(m)) + return MassFunction( + z=z, Mmin=Mmin, Mmax=Mmax, dlog10m=(Mmax - Mmin) / 49.5, hmf_model="SMT" + ).dndm + + +def b(m, z): + """Shet, Mo & Tormen 2001""" + nu = peaks.peakHeight(m, z) + delta_c = peaks.collapseOverdensity(corrections=True, z=z) + aa, bb, cc = 0.707, 0.5, 0.6 + return 1.0 + 1.0 / (np.sqrt(aa) * delta_c) * ( + np.sqrt(aa) * aa * nu**2 + + np.sqrt(aa) * bb * (aa * nu**2) ** (1 - cc) + - ((aa * nu**2) ** cc / (aa * nu**2) ** cc + bb * (1 - cc) * (1 - cc / 2)) + ) + + +def Plin(k, z): + lnk_min = np.log(np.min(k)) + lnk_max = np.log(np.max(k)) + dlnk = (lnk_max - lnk_min) / (49.5) + """Eisenstein & Hu (1998)""" + p = transfer.Transfer( + sigma_8=0.8344, + n=0.9624, + z=z, + lnk_min=lnk_min, + lnk_max=lnk_max, + dlnk=dlnk, + transfer_model="EH", + ) + return p.power + + +def rho_2h(r, m, z): + + # first compute P_2h (power spectrum) + m_array = np.logspace(np.log10(1.0e10), np.log10(1.0e15), 50, 10.0) + k_array = np.logspace(np.log10(1.0e-3), np.log10(1.0e3), 50, 10.0) + hmf_array = np.array([hmf(m_array, z)] * len(k_array)).reshape( + len(k_array), len(hmf(m_array, z)) + ) + bias_array = np.array([b(m_array, z)] * len(k_array)).reshape( + len(k_array), len(b(m_array, z)) + ) + + arr = [] + for i in range(len(k_array)): + arr.append( + np.trapz( + hmf_array[i, :] * bias_array[i, :] * rhoFourier(k_array[i], m_array, z), + m_array, + ) + ) + arr = np.array(arr) + P2h = np.array(arr * b(m, z) * Plin(k_array, z)) + + # then Fourier transform to get rho_2h + rcorr = 50.0 # Mpc/h + integrand = ( + lambda k: 1.0 + / (2 * np.pi**2.0) + * k**2 + * np.sin(k * r) + / (k * r) + * np.interp(k, k_array, P2h) + if k > 1.0 / rcorr + else 0.0 + ) + res = quad(integrand, 0.0, np.inf, epsabs=0.0, epsrel=1.0e-2, limit=1000)[0] + return res + + +# From Battaglia 2012, AGN Feedback Delta=200 +def Pth_gnfw(x, m, z): + P200c = G_cgs * m * Msol_cgs * 200.0 * rho_cz(z) * fb / (2.0 * r200(m, z)) + P0 = 18.1 * (m / 1e14) ** 0.154 * (1 + z) ** (-0.758) + al = 1.0 + bt = 4.35 * (m / 1e14) ** 0.0393 * (1 + z) ** 0.415 + xc = 0.497 * (m / 1e14) ** (-0.00865) * (1 + z) ** 0.731 + gm = -0.3 + ans = P0 * (x / xc) ** gm * (1 + (x / xc) ** al) ** (-bt) + ans *= P200c + return ans + + +def PthFourier(k, m, z): + ans = [] + for i in range(len(m)): + r200c = r200(m[i], z) / kpc_cgs / 1e3 + integrand = ( + lambda r: 4.0 + * np.pi + * r**2 + * Pth_gnfw(r / r200c, m[i], z) + * np.sin(k * r) + / (k * r) + ) + res = quad(integrand, 0.0, 10 * r200c, epsabs=0.0, epsrel=1.0e-4, limit=10000)[ + 0 + ] + ans.append(res) + ans = np.array(ans) + return ans + + +def Pth_2h(r, m, z): + # first compute P_2h (power spectrum) + m_array = np.logspace(np.log10(1.0e10), np.log10(1.0e15), 50, 10.0) + k_array = np.logspace(np.log10(1.0e-3), np.log10(1.0e3), 50, 10.0) + hmf_array = np.array([hmf(m_array, z)] * len(k_array)).reshape( + len(k_array), len(hmf(m_array, z)) + ) + bias_array = np.array([b(m_array, z)] * len(k_array)).reshape( + len(k_array), len(b(m_array, z)) + ) + + arr = [] + for i in range(len(k_array)): + arr.append( + np.trapz( + hmf_array[i, :] * bias_array[i, :] * PthFourier(k_array[i], m_array, z), + m_array, + ) + ) + arr = np.array(arr) + P2h = np.array(arr * b(m, z) * Plin(k_array, z)) + + # then Fourier transform to get Pth_2h + rcorr = 50.0 # Mpc/h + integrand = ( + lambda k: 1.0 + / (2 * np.pi**2.0) + * k**2 + * np.sin(k * r) + / (k * r) + * np.interp(k, k_array, P2h) + if k > 1.0 / rcorr + else 0.0 + ) + res = quad(integrand, 0.0, np.inf, epsabs=0.0, epsrel=1.0e-2, limit=1000)[0] + return res diff --git a/soliket/szlike/twohalo_cmass_average.txt b/soliket/szlike/twohalo_cmass_average.txt new file mode 100644 index 00000000..efb18b14 --- /dev/null +++ b/soliket/szlike/twohalo_cmass_average.txt @@ -0,0 +1,51 @@ +# R[Mpc/h] rho2h[g/cm3] Pth2h[g cm^-1 s^-2] +4.999999999999999237e-03 2.618270934242747938e-29 2.751490979892967716e-14 +5.838993118688259529e-03 2.617383733832416308e-29 2.749554969891851205e-14 +6.818768128017769353e-03 2.616394364464526138e-29 2.747341678347436946e-14 +7.962948035485325690e-03 2.614729847282972680e-29 2.744741507546696612e-14 +9.299119756734202866e-03 2.612742997193979042e-29 2.741399580629183163e-14 +1.085949925388581157e-02 2.610148224330089489e-29 2.736782781891133806e-14 +1.268170828316790913e-02 2.607321110310440424e-29 2.731645739095418349e-14 +1.480968147972586656e-02 2.603586669167293118e-29 2.724806351036433444e-14 +1.729472565001686016e-02 2.598714872296284476e-29 2.716307332677212162e-14 +2.019675681200996942e-02 2.592466273056603281e-29 2.706000585811806428e-14 +2.358574480902929024e-02 2.584806579527857772e-29 2.692004652265542090e-14 +2.754340032781188788e-02 2.575182911247127929e-29 2.674763578770746694e-14 +3.216514499587391773e-02 2.563347915176872614e-29 2.654120016556515949e-14 +3.756241205850358728e-02 2.549073885302699444e-29 2.628719952876999402e-14 +4.386533310618707243e-02 2.531580039566457857e-29 2.598453540984018233e-14 +5.122587563119895132e-02 2.510178316813152019e-29 2.562497544863164912e-14 +5.982150706187026185e-02 2.483966106988028524e-29 2.520284692708922264e-14 +6.985947361676432543e-02 2.452725155357683242e-29 2.471690045820780001e-14 +8.158179714469422639e-02 2.416377331657314424e-29 2.418123821017286189e-14 +9.527111042761822468e-02 2.373359060543518973e-29 2.353042319502274649e-14 +1.112574716393304325e-01 2.324564651281875056e-29 2.280637819578222734e-14 +1.299263222609410018e-01 2.268564873729186883e-29 2.202705273781179367e-14 +1.517277803236215528e-01 2.206217957048111193e-29 2.118276624497907384e-14 +1.771874930446940555e-01 2.136511223039484019e-29 2.028767490298870059e-14 +2.069193105211186157e-01 2.062098417641808487e-29 1.928924866550514023e-14 +2.416400860513061855e-01 1.976714177859083127e-29 1.825975692611017153e-14 +2.821869599305631415e-01 1.885493787164629467e-29 1.716348972683359832e-14 +3.295375434436236262e-01 1.789421230905669596e-29 1.604476268280970278e-14 +3.848334897033505730e-01 1.688181542329838820e-29 1.489986535587761221e-14 +4.494080196437306718e-01 1.581790287793085824e-29 1.374107168462487489e-14 +5.248180668366125845e-01 1.472336496194130193e-29 1.258067825694042138e-14 +6.128818161644512541e-01 1.361310943265311154e-29 1.144019516313607165e-14 +7.157225414306788158e-01 1.249077249288084029e-29 1.032925944111152182e-14 +8.358197988607614093e-01 1.137464909667880694e-29 9.251966611485096286e-15 +9.760692108022782643e-01 1.027606070764040262e-29 8.227280563030819913e-15 +1.139852281047596616e+00 9.212346267869875192e-30 7.269240224986001390e-15 +1.331117925071608088e+00 8.189669712163945643e-30 6.360424797785021405e-15 +1.554477680931142913e+00 7.218328515111470686e-30 5.536101749394299826e-15 +1.815316896422285664e+00 6.308114584067706573e-30 4.774224007802583545e-15 +2.119924573289651004e+00 5.469180596647709917e-30 4.089881947029861928e-15 +2.475644999115283795e+00 4.689615038357133966e-30 3.471546340910047734e-15 +2.891054822829829174e+00 3.980663869180243879e-30 2.927149441266113813e-15 +3.376169843250778957e+00 3.354875454348986625e-30 2.443075132559373557e-15 +3.942686496452824318e+00 2.791231382619242046e-30 2.019357333523417481e-15 +4.604263864386632932e+00 2.306812732981102104e-30 1.660170082722305835e-15 +5.376853004155713656e+00 1.872691355954030531e-30 1.341104704988826752e-15 +6.279081538292701659e+00 1.514658284834591463e-30 1.082700622540113486e-15 +7.332702778754716100e+00 1.196549718849060260e-30 8.528035872974859142e-16 +8.563120213307023221e+00 9.472798710328329015e-31 6.724522766896366154e-16 +1.000000000000000000e+01 7.214389710425451540e-31 5.131766200020638630e-16 diff --git a/tests/test_szlike.py b/tests/test_szlike.py new file mode 100644 index 00000000..161ef618 --- /dev/null +++ b/tests/test_szlike.py @@ -0,0 +1,84 @@ +import numpy as np +import os +from soliket.szlike import KSZLikelihood, TSZLikelihood +from cobaya.model import get_model + +sz_data_file = "soliket/szlike/gnfw_test_projections.txt" +cov_tsz_file = "soliket/szlike/cov_diskring_tsz_varweight_bootstrap.txt" +cov_ksz_file = "soliket/szlike/cov_diskring_ksz_varweight_bootstrap.txt" +twohalo_term = "soliket/szlike/twohalo_cmass_average.txt" +beam_file = "soliket/szlike/beam_f150_daynight.txt" +beam_response = "soliket/szlike/act_planck_s08_s18_cmb_f150_daynight_response_tsz.txt" + + +def test_ksz(request): + info = { + "params": { + "gnfw_rho0": 3.1864410705359125, + "gnfw_bt_ksz": 3.45493977635, + "gnfw_A2h_ksz": 1.0, + }, + "likelihood": { + "KSZLikelihood": { + "external": KSZLikelihood, + "sz_data_file": os.path.join(request.config.rootdir, sz_data_file), + "cov_ksz_file": os.path.join(request.config.rootdir, cov_ksz_file), + "twohalo_term": os.path.join(request.config.rootdir, twohalo_term), + "beam_file": os.path.join(request.config.rootdir, beam_file), + } + }, + } + + thta_arc, ksz_fft, tsz_fft, ksz_hank, tsz_hank = np.loadtxt( + os.path.join(request.config.rootdir, sz_data_file), + usecols=(0, 1, 2, 3, 4), + unpack=True, + ) + + model = get_model(info) + loglikes, derived = model.loglikes() + + lhood = model.likelihood["KSZLikelihood"] + + ksz_theory = lhood._get_theory(**info["params"]) + ksz_theory *= 3282.8 * 60.0**2 + print("ksz_theory", ksz_theory) + print("ksz_hank", ksz_hank) + assert np.allclose(ksz_hank, ksz_theory, atol=1.0e-3, rtol=1.0e-3) + + +def test_tsz(request): + info = { + "params": { + "gnfw_P0": 9.10766709141, + "gnfw_bt_tsz": 4.76625102519, + "gnfw_A2h_tsz": 1.0, + }, + "likelihood": { + "TSZLikelihood": { + "external": TSZLikelihood, + "sz_data_file": os.path.join(request.config.rootdir, sz_data_file), + "cov_tsz_file": os.path.join(request.config.rootdir, cov_tsz_file), + "twohalo_term": os.path.join(request.config.rootdir, twohalo_term), + "beam_file": os.path.join(request.config.rootdir, beam_file), + "beam_response": os.path.join(request.config.rootdir, beam_response), + } + }, + } + + thta_arc, ksz_fft, tsz_fft, ksz_hank, tsz_hank = np.loadtxt( + os.path.join(request.config.rootdir, sz_data_file), + usecols=(0, 1, 2, 3, 4), + unpack=True, + ) + + model = get_model(info) + loglikes, derived = model.loglikes() + + lhood = model.likelihood["TSZLikelihood"] + + tsz_theory = lhood._get_theory(**info["params"]) + tsz_theory *= 3282.8 * 60.0**2 + print("tsz_theory", tsz_theory) + print("tsz_hank", tsz_hank) + assert np.allclose(tsz_hank, tsz_theory, atol=1.0e-3, rtol=1.0e-3) diff --git a/tests/test_szlike.yaml b/tests/test_szlike.yaml new file mode 100644 index 00000000..123a0e88 --- /dev/null +++ b/tests/test_szlike.yaml @@ -0,0 +1,43 @@ +sz_data_file: gnfw_test_projections.txt +beam_file: beam_f150_daynight.txt +transform_type: Hankel #FFT or Hankel +#optional beam response file +beam_response: act_planck_s08_s18_cmb_f150_daynight_response_tsz.txt #False, or give file +#optional 2-halo file +twohalo_term: twohalo_cmass_average.txt + +redshift: 0.55 +mass_halo_mean_Msol: 10**13. +frequency_GHz: 150. +input_model: "gnfw" #or "obb" + +params: + Omega_m: + value: 0.25 + Omega_b: + value: 0.044 + hh: + value: 0.7 + Omega_L: + value: 0.75 + rhoc_0: + value: 2.77525e2 + C_OVER_HUBBLE: + value: 2997.9 + XH: + value: 0.76 + v_rms: + value: 1.06e-3 + + gnfw_rho0: + value: 3.1864410705359125 + gnfw_bt_ksz: + value: 3.45493977635 + gnfw_A2h_ksz: + value: 1.0 + gnfw_P0: + value: 9.10766709141 + gnfw_bt_tsz: + value: 4.76625102519 + gnfw_A2h_tsz: + value: 1.0