diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 3ea10d3e5..ecd65ce66 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -273,6 +273,7 @@ jobs: ceci --dry-run examples/buzzard/pipeline.yml ceci --dry-run examples/cosmodc2/pipeline.yml ceci --dry-run examples/skysim/pipeline.yml + ceci --dry-run examples/desy3/pipeline_Bmodes.yml Test_Auto_Installer: diff --git a/examples/desy1/pipeline_Bmodes.yml b/examples/desy1/pipeline_Bmodes.yml new file mode 100644 index 000000000..0f6890305 --- /dev/null +++ b/examples/desy1/pipeline_Bmodes.yml @@ -0,0 +1,116 @@ + +# Stages to run +stages: + + - name: TXSourceSelectorMetacal # select and split objects into source bins + - name: TXShearCalibration # Calibrate and split the source sample tomographically + #- name: TXTruthLensCatalogSplitter + #- name: TXStarCatalogSplitter # Split the star catalog into separate bins (psf/non-psf) + #- name: TXPhotozPlot + - name: TXSourceMaps + #nprocess: 6 + #threads_per_process: 1 + #nodes: 1 + #- name: TXLensMaps + #- name: TXAuxiliarySourceMaps # make PSF and flag maps + #- name: TXAuxiliaryLensMaps # make depth and bright object maps + #- name: TXSourceNoiseMaps # noise maps for sources using random rotations + # stage below is not working because TXAuxiliaryLensMaps is broken + #- name: TXLensNoiseMaps # noise maps for lenses using density splits + #- name: TXDensityMaps # turn mask and ngal maps into overdensity maps + #- name: TXMapPlots # make pictures of all the maps + #- name: TXTracerMetadata # collate metadata + #- name: TXJackknifeCenters # Split the area into jackknife regions + #- name: TXRandomCat + #- name: TXTwoPoint # Compute real-space 2-point correlations + # threads_per_process: 128 + #- name: TXBlinding # Blind the data following Muir et al + #threads_per_process: 2 + #- name: TXTwoPointTheoryReal # compute theory using CCL to save in sacc file and plot later + # - name: TXTwoPointPlots # Make plots of 2pt correlations + #- name: TXSourceHistogramPlots + #- name: TXSourceDiagnosticPlots # Make a suite of diagnostic plots + #nprocess: 8 + #nodes: 2 + #- name: TXLensDiagnosticPlots # Make a suite of diagnostic plots + #- name: TXGammaTFieldCenters # Compute and plot gamma_t around center points # ERROR: EXPOSURE FILE + ##threads_per_process: 2 + #- name: TXGammaTStars # Compute and plot gamma_t around bright stars + #threads_per_process: 2 + #- name: TXGammaTRandoms # Compute and plot gamma_t around randoms + #threads_per_process: 2 + #- name: TXRoweStatistics # Compute and plot Rowe statistics + # threads_per_process: 128 + #- name: TXFocalPlanePlot # mean PSF ellipticity and mean residual on the focal plane + #- name: TXGalaxyStarDensity # Compute and plot the star-galaxy density cross-correlation + #- name: TXGalaxyStarShear # Compute and plot the star-galaxy shear cross-correlation + #- name: TXPSFDiagnostics # Compute and plots other PSF diagnostics + #- name: TXBrighterFatterPlot # Make plots tracking the brighter-fatter effect + #- name: TXConvergenceMaps # Make convergence kappa maps from g1, g2 maps #ERROR: WLMASSMAP + #- name: TXConvergenceMapPlots # Plot the convergence map + #- name: TXMapCorrelations # plot the correlations between systematics and data # ERROR: NEEDS CONVERGENCE MAP + #- name: TXApertureMass # Compute aperture-mass statistics + ##threads_per_process: 2 + # Disabled since not yet synchronised with new Treecorr MPI + # - name: TXSelfCalibrationIA # Self-calibration intrinsic alignments of galaxies + # + # Disabling these as they takes too long for a quick test. + # The latter three need NaMaster + + ##- name: TXRealGaussianCovariance # Compute covariance of real-space correlations + - name: TXTwoPointFourier # Compute power spectra C_ell + ##- name: TXFourierNamasterCovariance # Compute the C_ell covariance + ##- name: TXFourierTJPCovariance # Compute the C_ell covariance with TJPCov + +#=========== + +launcher: + name: mini + interval: 1.0 + +site: + name: local + max_threads: 128 + +modules: txpipe + +# where to find any modules that are not in this repo, +# and any other code we need. +python_paths: + - submodules/WLMassMap/python/desc/ + - submodules/TJPCov + +# where to put outputs +output_dir: data/desy1/outputs + +# configuration settings +config: examples/desy1/config.yml + +# These are overall inputs to the whole pipeline, not generated within it +inputs: + shear_catalog: /global/cfs/cdirs/lsst/groups/WL/projects/txpipe-sys-tests/des-y1/shear_catalog_desy1_masked.h5 + # subsampled catalog for debugging: + #shear_catalog: /global/cscratch1/sd/jjeffers/TXPipe/data/desy1/inputs/small_DESY1_shear_catalog.h5 + photometry_catalog: /global/cfs/cdirs/lsst/groups/WL/projects/txpipe-sys-tests/des-y1/photometry_catalog_desy1_RM.h5 + lens_tomography_catalog: /global/cfs/cdirs/lsst/groups/WL/projects/txpipe-sys-tests/des-y1/lens_tomography_catalog_desy1_RM_091423.h5 + shear_photoz_stack: /global/cfs/cdirs/lsst/groups/WL/projects/txpipe-sys-tests/des-y1/shear_photoz_stack.hdf5 + lens_photoz_stack: /global/cfs/cdirs/lsst/groups/WL/projects/txpipe-sys-tests/des-y1/lens_photoz_stack.hdf5 + calibration_table: data/example/inputs/sample_cosmodc2_w10year_errors.dat + exposures: data/example/inputs/exposures.hdf5 + star_catalog: /global/cfs/cdirs/lsst/groups/WL/projects/txpipe-sys-tests/des-y1/DES_psf_catalog.hdf5 + fiducial_cosmology: data/fiducial_cosmology.yml + # For the self-calibration extension we are not using Random_cat_source for now + # So it is set to Null, so the yaml intepreter returns a None value to python. + random_cats_source: Null + mask: /global/cfs/cdirs/lsst/groups/WL/projects/txpipe-sys-tests/des-y1/mask_desy1.h5 + #random_cats: /global/cfs/cdirs/lsst/groups/WL/projects/txpipe-sys-tests/des-y1/randoms_desy1_RM.hdf5 + #binned_random_catalog: /global/cfs/cdirs/lsst/groups/WL/projects/txpipe-sys-tests/des-y1/binned_randoms_desy1_RM.hdf5 +# if supported by the launcher, restart the pipeline where it left off +# if interrupted +resume: True + +# where to put output logs for individual stages +log_dir: data/desy1/logs + +# where to put an overall parsl pipeline log +pipeline_log: data/desy1/log.txt diff --git a/examples/desy3/config.yml b/examples/desy3/config.yml index 937905ab0..c188a1194 100644 --- a/examples/desy3/config.yml +++ b/examples/desy3/config.yml @@ -60,7 +60,7 @@ TXSourceSelectorMetacal: delta_gamma: 0.02 source_zbin_edges: [0.0, 0.2, 0.4, 0.6, 0.8] shear_prefix: mcal_ - resp_mean_diag : True + use_diagonal_response : True TXSourceSelectorLensfit: bands: ri #used for selection diff --git a/examples/desy3/config_Bmodes.yml b/examples/desy3/config_Bmodes.yml new file mode 100644 index 000000000..329dd4312 --- /dev/null +++ b/examples/desy3/config_Bmodes.yml @@ -0,0 +1,219 @@ +# Values in this section are accessible to all the different stages. +# They can be overridden by individual stages though. +global: + # This is read by many stages that read complete + # catalog data, and tells them how many rows to read + # at once + chunk_rows: 100000 + # These mapping options are also read by a range of stages + pixelization: healpix + nside: 512 + sparse: True # Generate sparse maps - faster if using small areas + +PZPrepareEstimatorLens: + name: PZPrepareEstimatorLens + classname: Inform_BPZ_lite + aliases: + input: spectroscopic_catalog + model: lens_photoz_model + zmin: 0.0 + zmax: 3.0 + nzbins: 301 + columns_file: ./data/bpz_ugrizy.columns + data_path: ./data/example/rail-bpz-inputs + spectra_file: CWWSB4.list + prior_band: i + ref_band: i + # Not sure about this + prior_file: hdfn_gen + p_min: 0.005 + gauss_kernel: 0.0 + mag_err_min: 0.005 + inform_options: {'save_train': False, 'load_model': False, 'modelfile': 'BPZpriormodel.out'} + madau_reddening: no + bands: riz + zp_errors: [0.01, 0.01, 0.01, 0.01, 0.01, 0.01] + hdf5_groupname: photometry + + +PZEstimatorLens: + name: PZEstimatorLens + classname: BPZ_lite + aliases: + model: lens_photoz_model + input: photometry_catalog + output: lens_photoz_pdfs + zmin: 0.0 + zmax: 3.0 + dz: 0.01 + nzbins: 301 + band_names: [mag_u, mag_g, mag_r, mag_i, mag_z, mag_y] + bands: [mag_u, mag_g, mag_r, mag_i, mag_z, mag_y] + err_bands: [mag_err_u, mag_err_g, mag_err_r, mag_err_i, mag_err_z, mag_err_y] + hdf5_groupname: photometry + nondetect_val: .inf + columns_file: ./data/bpz_ugrizy.columns + data_path: ./data/example/rail-bpz-inputs + spectra_file: CWWSB4.list + prior_band: mag_i + ref_band: mag_i + prior_file: hdfn_gen + p_min: 0.005 + gauss_kernel: 0.0 + zp_errors: [0.01, 0.01, 0.01, 0.01, 0.01, 0.01] + mag_err_min: 0.005 + madau_reddening: false + mag_limits: + mag_u: 27.79 + mag_g: 29.04 + mag_r: 29.06 + mag_i: 28.62 + mag_z: 27.98 + mag_y: 27.05 + +TXSourceSelectorMetacal: + input_pz: True + true_z: False + bands: riz #used for selection + T_cut: 0.5 + s2n_cut: 10.0 + max_rows: 10000000000 + delta_gamma: 0.02 + source_zbin_edges: [0.0, 0.2, 0.4, 0.6, 0.8] + shear_prefix: mcal_ + use_diagonal_response : True + +TXMeanLensSelector: + # Mag cuts + lens_zbin_edges: [0.0,0.2,0.4] + cperp_cut: 0.2 + r_cpar_cut: 13.5 + r_lo_cut: 16.0 + r_hi_cut: 21.6 + i_lo_cut: 17.5 + i_hi_cut: 21.9 + r_i_cut: 2.0 + +TXModeLensSelector: + # Mag cuts + lens_zbin_edges: [0.0,0.2,0.4] + cperp_cut: 0.2 + r_cpar_cut: 13.5 + r_lo_cut: 16.0 + r_hi_cut: 21.6 + i_lo_cut: 17.5 + i_hi_cut: 21.9 + r_i_cut: 2.0 + +TXRandomCat: + density: 10 # gals per sq arcmin + +TXTwoPoint: + bin_slop: 0.1 + delta_gamma: 0.02 + do_pos_pos: True + do_shear_shear: True + do_shear_pos: True + flip_g2: True # use true when using metacal shears + min_sep: 2.5 + max_sep: 60.0 + nbins: 1000 + verbose: 0 + subtract_mean_shear: True + var_method: "jackknife" + +TXSourceDiagnosticPlots: + psf_prefix: mcal_psf_ + shear_prefix: mcal_ + nbins: 20 + g_min: -0.01 + g_max: 0.01 + psfT_min: 0.2 + psfT_max: 0.28 + T_min: 0.01 + T_max: 2.1 + s2n_min: 1.25 + s2n_max: 300 + +TXSourceMaps: + sparse: True # Generate sparse maps - faster if using small areas + +TXLensMaps: + sparse: True # Generate sparse maps - faster if using small areas + nside: 128 + +TXDensityMaps: + nside: 128 + +TXAuxiliaryLensMaps: + flag_exponent_max: 8 + bright_obj_threshold: 22.0 # The magnitude threshold for a object to be counted as bright + depth_band : i + snr_threshold: 10.0 # The S/N value to generate maps for (e.g. 5 for 5-sigma depth) + snr_delta: 1.0 # The range threshold +/- delta is used for finding objects at the boundary + +TXJackknifeCenters: + npatch: 5 + +TXSimpleMask: + depth_cut : 23.5 + bright_object_max: 10.0 + +PZRailSummarizeLens: + leafsize: 20 + zmin: 0.0 + zmax: 3.0 + nzbins: 50 + name: PZRailSummarizeLens + catalog_group: "photometry" + tomography_name: "lens" + aliases: + tomography_catalog: lens_tomography_catalog + photometry_catalog: photometry_catalog + model: lens_direct_calibration_model + photoz_stack: lens_photoz_stack + model: None + +FlowCreator: + n_samples: 1000000 + seed: 5763248 + aliases: + # This input was generated using get_example_flow in pzflow, + # not something specific. + output: ideal_specz_catalog + model: flow + +TXParqetToHDF: + hdf_group: photometry + aliases: + input: specz_catalog_pq + output: spectroscopic_catalog + + +Inform_NZDirSource: + name: Inform_NZDirSource + usecols: [r, i, z] + hdf5_groupname: photometry + aliases: + input: spectroscopic_catalog + model: source_direct_calibration_model + +Inform_NZDirLens: + name: Inform_NZDirLens + usecols: [u, g, r, i, z, "y"] + hdf5_groupname: photometry + aliases: + input: spectroscopic_catalog + model: lens_direct_calibration_model + +TXPureB: + method : 'namaster' # 'namaster' or 'hybrideb' + Nell : 20 # Number of ell bins + lmin : 100 # ell min + lmax : 2000 # ell max + lspacing : 'log' # ell spacing for binning + bin_file : # Load predfined bin edges instead + theta_min : 2.5 # (only for hybrideb) theta min in arcmin + theta_max : 250 # (only for hybrideb) theta max in arcmin + Ntheta : 1000 # (only for hybrideb) Number of theta bins + Nsims : 3 # (only for hybrideb) Number of Gaussian draws to use to estimate errors diff --git a/examples/desy3/pipeline_Bmodes.yml b/examples/desy3/pipeline_Bmodes.yml new file mode 100644 index 000000000..e08a54956 --- /dev/null +++ b/examples/desy3/pipeline_Bmodes.yml @@ -0,0 +1,92 @@ + +# Stages to run +stages: + - name: TXPSFDiagnostics # Compute and plots other PSF diagnostics + - name: TXAuxiliaryLensMaps # make depth and bright object maps + - name: TXSimpleMask # combine maps to make a simple mask + - name: TXSourceSelectorMetacal # select and split objects into source bins + - name: TXSourceNoiseMaps # Compute shear noise using rotations + - name: TXShearCalibration # Calibrate and split the source sample tomographically + - name: TXSourceMaps # make source g1 and g2 maps + - name: TXConvergenceMaps # Make convergence kappa maps from g1, g2 maps + - name: FlowCreator # Simulate a spectroscopic population + - name: GridSelection # Simulate a spectroscopic sample + - name: TXParqetToHDF # Convert the spec sample format + - name: PZPrepareEstimatorLens # Prepare the p(z) estimator + classname: Inform_BPZ_lite + - name: PZEstimatorLens # Measure lens galaxy PDFs + classname: BPZ_lite + - name: TXMeanLensSelector # select objects for lens bins from the PDFs + - name: Inform_NZDirLens # Prepare the DIR method inputs for the lens sample + classname: Inform_NZDir + - name: Inform_NZDirSource # Prepare the DIR method inputs for the source sample + classname: Inform_NZDir + - name: PZRailSummarizeLens # Run the DIR method on the lens sample to find n(z) + classname: PZRailSummarize + - name: TXLensCatalogSplitter # Split the lens sample tomographically + - name: TXLSSWeightsUnit # add systematic weights to the lens sample (weight=1 for this example) + - name: TXTracerMetadata # collate metadata + - name: TXLensNoiseMaps # Compute lens noise using half-splits + - name: TXLensMaps # make source lens and n_gal maps + - name: TXDensityMaps # turn mask and ngal maps into overdensity maps + - name: TXTwoPointFourier # Compute power spectra C_ell + - name: TXTwoPoint # Compute real-space 2-point correlations + threads_per_process: 128 + - name: TXPureB # Compute pure BB + - name: TXJackknifeCenters # Split the area into jackknife regions + - name: TXSourceDiagnosticPlots # Make a suite of diagnostic plots + +# modules and packages to import that have pipeline +# stages defined in them +modules: > + txpipe + rail.creation.degradation.grid_selection + rail.creation.engines.flowEngine + rail.estimation.algos.NZDir + rail.estimation.algos.bpz_lite + +# where to find any modules that are not in this repo, +# and any other code we need. +python_paths: + - submodules/WLMassMap/python/desc/ + +# Where to put outputs +output_dir: data/des-y3/outputs + +# How to run the pipeline: mini, parsl, or cwl +launcher: + name: mini + interval: 1.0 + +# Where to run the pipeline: cori-interactive, cori-batch, or local +site: + name: local + max_threads: 128 + +# configuration settings +config: examples/desy3/config_Bmodes.yml + +# These are overall inputs to the whole pipeline, not generated within it +inputs: + # See README for paths to download these files + shear_catalog : /global/cfs/cdirs/lsst/groups/WL/projects/txpipe-sys-tests/des-y3/shear_catalog_desy3_unmasked_withfakez_v2.h5 + star_catalog : /global/cfs/cdirs/lsst/groups/WL/projects/txpipe-sys-tests/des-y3/DES_psf_y3_catalog.hdf5 + shear_photoz_stack : /global/cfs/cdirs/lsst/groups/WL/projects/txpipe-sys-tests/des-y3/shear_photoz_stack.hdf5 + random_cats : /global/cfs/cdirs/lsst/groups/WL/projects/txpipe-sys-tests/des-y3/randoms_desy3_RM.hdf5 + binned_random_catalog : /global/cfs/cdirs/lsst/groups/WL/projects/txpipe-sys-tests/des-y3/binned_randoms_desy3_RM.hdf5 + binned_random_catalog_sub: /global/cfs/cdirs/lsst/groups/WL/projects/txpipe-sys-tests/des-y3/binned_randoms_desy3_RM.hdf5 + fiducial_cosmology : data/fiducial_cosmology.yml + photometry_catalog : data/example/inputs/photometry_catalog.hdf5 + calibration_table : data/example/inputs/sample_cosmodc2_w10year_errors.dat + exposures : data/example/inputs/exposures.hdf5 + flow : data/example/inputs/example_flow.pkl + random_cats_source : Null + +# if supported by the launcher, restart the pipeline where it left off +# if interrupted +resume: True +# where to put output logs for individual stages +log_dir: data/des-y3/logs +# where to put an overall parsl pipeline log +pipeline_log: data/des-y3/log.txt + diff --git a/txpipe/data_types/types.py b/txpipe/data_types/types.py index 4fe3d879a..9c838bb56 100755 --- a/txpipe/data_types/types.py +++ b/txpipe/data_types/types.py @@ -528,7 +528,7 @@ def get_2d_n_of_z(self, zmax=3.0, nz=301): def get_bin_n_of_z(self, bin_index, zmax=3.0, nz=301): ensemble = self.read_ensemble() npdf = ensemble.npdf - if bin_index >= npdf - 1: + if bin_index > npdf - 1: raise ValueError(f"Requested n(z) for bin {bin_index} but only {npdf-1} bins available. For the 2D bin use get_2d_n_of_z.") z = np.linspace(0, zmax, nz) return z, ensemble.pdf(z)[bin_index] diff --git a/txpipe/diagnostics.py b/txpipe/diagnostics.py index f8de05ab2..6d75dd45c 100644 --- a/txpipe/diagnostics.py +++ b/txpipe/diagnostics.py @@ -300,8 +300,11 @@ def plot_psf_shear(self): fig.close() f = self.open_output("g_psf_g_out") - data =[mu1,mu2,mean11,mean12,mean21,mean22,std11,std12,std21,std22,line11,line12,line21,line22] - f.write(''.join([str(i) + '\n' for i in data])) + data = np.c_[mu1,mu2,mean11,mean12,mean21,mean22,std11,std12,std21,std22,line11,line12,line21,line22] + f.write('#[0]mu1 [1]mu2 [2]mean11 [3]mean12 [4]mean21 [5]mean22 [6]std11 [7]std12 [8]std21 [9]std22 [10]line11 [11]line12 [12]line21 [13]line22\n') + for row in data: + line = ' '.join(map(str, row)) + '\n' + f.write(line) f.close() def plot_psf_size_shear(self): @@ -372,9 +375,14 @@ def plot_psf_size_shear(self): fig.close() f = self.open_output("g_psf_T_out") - data =[mu,mean1,mean2,std1,std2,line1,line2] - f.write(''.join([str(i) + '\n' for i in data])) + f.write('#[0]mu [1]mean1 [2]mean2 [3]std1 [4]std2 [5]line1 [6]line2\n') + data = np.c_[mu,mean1,mean2,std1,std2,line1,line2] + for row in data: + line = ' '.join(map(str, row)) + '\n' + f.write(line) f.close() + + def plot_snr_shear(self): @@ -445,11 +453,17 @@ def plot_snr_shear(self): plt.tight_layout() fig.close() + f = self.open_output("g_snr_out") - data =[mu,mean1,mean2,std1,std2,line1,line2] - f.write(''.join([str(i) + '\n' for i in data])) + f.write('#[0]mu [1]mean1 [2]mean2 [3]std1 [4]std2 [5]line1 [6]line2\n') + data = np.c_[mu,mean1,mean2,std1,std2,line1,line2] + for row in data: + line = ' '.join(map(str, row)) + '\n' + f.write(line) f.close() - + + + def plot_size_shear(self): # mean shear in bins of galaxy size print("Making mean shear galaxy size plot") @@ -519,10 +533,15 @@ def plot_size_shear(self): plt.tight_layout() fig.close() + f = self.open_output("g_T_out") - data =[mu,mean1,mean2,std1,std2,line1,line2] - f.write(''.join([str(i) + '\n' for i in data])) + f.write('#[0]mu [1]mean1 [2]mean2 [3]std1 [4]std2 [5]line1 [6]line2\n') + data = np.c_[mu,mean1,mean2,std1,std2,line1,line2] + for row in data: + line = ' '.join(map(str, row)) + '\n' + f.write(line) f.close() + def plot_mag_shear(self): # mean shear in bins of magnitude diff --git a/txpipe/twopoint_null_tests.py b/txpipe/twopoint_null_tests.py index b1e8f1191..63260e881 100644 --- a/txpipe/twopoint_null_tests.py +++ b/txpipe/twopoint_null_tests.py @@ -8,6 +8,7 @@ PNGFile, TextFile, QPNOfZFile, + MapsFile ) import numpy as np from .twopoint import TXTwoPoint, SHEAR_SHEAR, SHEAR_POS, POS_POS @@ -766,3 +767,365 @@ def write_output(self, source_list, lens_list, meta, results): # Finally, save the output to Sacc file S.save_fits(self.get_output("aperture_mass_data"), overwrite=True) + + + +class TXShearBmode(PipelineStage): + ''' + Make shear B-mode measurements. + This stage computes xip and xim in very narrow angular bins and transforms them into + Fourier B-modes using narrow window functions. This allows us to compute B-modes without + being affected by (a) the mask (which is an issue if we use https://arxiv.org/abs/astro-ph/0511629, + which is implemented in NaMaster, since masks for LSS is more complicated than CMB) and + (b) noise bias (since the measurements we make here are real-space measurements). + ''' + name = "TXShearBmode" + parallel = False + + inputs = [ + ("shear_photoz_stack" , QPNOfZFile), + ("source_maps" , MapsFile), + ("mask" , MapsFile), + ("twopoint_data_real_raw", SACCFile), + ] + + outputs = [("twopoint_data_fourier_shearbmode", SACCFile), + ] + + config_options = { + 'method' : 'hybrideb', + 'purify_b' : False, + 'Nell' : 20, + 'lmin' : 200, + 'lmax' : 2000, + 'lspacing' : 'log', + 'bin_file' : "", + 'theta_min' : 2.5, + 'theta_max' : 250, + 'Ntheta' : 1000, + 'Nsim' : 1000 + } + + def run(self): + import os,sys + import pymaster as nmt + import healpy as hp + import sacc + import pyccl + import hybrideb + import datetime + import pickle + from tqdm import tqdm + + if self.config["method"] == 'namaster': + self.run_namaster(self.config["purify_b"]) + + elif self.config["method"] == 'hybrideb': + self.run_hybrideb() + + else: + raise ValueError("Method must 'namaster' or 'hybrideb'") + + + + def run_namaster(self,purify_b): + ''' + B-mode calculation that is already implemented in NaMaster + ''' + print('running namaster') + import pymaster as nmt + import healpy as hp + from tqdm import tqdm + import pickle + + + lmin = self.config['lmin'] + lmax = self.config['lmax'] + Nell = self.config['Nell'] + Nsims = self.config["Nsims"] + lspacing = self.config['lspacing'] + + if purify_b: + print("WARNING: Namaster's B-mode purification requires the mask to be heavily apodized.") + print("For a realistic LSS mask, this will most likely not work, but this function is ") + print("implemented for completeness nonetheless.") + + # Open source maps (g1,g2,weights) + with self.open_input("source_maps", wrapper=True) as f: + + # Get the number of tomographic bins + # +1 comes from also loading the non-tomographic sample + nbin_source = f.file["maps"].attrs["nbin_source"]+1 + + # Load the tomographic samples + g1_maps = [f.read_map(f"g1_{b}") for b in range(nbin_source-1)] + g2_maps = [f.read_map(f"g2_{b}") for b in range(nbin_source-1)] + + # Load non-tomographic sample + g1_maps.append(f.read_map(f"g1_2D")) + g2_maps.append(f.read_map(f"g2_2D")) + + if lmax>3*hp.npix2nside(len(g1_maps[0])): + raise ValueError("lmax must be smaller than 3*nside") + + # Open mask + with self.open_input("mask", wrapper=True) as f: + mask = f.read_map("mask") + + # Define binning + if lspacing=='lin': + bine = np.linspace(lmin+1, lmax+1, Nell+1, dtype=np.int64) + elif lspacing=='log': + bine = np.geomspace(lmin+1, lmax+1, Nell+1, dtype=np.int64) + else: + raise ValueError("lspacing must be either 'log' or 'lin'") + + b = nmt.NmtBin.from_edges(bine[:-1], bine[1:]) + + # Initialize the fields + fields = {} + for i in range(nbin_source): + g1_maps[i][g1_maps[i]==hp.UNSEEN]=0 + g2_maps[i][g2_maps[i]==hp.UNSEEN]=0 + fields[i] = nmt.NmtField(mask, [g1_maps[i], g2_maps[i]], purify_e=False, purify_b=purify_b,lmax=lmax) + + # To speed the covariance calculation up, one can use multiple nodes to do this calculation externally + products_dict = {'mask':mask,'Nell':Nell,'nbin_source':nbin_source, 'g1_maps': g1_maps, 'g2_maps': g2_maps} + + file_namaster_intermediate = self.get_output("twopoint_data_fourier_shearbmode")[:-5]+ f"_namaster.pkl" + + with open(file_namaster_intermediate, 'wb') as handle: + pickle.dump(products_dict, handle, protocol=pickle.HIGHEST_PROTOCOL) + + # Compute Cls and store them in dictionrary + ret=np.zeros((nbin_source,nbin_source,len(b.get_effective_ells()))) + + for zi in range(nbin_source): + for zj in range(zi,nbin_source): + + if zi!=zj and (zi==nbin_source-1 or zj==nbin_source-1): + # No need to take cross-correlation between tomographic and non-tomographic sample + continue + else: + field1 = fields[zi] + field2 = fields[zj] + w_yp = nmt.NmtWorkspace() + w_yp.compute_coupling_matrix(field1, field2, b) + + cl_coupled = nmt.compute_coupled_cell(field1,field2) + cl_decoupled = w_yp.decouple_cell(cl_coupled) + ret[zj,zi,:]= cl_decoupled[3] + + # Compute covariance by randomly rotating shear values in pixels. + tmparr = np.zeros((int(nbin_source*(nbin_source+1)/2*len(b.get_effective_ells())),Nsims)) + + for k in tqdm(range(Nsims)): + fields = {} + for i in range(nbin_source): + # Rotate shear maps randomly to create noise maps + idx = np.where(g1_maps[i]!=0)[0] + psi = np.random.uniform(0, 2 * np.pi,size=len(idx)) + g1_maps[i][idx] = g1_maps[i][idx]*np.cos(2*psi) + g2_maps[i][idx]*np.sin(2*psi) + g2_maps[i][idx] = -g1_maps[i][idx]*np.sin(2*psi) + g2_maps[i][idx]*np.cos(2*psi) + + fields[i] = nmt.NmtField(mask, [g1_maps[i], g2_maps[i]], purify_e=False, purify_b=purify_b, lmax=lmax) + + tmp = np.array([]) + for zi in range(nbin_source): + for zj in range(zi,nbin_source): + + field1 = fields[zi] + field2 = fields[zj] + w_yp = nmt.NmtWorkspace() + w_yp.compute_coupling_matrix(field1, field2, b) + + cl_coupled = nmt.compute_coupled_cell(field1,field2) + cl_decoupled = w_yp.decouple_cell(cl_coupled) + tmp = np.concatenate([tmp,cl_decoupled[3]]) + + tmparr[:,k] = tmp + + n_bins = b.get_n_bands() + bin_weights = np.zeros([n_bins, b.lmax+1]) + + for i in range(n_bins): + bin_weights[i, b.get_ell_list(i)] = b.get_weight_list(i) + + ell = b.get_effective_ells() + bin_weights = bin_weights + results = ret + cov = np.cov(tmparr) + + print('Saving ShearBmode Cls in sacc file') + self.save_power_spectra(nbin_source, ell, results, cov) + + + def run_hybrideb(self): + ''' + B-mode method of Becker and Rozo 2015 http://arxiv.org/abs/1412.3851 + ''' + import hybrideb + import pickle + import sacc + import os + + print('running hybrideb') + Nell = self.config['Nell'] + Nsims = self.config["Nsims"] + Ntheta = self.config["Ntheta"] + + # Check if nbins is less than 1000, and throw warning. + if self.config['Ntheta']<1000: + print("WARNING: Calculating hybridEB using Ntheta<1000, which may lead to inaccurate results.") + + # Computing the weights takes a few minutes so its a lot faster + # to precompute them and load them again in subsequent runs. + file_precomputed_weights = self.get_output("twopoint_data_fourier_pureB")+ f".precomputedweights.pkl" + + if os.path.exists(file_precomputed_weights): + print(f"{self.rank} WARNING: Using precomputed weights from revious run: {file_precomputed_weights}") + with open(file_precomputed_weights, "rb") as f: + geb_dict = pickle.load(f) + + else: + # To run this method we need to first compute the Fourier and real-space windows + # hybrideb.HybridEB: Sets up the estimator by creating top hat filters in real-space + # hybrideb.BinEB : Calculates the intermediate filter quantities, namely fa, fb, M+, M- [Eqs 6-11] + # hybrideb.GaussEB : Creates Gaussian windows in Fourier space + heb = hybrideb.HybridEB(self.config['theta_min'], self.config['theta_max'], Ntheta) + beb = hybrideb.BinEB(self.config['theta_min'], self.config['theta_max'], Ntheta) + geb = hybrideb.GaussEB(beb, heb) + + # geb is in a fromat that can not be naturally saved so we convert it to a regular dictionary. + geb_dict = {} + for i in range(0,Nell): + geb_dict[f"{i+1}"] = {} + for j in range(0,6): + geb_dict[f"{i+1}"][f"{j+1}"]=geb(i)[j] + + # geb_dict stores the following information in a regular dictionary format + # 1: theta in radians + # 2: Fp plus component of the real-space window function + # 3: Fm minus component of the real-space window function + # 4: ell + # 5: Wp plus component of Fourier-space window function + # 6: Wm minus component of Fourier-space window function + + # 2 and 3 are used to get Xp/Xm = sum((Fp*xip_hat +/- Fm*xim_hat)/2) [eq1] + # 5 and 6 are used to get / = \int ell factors(ell) (Wp*Pe + Wm*Pb) [eq5] + + # Save this for subsequent runs since the windows don't change. + + with open(file_precomputed_weights, 'wb') as handle: + pickle.dump(geb_dict, handle, protocol=pickle.HIGHEST_PROTOCOL) + + En0 = np.zeros(Nell) + Bn0 = np.zeros(Nell) + En = np.zeros((Nell,Nsims)) + Bn = np.zeros((Nell,Nsims)) + ell = np.zeros(Nell) + + # Load xip/xim that were computed in the TXTwopoint stage and check length required for covariance + filename = self.get_input("twopoint_data_real_raw") + data_twopt = sacc.Sacc.load_fits(filename) + Qxip = sacc.standard_types.galaxy_shear_xi_plus + + source_tracers = set() + for b1, b2 in data_twopt.get_tracer_combinations(Qxip): + source_tracers.add(b1) + source_tracers.add(b2) + + nbin_source = max(len(source_tracers), 1) + + corvarr = np.zeros((int(nbin_source*(nbin_source+1)/2*Nell), Nsims)) + results = np.zeros((nbin_source,nbin_source,Nell)) + + # Loop over all bin combinations + c = 0 + for zi in range(0,nbin_source): + for zj in range(zi,nbin_source): + + # Initial and final row indices of the resulting array (covarr) to fill. + # For each tomographic bin we have Nell bins to fill. + ii = Nell*c + ff = Nell*(c+1) + + tmp = data_twopt.copy() + + # From the sacc file only load relevant tracer combinations + tmp.keep_selection(tracers=(f'source_{zj}', f'source_{zi}')) + dvec = tmp.mean + xip = dvec[:int(Ntheta)] + xim = dvec[int(Ntheta):] + + # Make random draws based on mean and covariance of xip/xim measurements + x = np.random.multivariate_normal(mean = tmp.mean, cov =tmp.covariance.covmat , size = Nsims) + Rxip = x[:,:int(Ntheta)] + Rxim = x[:,int(Ntheta):] + + # Convert each random draw into B-mode measurement and compute covariance + for n in range(Nsims): + for i in range(int(Nell)): + res = geb_dict[f"{i+1}"] + Fp = res['2'] + Fm = res['3'] + En[i,n] = np.sum(Fp*Rxip[n,:] + Fm*Rxim[n,:])/2 + Bn[i,n] = np.sum(Fp*Rxip[n,:] - Fm*Rxim[n,:])/2 + + corvarr[ii:ff,:] = Bn + + # Compute actual data vector + for i in range(int(Nell)): + res = geb_dict[f"{i+1}"] + Fp = res['2'] + Fm = res['3'] + En0[i] = np.sum(Fp*xip + Fm*xim)/2 + Bn0[i] = np.sum(Fp*xip - Fm*xim)/2 + ell[i] = res['4'][np.argmax(res['5'])] # setting ell to the peak of the Gaussian window function + + results[zj,zi,:]= Bn0[:] + c+=1 + + cov = np.cov(corvarr) + + print('Saving pureB Cls in sacc file') + self.save_power_spectra(nbin_source, ell, results, cov) + + + def save_power_spectra(self, nbin_source, ell, results, cov): + import sacc + import datetime + + S = sacc.Sacc() + S.metadata['nbin_source'] = nbin_source + S.metadata['creation'] = datetime.datetime.now().isoformat() + S.metadata['method'] = self.config['method'] + S.metadata['info'] = 'ClBB' + + CBB = sacc.standard_types.galaxy_shear_cl_bb + + with self.open_input("shear_photoz_stack", wrapper=True) as f: + for i in range(nbin_source): + if i==nbin_source-1: + z,Nz = f.get_2d_n_of_z() + S.add_tracer("NZ", f"source_{i}", z, Nz) + else: + z, Nz = f.get_bin_n_of_z(i) + S.add_tracer("NZ", f"source_{i}", z, Nz) + + for zi in range(0,nbin_source): + for zj in range(zi,nbin_source): + tracer1 = f"source_{zj}" + tracer2 = f"source_{zi}" + val = results[zj,zi,:] + print(val) + ell = ell + for k in range(0,len(val)): + S.add_data_point(CBB, (tracer1, tracer2), value=val[k], ell=ell[k]) + + S.add_covariance(cov) + + output_filename = self.get_output("twopoint_data_fourier_shearbmode") + S.save_fits(output_filename, overwrite=True) + +