commit 131567e9b16da54bbfa07a934c814fa5de33732c Author: Tibeuleu Date: Fri Feb 13 14:45:04 2026 +0100 Initial commit diff --git a/.gitignore b/.gitignore new file mode 100755 index 0000000..24f2996 --- /dev/null +++ b/.gitignore @@ -0,0 +1,21 @@ +# Byte-compiled / optimized / DLL files +__pycache__/ + +# Jupyter Notebook +.ipynb_checkpoints + +# IPython +profile_default/ +ipython_config.py + +# Ignore any .fits/.npy/.par data file +*.fits +*.npy +*.par +*.zip +*.tar +*.tar.xz + +# Remove local data and plots folders +data/ +plots/ \ No newline at end of file diff --git a/README.md b/README.md new file mode 100755 index 0000000..37e3998 --- /dev/null +++ b/README.md @@ -0,0 +1,2 @@ +# Project 3 : SZ signal recovery from stacking and characterization of signals +users : t.barnouin / k.zhang diff --git a/report/P3_SZ_Signal_Using_Stacking.pdf b/report/P3_SZ_Signal_Using_Stacking.pdf new file mode 100755 index 0000000..8455b15 Binary files /dev/null and b/report/P3_SZ_Signal_Using_Stacking.pdf differ diff --git a/report/README.md b/report/README.md new file mode 100755 index 0000000..b67228a --- /dev/null +++ b/report/README.md @@ -0,0 +1 @@ +# t.barnouin : directory where we will put our report diff --git a/report/Spectral Energy Distribution of SZ Catalogue.pdf b/report/Spectral Energy Distribution of SZ Catalogue.pdf new file mode 100755 index 0000000..85df6bc Binary files /dev/null and b/report/Spectral Energy Distribution of SZ Catalogue.pdf differ diff --git a/src/README.md b/src/README.md new file mode 100755 index 0000000..44b260d --- /dev/null +++ b/src/README.md @@ -0,0 +1 @@ +# t.barnouin : directory where all our code will be diff --git a/src/SED.py b/src/SED.py new file mode 100755 index 0000000..2535308 --- /dev/null +++ b/src/SED.py @@ -0,0 +1,208 @@ +#!/usr/bin/ipython3 +#-*- coding:utf-8 -*- + +import sys, time +import concurrent.futures +import numpy as np +import stack_maps as stack +from astropy.io import fits +import matplotlib.pyplot as plt +from scipy.optimize import curve_fit + +def aperture(stack, px_size, center, angle): + """ + Create a boolean mask for the region to work on the stack + """ + yy, xx = np.indices(stack.shape) + radius = (angle/60.)/px_size + mask = np.sqrt(np.abs(xx-center[0])**2 + np.abs(yy-center[1])**2) <= radius + + return mask + + +def stack_flux(freq, stacked_data, px_size, center, angle): + """ + Compute integrated flux in Mjy/sr of stacked image given an aperture angle + """ + conversion = {100: 244.1, 143: 371.74, 217: 483.690, 353: 287.450, 545: 1., 857: 1.} # Conversion coefficient in Mjy/sr/K + mask = aperture(stacked_data, px_size, center, angle) + region = stacked_data[mask] + background = np.mean(stacked_data[~mask]) + flux = np.sum(region - background)*conversion[freq] + + return flux + + +def bootstrap_flux(data_sample, wcs_sample, freq, angle, n_it, savename=None): + """ + Implement the 'bootstrap' technic to compute flux and error given a sample. + """ + n_obj = data_sample.shape[0] + flux = np.zeros(n_it) + for i in range(n_it): + ind_it = np.random.randint(0, n_obj, n_obj) + data_sample_it = data_sample[ind_it] + wcs_sample_it = wcs_sample[ind_it] + + stack_it, px_size, center = stack.stack_sample(data_sample_it, wcs_sample_it) + flux[i] = stack_flux(freq, stack_it, px_size, center, angle) + + if savename != None: + np.save("../data/flux_{0:s}_{1:d}.npy".format(savename,freq),np.array([flux.mean(), flux.std()])) + + return flux.mean(), flux.std() + + +def parallel_bootstrap(data_sample, wcs_sample, freq, angle, n_it, savename=None): + """ + Implement the 'bootstrap' technic with multithreading for multiple frequencies + """ + nfreq = len(freq) + exe = concurrent.futures.ThreadPoolExecutor() + future_flux = [] + flux_sample = np.zeros((nfreq, 2)) + for i in range(nfreq): + future_flux.append(exe.submit(bootstrap_flux, data_sample[i], wcs_sample[i], freq[i], angle, n_it)) + for i in range(nfreq): + flux_sample[i] = future_flux[i].result() + + if savename != None: + np.save("../data/flux_{0:s}.npy".format(savename),flux_sample) + + return flux_sample + + +def s_SZ(nu, y_SZ): + """ + SED law for thermal SZ effect, in Mjy/sr + """ + h = 6.62607015e-34 #SI units + k_b = 1.380649e-23 + T_cmb = 2.72548 + x = h/(k_b*T_cmb)*nu #reduced frequency + + g = x*(np.exp(x)+1.)/(np.exp(x)-1.)-4. + return y_SZ*g*x**4*np.exp(x)/(np.exp(x)-1.)**2 + + +def s_d(nu, A_d, beta_d, T_d): + """ + SED law for thermal dust, in Mjy/sr (appear in IR) + """ + h = 6.62607015e-34 #SI units + k_b = 1.380649e-23 + nu_0 = 545e9 + gamma = h/(k_b*T_d) + return A_d*(nu/nu_0)**(beta_d+1.)*(np.exp(gamma*nu_0)-1.)/(np.exp(gamma*nu)-1.) + + +def flux_fit(nu, y_SZ, A_d, beta_d, T_d): + """ + SED law containing thermal SZ effect and IR contamination by thermal dust + """ + return s_SZ(nu, y_SZ)+s_d(nu, A_d, beta_d, T_d) + + +def plot_SED(freq, flux, samplename, save=False): + """ + Nicely plot SED and contributions from SZ and dust for a given sample + """ + freq_lin = np.linspace(freq.min(), freq.max(), 300) + popt, _ = curve_fit(flux_fit, xdata=freq*1e9, ydata=flux[:,0], p0=[1., 10., 2.3, 16.]) + chisq = np.sum(((flux[:,0]-flux_fit(freq,*popt))/flux[:,1])**2) + + plt.figure(figsize=(10,5)) + plt.errorbar(freq,flux[:,0], yerr=flux[:,1], fmt='ok', label=samplename) + plt.plot(freq_lin,s_SZ(freq_lin*1e9,popt[0]), 'b--', label=r"$S_{SZ}$ model fit with $y_{SZ}$ = "+"{0:.5f}".format(popt[0])) + plt.plot(freq_lin,s_d(freq_lin*1e9,*popt[1:]), 'r--', label=r"$S_{d}$ model fit with $A_{d}$ = "+"{0:.3e}".format(popt[1])+r"; $\beta_{d}$ = "+"{0:.5f}".format(popt[2])+r"; $T_{d}$ = "+"{0:.5f}\n".format(popt[3])) + plt.plot(freq_lin,flux_fit(freq_lin*1e9, *popt), 'g--', label=r"$S_{SZ}+S_{d}$ model fit with $\chi^2$/ndf = "+"{0:.5f}".format(chisq/(len(freq)-len(popt)))) + plt.plot(freq,[0.]*len(freq), 'k--') + + plt.title("SED of the SZ effect for '{0:s}' sample.".format(samplename)) + plt.xlabel("Frequency [GHz]") + plt.ylabel("Flux density [Mjy/sr]") + plt.legend() + if save: + plt.savefig("../plots/SED_{0:s}.png".format(samplename), bbox_inches='tight') + plt.show() + + return 0 + + +def compute_SED(sample, savename, ROI=[0.5, 128], angle=None, n_it=500, plot=False): + """ + Automatically compute and save SED + """ + freq = np.array([100, 143, 217, 353, 545, 857]) + + if angle == None: + rimo = fits.getdata("/home/abeelen/Planck/products_2015/HFI_RIMO_R2.00.fits", 2) + true_FWHM = rimo[:]['FWHM'] + angle = 1.*true_FWHM[freq == 143] + + print("Getting sample data for {0:s}...".format(savename)) + t0 = time.time() + wcs_sample, data_sample = stack.get_data_sample(freq, sample, ROI=ROI, savename=savename) + t1 = time.time() + print("Elapsed time : {0:f}sec".format(t1-t0)) + + stack_freq = np.zeros((len(freq), data_sample.shape[2], data_sample.shape[3])) + for i in range(len(freq)): + stack_freq[i] = stack.stack_sample(data_sample[i], wcs_sample[i])[0] + np.save("../data/stack_{0:s}.npy".format(savename),stack_freq) + + print("Computing flux values...") + t2 = time.time() + flux_sample = parallel_bootstrap(data_sample, wcs_sample, freq, angle, n_it, savename=savename) + t3 = time.time() + print("Elapsed time : {0:f}sec".format(t3-t2)) + + if plot: + plot_SED(freq, flux_sample, savename, save=True) + + return 0 + + +def main(): + # General information on the Planck data channels + rimo = fits.getdata("/home/abeelen/Planck/products_2015/HFI_RIMO_R2.00.fits", 2) + true_FWHM = rimo[:]['FWHM'] + + frequencies = np.array([100, 143, 217, 353, 545, 857]) + + # Get objects sample from PSZ catalog + with fits.open("../data/PSZ2v1.fits") as hdu: + PSZ_catalog = hdu[1].data + + clusters = PSZ_catalog[PSZ_catalog['redshift'] > 0] + candidates = PSZ_catalog[(PSZ_catalog['redshift'] < 0) & (PSZ_catalog['validation'] < 0)] +# redshift_sup0_4 = PSZ_catalog[PSZ_catalog['redshift'] > 0.4] +# redshift_inf0_4 = PSZ_catalog[np.abs(PSZ_catalog['redshift']) < 0.4] +# QN_inf0_4 = PSZ_catalog[PSZ_catalog['q_neural'] < 0.4] + QN_sup0_4 = PSZ_catalog[PSZ_catalog['q_neural'] > 0.4] + + # Compute SED and plot + angle = 1.*true_FWHM[frequencies == 143] + + #compute_SED(clusters, "clusters", ROI=[0.5, 128], angle=angle, n_it=500, plot=True) + #compute_SED(candidates, "candidates", ROI=[0.5, 128], angle=angle, n_it=500, plot=True) + #compute_SED(QN_sup0_4, "QN_sup0_4", ROI=[0.5, 128], angle=angle, n_it=500, plot=True) + flux_clusters = np.load("../data/flux_clusters.npy") + flux_candidates = np.load("../data/flux_candidates.npy") + flux_QN_sup0_4 = np.load("../data/flux_QN_sup0_4.npy") + + plot_SED(frequencies, flux_clusters, "Known clusters", save=True) + plot_SED(frequencies, flux_candidates, "Candidate clusters", save=True) + plot_SED(frequencies, flux_QN_sup0_4, "High-reliability clusters", save=True) + + stack_clusters = np.load("../data/stack_clusters.npy") + stack_candidates = np.load("../data/stack_candidates.npy") + stack_QN_sup0_4 = np.load("../data/stack_QN_sup0_4.npy") + + stack.plot_stacks(np.array([stack_clusters, stack_candidates, stack_QN_sup0_4]), ["Known clusters", "Candidate clusters", "High-reliability clusters"], save=True) + + return 0 + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/src/stack_maps.py b/src/stack_maps.py new file mode 100755 index 0000000..d5dc336 --- /dev/null +++ b/src/stack_maps.py @@ -0,0 +1,246 @@ +#!/usr/bin/ipython3 +#-*- coding:utf-8 -*- + +import sys +import matplotlib.pyplot as plt +from astropy.io import fits +from astropy import wcs +import numpy as np +import healpy as hp + + +def create_WCS(sky_pos, ROI=[0.5, 128], radec=False): + """ + Create a WCS object centered at given sky position. + """ + w = wcs.WCS(naxis=2) + center = np.array([int(ROI[1]/2.), int(ROI[1]/2.)]) + w.wcs.crpix = center + w.wcs.crval = sky_pos + w.wcs.cdelt = [-ROI[0]/60., ROI[0]/60.] + + if radec: + w.wcs.ctype = ['RA--TAN','DEC--TAN'] + else: + w.wcs.ctype = ['GLON-TAN','GLAT-TAN'] + + return w + + +def create_ROI(hpmap, sky_pos, ROI=[0.5, 128], radec=False): + """ + Create a projected map given a Healpix map and a region of interest. + """ + nside = hp.get_nside(hpmap) + w = create_WCS(sky_pos, ROI, radec) + + data = np.zeros((ROI[1], ROI[1])) + yy, xx = np.indices(data.shape) + theta, phi = w.wcs_pix2world(xx, yy, 0) + ipx = hp.ang2pix(nside, theta, phi, lonlat=True) + data = hpmap[ipx] + + return w, data + + +def plot_ROI(w, data, name=None, freq=None, pos=None, saveplot=False): + """ + Display a nice plot of the given data. + """ + fig = plt.figure(figsize=(10,5)) + ax = plt.subplot(projection=w, label='overlays', figure=fig) + + ax.imshow(data, origin='lower') + + ax.coords.grid(True, color='white', ls='solid') + ax.coords[0].set_axislabel('Galactic Longitude') + ax.coords[1].set_axislabel('Galactic Latitude') + + overlay = ax.get_coords_overlay('fk5') + overlay.grid(color='white', ls='dotted') + overlay[0].set_axislabel('Right Ascension (J2000)') + overlay[1].set_axislabel('Declination (J2000)') + + if name != None: + ax.set_title("Object {0:s} at {1:d}GHz".format(name, freq), y=1.0, pad=+40) + if saveplot: + plt.savefig("../plots/{0:s}_G{1:.2f}±{2:.2f}_{2:d}GHz.png".format(name, *pos, freq), bbox_inches='tight') + plt.show() + + return 0 + + +def save_file(map_name,sky_pos,header,data): + """ + Save given projected map to fits file with appropriate header. + """ + hdu = fits.PrimaryHDU(data,header) + filename = '../data/' + str(map_name) + "_G{0:.2f}±{1:.2f}".format(*sky_pos) + '.fits' + + hdu.writeto(filename, overwrite=True) + + return 0 + + +def get_data_sample(freq, sample, ROI=[0.5, 128], savename=None): + """ + Construct on array (nobj, npix, npix) at a given frequency for a given sample and ROI parameters. Allow local saving. + """ + if type(freq) == int: + freq = np.array([freq]) + nfreq = len(freq) + nobj = len(sample) + + # Get Planck data + dirf = "../data/Planck2015/" + fits_basename = "HFI_SkyMap_{0:d}_2048_R2.00_full.fits" + + data_sample = np.zeros((nfreq, nobj, ROI[1], ROI[1])) + wcs_sample = np.zeros((nfreq, nobj), dtype=object) + for i in range(nfreq): + hpmap = hp.read_map(dirf+fits_basename.format(freq[i]), verbose=False) + for j in range(nobj): + pos_obj = [sample['GLON'][j], sample['GLAT'][j]] + wcs_obj, data_obj = create_ROI(hpmap, pos_obj, ROI, radec=False) + data_sample[i,j] = data_obj + wcs_sample[i,j] = wcs_obj + + if nfreq == 1: + wcs_sample, data_sample = wcs_sample[0], data_sample[0] + + if savename != None: + np.save("../data/wcs_{0:s}.npy".format(savename), wcs_sample) + np.save("../data/data_{0:s}.npy".format(savename), data_sample) + + return wcs_sample, data_sample + + +def save_sample(sample, sample_name, freq, ROI=[0.5, 128]): #sample from PSZ2v1.fits, freq from planck full sky map + # Get Planck data + dirf = "../data/Planck2015/" + fits_basename = "HFI_SkyMap_{0:d}_2048_R2.00_full.fits" + hpmap = hp.read_map(dirf+fits_basename.format(freq), verbose=False) + + hdul = fits.HDUList([]) + primary_hdu = fits.PrimaryHDU() + hdul.append(primary_hdu) + hdr = hdul[0].header + hdr['sample'] = (sample_name, 'Same of the sample') + hdr['freq'] = (freq, 'Detection frequency in GHz') + hdr['nobj'] = len(sample) + hdr['ROI_size'] = (ROI[1], 'Size in pixels of the ROI (square)') + hdr['px_size'] = (ROI[0], 'Size of a pixel in arcmin') + + for i in range(len(sample)): + sample_pos = [sample['GLON'][i], sample['GLAT'][i]] + sample_w, sample_data = create_ROI(hpmap, sample_pos, ROI, radec=False) + + sample_header = sample_w.to_header() + sample_header['targname'] = sample['NAME'][i] + + hdu_sample = fits.ImageHDU(data=sample_data,header=sample_header) + hdul.append(hdu_sample) + + filename = '../data/' + str(sample_name) + "_{0:d}GHz".format(freq) + '.fits' + hdul.writeto(filename, overwrite=True) + + return 0 + + +def data_sample_from_fits(sample): + """ + Construct an array (ncluster, npix, npix) of 'sample'. + Input 'sample' must be a valid path to a fits file. + """ + with fits.open(sample) as hdul: + header = hdul[0].header + + nobj = header['nobj'] + ROI_size = header['ROI_size'] + data_sample = np.zeros((nobj, ROI_size, ROI_size)) + wcs_sample = np.zeros(nobj,dtype=object) + + for i in range(nobj): + data_sample[i] = hdul[i+1].data + wcs_sample[i] = wcs.WCS(hdul[i+1].header) + + return wcs_sample, data_sample + + +def stack_sample(data_sample, wcs_sample=None): + """ + Compute the stack of a given sample. + Input 'sample' can either be an array (ncluster, npixel, npixel) or a path to fits file containing the sample. + """ + if str(type(data_sample)) != "": + wcs_sample, data_sample = data_sample_from_fits(sample) + + data_stack = data_sample.mean(axis=0) + px_size = np.array([np.abs(wcs_sample[i].wcs.cdelt) for i in range(wcs_sample.shape[0])]).mean() + center = np.array([wcs_sample[i].wcs.crpix for i in range(wcs_sample.shape[0])]).mean(axis=0) + + return data_stack, px_size, center + + +def plot_stacks(stack_list, stack_names, save=False): + """ + Plot together stacks at each frequencies for different samples. + """ + freq = np.array([100, 143, 217, 353, 545, 857]) + nfreq, nsample = len(freq), len(stack_list) + + grid = plt.GridSpec(nsample, nfreq, hspace=0, wspace=0) + fig = plt.figure(figsize=(16,9)) + for i in range(nsample): + for j in range(nfreq): + ax = fig.add_subplot(grid[i,j], xticks=[], yticks=[]) + ax.imshow(stack_list[i,j], origin='lower') + if i == nsample-1: + ax.set_xlabel("{0:d}GHz".format(freq[j])) + if j == 0: + ax.set_ylabel(stack_names[i]) + if save: + filename = "../plots/Stacks" + for name in stack_names: + filename += "_"+name + plt.savefig(filename+".png", bbox_inches='tight') + plt.show() + return 0 + + +def main(): + # Get Planck data + dirf = "../data/Planck2015/" + ext = ".fits" + fits_basename = "HFI_SkyMap_{0:d}_2048_R2.00_{1:s}" + mission = "full" #Can be 'full', 'halfmission-1', 'halfmission-2' + frequencies = np.array([100, 143, 217, 353, 545, 857]) + + hfi_maps = np.array([fits_basename.format(freq,mission) for freq in frequencies]) + + ## Get known clusters and candidates + with fits.open("../data/PSZ2v1.fits") as hdu: + PSZ_catalog = hdu[1].data + + clusters = PSZ_catalog[PSZ_catalog['redshift']>0] + candidates = PSZ_catalog[(PSZ_catalog['redshift']<0) & (PSZ_catalog['validation']<0)] + +# wcs_clusters, data_clusters = get_data_sample(frequencies, clusters, ROI=[0.5, 128], savename="clusters") +# wcs_candidates, data_candidates = get_data_sample(frequencies, candidates, ROI=[0.5, 128], savename="candidates") + wcs_clusters, data_clusters = np.load("../data/wcs_clusters.npy"), np.load("../data/data_clusters.npy") + wcs_candidates, data_candidates = np.load("../data/wcs_candidates.npy"), np.load("../data/data_candidates.npy") + + stacks_clusters = np.zeros((len(frequencies), data_clusters.shape[2], data_clusters.shape[3])) + stacks_candidates = np.zeros((len(frequencies), data_candidates.shape[2], data_candidates.shape[3])) + + for i in range(len(frequencies)): + stacks_clusters[i] = stack_sample(data_clusters[i], wcs_clusters[i])[0] + stacks_candidates[i] = stack_sample(data_candidates[i], wcs_candidates[i])[0] + + plot_stacks(np.array([stacks_clusters, stacks_candidates]), ['clusters', 'candidates'], save=True) + + return 0 + + +if __name__ == "__main__": + sys.exit(main())