diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml new file mode 100644 index 0000000..d3b6b70 --- /dev/null +++ b/.github/workflows/CI.yml @@ -0,0 +1,56 @@ +# Author - Stuart McAlpine - stuart.mcalpine@fysik.su.se - Jan 2023 +name: DESC Legacy Blinding for Cosmosis + +# How does the workflow get triggered? +on: + # Triggers when push/pull-request made to the main branch. + pull_request: + branches: + - main + push: + branches: + - main + - develop + +# List of jobs for this workflow. + +jobs: + # Our pytest job. + legacy_blinding-pytest: + # Our strategy lists the OS and Python versions we want to test on. + strategy: + # Don't quit all jobs if only one job fails. + fail-fast: false + matrix: + python-version: ["3.9"] + os: [ubuntu-20.04, macos-latest] + include: + - os: ubuntu-20.04 + INSTALL_DEPS: sudo apt-get update && sudo apt-get -y install gfortran-7 swig libopenmpi-dev openmpi-bin libopenblas-dev && sudo ln -s `which gfortran-7` /usr/local/bin/gfortran + - os: macos-latest + INSTALL_DEPS: brew update-reset && HOMEBREW_NO_AUTO_UPDATE=1 HOMEBREW_NO_INSTALLED_DEPENDENTS_CHECK=1 brew install gcc swig libomp open-mpi openblas && if [ ! -f /usr/local/bin/gfortran ]; then ln -s /usr/local/bin/gfortran-$(brew list --versions gcc | awk '{print $2}' | cut -d. -f1) /usr/local/bin/gfortran; fi + + # What operating system is this job running on? + runs-on: ${{ matrix.os }} + + # Our CI steps for this job. + steps: + # Check out this repository code. + - name: Check out repository code + uses: actions/checkout@v3 + + # Install Python. + - name: Set up Python ${{ matrix.python-version }} + uses: actions/setup-python@v4 + with: + python-version: ${{ matrix.python-version }} + + - name: Install dependencies + run: ${{ matrix.INSTALL_DEPS }} + + # Install our package. + - name: Install legacy_blinding + run: | + python -m pip install --upgrade pip + python -m pip install . + pytest ./tests diff --git a/.gitignore b/.gitignore index 68bc17f..72b15a5 100644 --- a/.gitignore +++ b/.gitignore @@ -2,6 +2,7 @@ __pycache__/ *.py[cod] *$py.class +*.DS_Store # C extensions *.so diff --git a/README.md b/README.md index 7cf7a0a..227688a 100644 --- a/README.md +++ b/README.md @@ -1,2 +1,40 @@ -# legacy_blinding -Repository with Legacy Cosmology Blinding for 3x2pt data-vector blinding using Cosmosis v2 +# Legacy Cosmosis Blinding for 2pt Data-Vectors using Cosmosis + +This repository contains an implementation of [Muir et al. 2020](https://arxiv.org/abs/1911.05929) data-vector blinding strategy using Cosmosis V2. +The package in this repository is a direct adaptation of Jessie Muir's DES Y3 blinding scripts to the DESC context. + +For the DESC Blinding Library, see [here](https://github.com/LSSTDESC/Blinding). + +## Installation +- Install cosmosis via `conda install -c conda-forge cosmosis` +- Activate the cosmosis configuration via `source cosmosis-configure` +- Install the cosmosis standard library via `conda install -c conda-forge cosmosis-build-standard-library` +- Install the standard library in the correct cosmosis folder `cosmosis-build-standard-library -i` + +### Legacy Blinding installation +Simply install with +```bash +python -m pip install -e . +``` + +### Test the installation +To run the unit tests: +```bash +cd tests +pytest . +``` + +## Usage +> Don't forget to `source cosmosis-configure` + +To get the docstring on the command-line interface: +``` +python -m blind_2pt_cosmosis --help +``` + +To run an example: +```bash +python -m blind_2pt_cosmosis -u src/blind_2pt_cosmosis/cosmosis_files/sim_fiducial.fits [--log-level DEBUG] +``` + +[!] Note: You can also include all the command line arguments in a `.config` file and call `python -m blind_2pt_cosmosis @example.config` instead. This will be adapted to use `yaml` files in the near future. \ No newline at end of file diff --git a/pyproject.toml b/pyproject.toml new file mode 100644 index 0000000..c277cf7 --- /dev/null +++ b/pyproject.toml @@ -0,0 +1,37 @@ +[build-system] +requires = ["setuptools", "wheel"] +build-backend = "setuptools.build_meta" + +[project] +name = "legacy_blinding" +description = "Legacy Cosmosis Blinding for DESC. This is a Cosmosis V2 port of Jessie Muir's legacy_blinding scripts for DES-Y3." +readme = "README.md" +requires-python = ">=3.8" +license = {file = "LICENSE"} +maintainers = [ + {name = "Arthur Loureiro", email = "arthur.loureiro@fysik.su.se"}, + {name = "Jessie Muir", email = "jmuir@perimeterinstitute.ca"}, +] +classifiers = [ + "Programming Language :: Python :: 3", + "License :: OSI Approved :: MIT License", + "Operating System :: OS Independent", +] +dependencies = [ + "numpy>=1.24.0", + "scipy>=1.9.0", + "cosmosis>=2.5", + "astropy>=5.1", + "pytest>=7.2.5", + "argparse>=1.0", + "sacc>=0.12", +] +keywords = ["desc", "python", "blinding", "cosmosis"] +dynamic = ["version"] + +[tool.setuptools.packages.find] +where = ["src"] + +[tool.setuptools.package-data] +"blind_2pt_cosmosis" = ["cosmosis_files/*.txt", "cosmosis_files/*.ini"] + diff --git a/src/.DS_Store b/src/.DS_Store new file mode 100644 index 0000000..2e9041f Binary files /dev/null and b/src/.DS_Store differ diff --git a/src/blind_2pt_cosmosis/.DS_Store b/src/blind_2pt_cosmosis/.DS_Store new file mode 100644 index 0000000..778d7fa Binary files /dev/null and b/src/blind_2pt_cosmosis/.DS_Store differ diff --git a/src/blind_2pt_cosmosis/__init__.py b/src/blind_2pt_cosmosis/__init__.py new file mode 100644 index 0000000..d4def22 --- /dev/null +++ b/src/blind_2pt_cosmosis/__init__.py @@ -0,0 +1,2 @@ +# Author: Arthur Loureiro & Jessie Muir +from ._version import __version__ \ No newline at end of file diff --git a/src/blind_2pt_cosmosis/__main__.py b/src/blind_2pt_cosmosis/__main__.py new file mode 100644 index 0000000..2174bd4 --- /dev/null +++ b/src/blind_2pt_cosmosis/__main__.py @@ -0,0 +1,49 @@ +import logging +from .param_shifts import draw_flat_param_shift, get_factordict +from .io import get_parser, get_stored_seed_and_tag +from .run_cosmosis_2pt import run_cosmosis_togen_2ptdict +from .twopt_utils import apply_2pt_blinding_and_save_fits + + +def main(): + """ + main function to be called from command line + """ + # gets the parser from the io module + parser = get_parser() + args = parser.parse_args() + + # Configure the logger level based on the input argument + log_level = getattr(logging, args.log_level) + logging.basicConfig(level=log_level, format='%(asctime)s - %(name)s - %(levelname)s - %(message)s') + + # Create a logger for the __main__ module + logger = logging.getLogger("2pt_blinding") + logger.debug(args) + + params_shifts = draw_flat_param_shift(args.seed, args.paramshifts) + logger.debug(f"Parameters shifts are: {params_shifts}") + + #get blinding factors + reff_dict = run_cosmosis_togen_2ptdict(inifile=args.ini) #args.ini is a template + logger.debug(f"Calculated Reference Dict") + # FIXME: Add nz_file and angles_file to args + shift_dict = run_cosmosis_togen_2ptdict(inifile=args.ini, pdict=params_shifts, + nz_file=None, angles_file=None) + logger.debug(f"Calculated Shifted Dict") + #gets the blinding factors data vectors + factor_dict = get_factordict(reff_dict, shift_dict, bftype = args.bftype) + logger.debug(f"Calculated Blinded Data Vectors") + + # Gets some I/O information + storeseed, tagstr = get_stored_seed_and_tag(args) + + # applies the shift to the 2pt data-vector: + #FIXME: This function is also saving the file, we want to split it! + apply_2pt_blinding_and_save_fits(factor_dict, origfitsfile=args.origfits, + outfname=args.outfname, outftag=tagstr, + bftype=args.bftype, storeseed = storeseed) + + +if __name__ == '__main__': + main() diff --git a/src/blind_2pt_cosmosis/_version.py b/src/blind_2pt_cosmosis/_version.py new file mode 100644 index 0000000..d368745 --- /dev/null +++ b/src/blind_2pt_cosmosis/_version.py @@ -0,0 +1 @@ +__version__ = "2023.10.11" \ No newline at end of file diff --git a/src/blind_2pt_cosmosis/cosmosis_files/blinding_values_template.ini b/src/blind_2pt_cosmosis/cosmosis_files/blinding_values_template.ini new file mode 100644 index 0000000..354243a --- /dev/null +++ b/src/blind_2pt_cosmosis/cosmosis_files/blinding_values_template.ini @@ -0,0 +1,69 @@ +[cosmological_parameters] +omega_m = 0.1 0.3 0.9 +h0 = 0.55 0.69 0.91 +omega_b = 0.03 0.048 0.12 +n_s = 0.87 0.97 1.07 +sigma_8_input = 0.2 0.843 1.6 +A_s = 0.5e-09 2.19e-9 5.0e-09 +; omnuh2 = 0.0006 0.00083 0.01 +omnuh2 = 0.0 +massive_nu = 0 +massless_nu = 3.046 +omega_k = 0.0 +tau = 0.0697186 +w = -1.0 +;Helium mass fraction. Needed for Planck +yhe = 0.245341 + +; %include ${METHODSDIR}/cosmosis/de_models/${DEMODEL}.ini + +[shear_calibration_parameters] +m1 = -0.1 0.0 0.1 +m2 = -0.1 0.0 0.1 +m3 = -0.1 0.0 0.1 +m4 = -0.1 0.0 0.1 + +[wl_photoz_errors] +bias_1 = -0.1 0.0 0.1 +bias_2 = -0.1 0.0 0.1 +bias_3 = -0.1 0.0 0.1 +bias_4 = -0.1 0.0 0.1 + +[lens_photoz_errors] +bias_1 = -0.05 0.0 0.05 +bias_2 = -0.05 0.0 0.05 +bias_3 = -0.05 0.0 0.05 +bias_4 = -0.05 0.0 0.05 +bias_5 = -0.05 0.0 0.05 + +[bias_lens] +b1 = 0.8 1.7 3.0 +b2 = 0.8 1.7 3.0 +b3 = 0.8 1.7 3.0 +b4 = 0.8 2.0 3.0 +b5 = 0.8 2.0 3.0 + +[mag_alpha_lens] +; alpha_1 = -4. 0.903125 4. +; alpha_2 = -4. 0.68572965 4. +; alpha_3 = -4. 0.65340057 4. +; alpha_4 = -4. 1.58867862 4. +; alpha_5 = -4. 1.93754879 4. +alpha_1 = 0.903125 +alpha_2 = 0.68572965 +alpha_3 = 0.65340057 +alpha_4 = 1.58867862 +alpha_5 = 1.93754879 + +[intrinsic_alignment_parameters] +z_piv = 0.62 +A1 = -5.0 0.7 5.0 +A2 = -5.0 -1.36 5.0 +alpha1=-5.0 -1.7 5.0 +alpha2=-5.0 -2.5 5.0 +bias_ta = 0.0 1.0 2.0 +#bias_tt = 1.0 + + + + diff --git a/src/blind_2pt_cosmosis/cosmosis_files/default_blinding_template.ini b/src/blind_2pt_cosmosis/cosmosis_files/default_blinding_template.ini new file mode 100644 index 0000000..3799488 --- /dev/null +++ b/src/blind_2pt_cosmosis/cosmosis_files/default_blinding_template.ini @@ -0,0 +1,51 @@ +; ############################################################ +; This is a template file that is referenced by blind_2pt_usingcosmosis.py +; It is used as a template for running cosmosis to generate +; blinding factors to be applied to 2pt data vectors. To adapt +; for a specific pipeline: +; 1. Edit the %include params.ini line to point towards your baseline pipeline file. +; 2. The values file should be edited to match your fiducial run. Since we want +; to blind sigma8, sigma8_input needs to be in the cosmological parameters, +; so it may need to be a values file written specifically for blinding. You could +; always set it up to %include your standard values file though. +; >> the start values for each parameter is what will be used as reference vals +; when generating the blinding factors. +; 3. Make sure the file in 2PT_FILE points to a file with the format expected by +; your pipeline (with the same binning, etc.) +; 4. The modules list in [pipeline] should match your standard pipeline with a one +; addition: you need sigma8_rescale to be able to include sigma8_input as a parameter +; ############################################################ + + +; import baseline parameter file; change this to point towards the params.ini file for +; whatever main pipeline you want to mimic +%include ./src/blind_2pt_cosmosis/cosmosis_files/params.ini + +[DEFAULT] +; this 2PT file is where we'll pull the n(z) info for computing blinding factors +2PT_FILE = ./src/blind_2pt_cosmosis/cosmosis_files/sim_fiducial.fits +;2PT_FILE = pipeline/2pt_sim_1110_baseline_Y3cov.fits +;should this be my 2pt measurement DV? + +[pipeline] +; THIS MODULE LIST SHOULD MATCH YOUR SETTINGS FOR YOUR STANDARD ANALYSIS PIPELINE +quiet=T +timing=F +debug=F +priors = ./src/blind_2pt_cosmosis/cosmosis_files/priors.ini +; sigma8_rescale included so we can set sigma8 as input +modules = consistency camb extrapolate fits_nz lens_photoz_bias fast_pt source_photoz_bias unbiased_galaxies IA pk_to_cl_gg pk_to_cl add_magnification add_intrinsic 2pt_shear 2pt_gal 2pt_gal_shear shear_m_bias 2pt_like +; if you want to blind sigma_8, set up values file to have sigma8_input as a parameter rather than A_s +values = ./src/blind_2pt_cosmosis/cosmosis_files/blinding_values_template.ini + + +[sigma8_rescale] +file = cosmosis-standard-library/utility/sample_sigma8/sigma8_rescale.py + +[test] +; if the sampler happens to be set to test, make sure it doesn't +; save any output when we do blinding. We don't want to see the cosmology +; and 2pt functions that go into making the blinding factor! +save_dir= + + diff --git a/src/blind_2pt_cosmosis/cosmosis_files/params.ini b/src/blind_2pt_cosmosis/cosmosis_files/params.ini new file mode 100644 index 0000000..06b575b --- /dev/null +++ b/src/blind_2pt_cosmosis/cosmosis_files/params.ini @@ -0,0 +1,311 @@ +;params.ini + +[runtime] +sampler = test +root = ${COSMOSIS_SRC_DIR} +; pre_script=./pre_script.sh + +[DEFAULT] +BASELINE_DIR=. +DATAVECTOR_DIR=. +2PT_FILE=./sim_fiducial.fits +2PT_DATA_SETS = xip xim gammat wtheta +RUN_NAME = fiducial +; use as the last item in the pipeline to switch between 2pt_like and save_2pt +ACTION = save_2pt +DATAFILE = fiducial + + +[multinest] +max_iterations=50000 +multinest_outfile_root=%(RUN_NAME)s/mn_${DATAFILE} +resume=T +# Testing with Posteriors only: These settings are very fast (~6h, 128 cores) +# and produce reasonably good posteriors, but completely unreliable evidences. +# Useful for very preliminary results and debugging: +live_points=250 +efficiency=0.3 +tolerance=0.01 +constant_efficiency=F +#this is usually known as the "suggested standard run" +#(should give decent results for posteriors and evidence in a Y1KP-like analysis) +#depending on your applications, it could be inappropriate +; live_points=500 +; efficiency=0.3 +; tolerance=0.1 +; constant_efficiency=F + +[output] +filename= %(RUN_NAME)s/chain_%(RUN_NAME)s_${DATAFILE}.txt +format=text +lock=F +privacy=F + +[emcee] +walkers = 100 +samples = 1000 +nsteps = 5 + +[star] +nsample_dimension = 50 + + +[test] +save_dir=%(RUN_NAME)s/test_%(RUN_NAME)s +fatal_errors=T + +[pipeline] +modules = consistency camb extrapolate fits_nz lens_photoz_bias fast_pt source_photoz_bias unbiased_galaxies IA pk_to_cl_gg pk_to_cl add_magnification add_intrinsic 2pt_shear 2pt_gal 2pt_gal_shear shear_m_bias %(ACTION)s + +quiet=T +timing=F +debug=F +priors = priors.ini +values = values.ini +extra_output = cosmological_parameters/sigma_8 + +; enable this if your sampler supports it (MN does not) +fast_slow=F +first_fast_module=shear_m_bias + + +[bin_bias] +file = cosmosis-standard-library/bias/binwise_bias/bin_bias.py +perbin=T + +[consistency] +file = cosmosis-standard-library/utility/consistency/consistency_interface.py + +[camb] +file = cosmosis-standard-library/boltzmann/camb/camb_interface.py +mode=all +lmax=2500 +feedback=0 +kmin=1e-5 +kmax=10.0 +nk=400 + +;[halofit] +;file = cosmosis-standard-library/boltzmann/halofit_takahashi/halofit_interface.so +;nk=700 + +[extrapolate] +file = cosmosis-standard-library/boltzmann/extrapolate/extrapolate_power.py +kmax = 500. + +[fits_nz] +file = cosmosis-standard-library/number_density/load_nz_fits/load_nz_fits.py +nz_file = %(2PT_FILE)s +data_sets = lens source +prefix_section = T +prefix_extension = T + +#the two modules below might be replaced in the future by a different parametrization +[lens_photoz_bias] +file = cosmosis-standard-library/number_density/photoz_bias/photoz_bias.py +mode = additive +sample = nz_lens +bias_section = lens_photoz_errors +interpolation = linear + +[source_photoz_bias] +file = cosmosis-standard-library/number_density/photoz_bias/photoz_bias.py +mode = additive +sample = nz_source +bias_section = wl_photoz_errors +interpolation = linear + +[unbiased_galaxies] +file = cosmosis-standard-library/bias/no_bias/no_bias.py +use_lin_power = False + + +[IA_old] +#linear_alignments_interface.py +file=cosmosis-standard-library/intrinsic_alignments/la_model/linear_alignments_interface.py +method=bk_corrected +do_galaxy_intrinsic=T + +[IA] +file=cosmosis-standard-library/intrinsic_alignments/tatt/tatt_interface.py +mode=all +sub_lowk=F +sub_const=F +include_s2_terms=F +do_galaxy_intrinsic=F +ia_model=tatt +Asigma8=F + + +[ia_z_field] +file = ${COSMOSIS_SRC_DIR}/cosmosis-standard-library/intrinsic_alignments/z_powerlaw/ia_z_powerlaw.py +do_galaxy_intrinsic = T + +[pk_to_cl] +file = cosmosis-standard-library/structure/projection/project_2d.py +ell_min_logspaced = 0.1 +ell_max_logspaced = 5.0e5 +; n_ell_logspaced = 400 +n_ell_logspaced = 100 +shear-shear = source-source +shear-intrinsic = source-source +intrinsic-intrinsic = source-source +intrinsicb-intrinsicb=source-source +lingal-shear = lens-source +lingal-intrinsic = lens-source +lingal-magnification = lens-lens +magnification-shear = lens-source +magnification-magnification = lens-lens +magnification-intrinsic = lens-source + +verbose = F +get_kernel_peaks = F +sig_over_dchi = 20. +shear_kernel_dchi = 10. + +[pk_to_cl_gg] +file = cosmosis-standard-library/structure/projection/project_2d.py +lingal-lingal = lens-lens +do_exact = lingal-lingal +do_rsd = True +ell_min_linspaced = 1 +; ell_min_linspaced=1 +ell_max_linspaced = 4 +; ell_max_linspaced=9 +n_ell_linspaced = 5 +; n_ell_linspaced=10 +; ell_min_logspaced=10. +; ell_max_logspaced=1.e5 +ell_min_logspaced = 5. +ell_max_logspaced = 5.e5 +n_ell_logspaced = 80 +; n_ell_logspaced=150 +limber_ell_start = 200 +; limber_ell_start=201 +ell_max_logspaced=1.e5 +auto_only=lingal-lingal +sig_over_dchi_exact = 3.5 + +[add_magnification] +file = cosmosis-standard-library/structure/magnification/add_magnification.py +include_intrinsic=T + +[add_intrinsic] +file=cosmosis-standard-library/shear/add_intrinsic/add_intrinsic.py +shear-shear=T +position-shear=T +perbin=F + +;[add_eb] +;file = cosmosis-des-library/IAs/add_bmode_cl/add_bmode_cl.py + + +[choose_xip] +file = cosmosis-standard-library/utility/copy/copy_section.py +# copy shear_xi_eplusb_plus->shear_xi_plus +# and shear_xi_eminusb_minus->shear_xi_minus +;source = shear_xi_eplusb_plus shear_xi_eminusb_minus +;source = xiplus_i_j ximinus_i_j +dest = shear_xi_plus shear_xi_minus + +[fast_pt] +file = cosmosis-standard-library/structure/fast_pt/fast_pt_interface.py +do_ia = T +k_res_fac = 0.5 +verbose = F + +[shear_m_bias] +file = cosmosis-standard-library/shear/shear_bias/shear_m_bias.py +m_per_bin = True +; Despite the parameter name, this can operate on xi as well as C_ell. +cl_section = shear_xi_plus shear_xi_minus +cross_section = galaxy_shear_xi +verbose = F + +;[shear_2pt_eplusb] +;file = cosmosis-standard-library/shear/cl_to_xi_fullsky/cl_to_xi_interface.py +;ell_max = 40000 +;xi_type='22' +;theta_file=%(2PT_FILE)s +;bin_avg = T +;input_section_name = shear_cl_eplusb +;output_section_name = shear_xi_eplusb + +;[shear_2pt_eminusb] +;file = cosmosis-standard-library/shear/cl_to_xi_fullsky/cl_to_xi_interface.py +;ell_max = 40000 +;xi_type='22' +;theta_file=%(2PT_FILE)s +;bin_avg = T +;input_section_name = shear_cl_eminusb +;output_section_name = shear_xi_eminusb + +[2pt_shear] +file = cosmosis-standard-library/shear/cl_to_xi_fullsky/cl_to_xi_interface.py +ell_max = 40000 +xi_type = EB +theta_file=%(2PT_FILE)s +bin_avg = T +; these get +input_section_name = shear_cl shear_cl_bb +output_section_name = shear_xi_plus shear_xi_minus + +[2pt_gal] +file = cosmosis-standard-library/shear/cl_to_xi_fullsky/cl_to_xi_interface.py +ell_max = 40000 +xi_type='00' +theta_file=%(2PT_FILE)s +bin_avg = T + +[2pt_gal_shear] +file = cosmosis-standard-library/shear/cl_to_xi_fullsky/cl_to_xi_interface.py +ell_max = 40000 +xi_type='02' +theta_file=%(2PT_FILE)s +bin_avg = T + + +[2pt_like] +file = cosmosis-standard-library/likelihood/2pt/2pt_like.py +do_pm_marg = True +data_file = %(2PT_FILE)s +data_sets = %(2PT_DATA_SETS)s +make_covariance=F +covmat_name=COVMAT +;%include scale_cuts/${SCALE_CUTS} +;%include ../y3-3x2pt-methods/cosmosis/scale_cuts/scales_3x2pt_Y1.ini +;%include ${HOME}/workspace/des-repos/y3-3x2pt-methods/cosmosis/scale_cuts/scales_3x2pt_Y1.ini + +[save_2pt] +file = cosmosis-standard-library/likelihood/2pt/save_2pt.py +theta_min = 2.5 +theta_max = 250.0 +n_theta = 20 +real_space = T +make_covariance = F +shear_nz_name = nz_source +position_nz_name = nz_lens +filename = ./sim_%(RUN_NAME)s.fits +overwrite = T +auto_only = galaxy_xi +;cut_wtheta = 1,2 1,3 2,3 1,4 2,4 3,4 1,5 2,5 3,5 4,5 +spectrum_sections = shear_xi_plus shear_xi_minus galaxy_shear_xi galaxy_xi +output_extensions = xip xim gammat wtheta +two_thirds_midpoint = T +copy_covariance=./sim_fiducial.fits + + +[sim_fits_nz_lens] +file = cosmosis-standard-library/number_density/load_nz_fits/load_nz_fits.py +;nz_file = nz_inputs/nz_y3_redmagic_v6.4.22_v2_gold_2.2.1_combined_max_bin_edges.fits +nz_file = ./cosmosis_files/sim_fiducial.fits +data_sets = lens +prefix_section = T +prefix_extension = T + +[sim_fits_nz_source] +file = cosmosis-standard-library/number_density/load_nz_fits/load_nz_fits.py +nz_file = nz_inputs/barcelona_lens1.fits +data_sets = source +prefix_section = T +prefix_extension = T \ No newline at end of file diff --git a/src/blind_2pt_cosmosis/cosmosis_files/priors.ini b/src/blind_2pt_cosmosis/cosmosis_files/priors.ini new file mode 100644 index 0000000..fe6eec1 --- /dev/null +++ b/src/blind_2pt_cosmosis/cosmosis_files/priors.ini @@ -0,0 +1,29 @@ +; These lines specify Gaussian priors +; for the various parameters in the various sections, +; with the specified mean and std. dev + +[wl_photoz_errors] +bias_1 = gaussian -0.0037 0.0177 +bias_2 = gaussian -0.0171 0.0150 +bias_3 = gaussian 0.0200 0.0138 +bias_4 = gaussian 0.0224 0.0215 + + +[shear_calibration_parameters] +m1 = gaussian 0.0 0.021 +m2 = gaussian 0.0 0.021 +m3 = gaussian 0.0 0.021 +m4 = gaussian 0.0 0.021 + + +[lens_photoz_errors] +bias_1 = gaussian 0.0 0.01 +bias_2 = gaussian 0.0 0.01 +bias_3 = gaussian 0.0 0.01 +bias_4 = gaussian 0.0 0.01 +bias_5 = gaussian 0.0 0.01 + + +[planck] +a_planck = gaussian 1.0 0.0025 + diff --git a/src/blind_2pt_cosmosis/cosmosis_files/sim_fiducial.fits b/src/blind_2pt_cosmosis/cosmosis_files/sim_fiducial.fits new file mode 100644 index 0000000..57c6c83 Binary files /dev/null and b/src/blind_2pt_cosmosis/cosmosis_files/sim_fiducial.fits differ diff --git a/src/blind_2pt_cosmosis/cosmosis_files/sim_fiducial.fits_cov.png b/src/blind_2pt_cosmosis/cosmosis_files/sim_fiducial.fits_cov.png new file mode 100644 index 0000000..c053b7d Binary files /dev/null and b/src/blind_2pt_cosmosis/cosmosis_files/sim_fiducial.fits_cov.png differ diff --git a/src/blind_2pt_cosmosis/cosmosis_files/sim_fiducial.fits_gammat.png b/src/blind_2pt_cosmosis/cosmosis_files/sim_fiducial.fits_gammat.png new file mode 100644 index 0000000..6265db0 Binary files /dev/null and b/src/blind_2pt_cosmosis/cosmosis_files/sim_fiducial.fits_gammat.png differ diff --git a/src/blind_2pt_cosmosis/cosmosis_files/sim_fiducial.fits_nz_lens.png b/src/blind_2pt_cosmosis/cosmosis_files/sim_fiducial.fits_nz_lens.png new file mode 100644 index 0000000..6491040 Binary files /dev/null and b/src/blind_2pt_cosmosis/cosmosis_files/sim_fiducial.fits_nz_lens.png differ diff --git a/src/blind_2pt_cosmosis/cosmosis_files/sim_fiducial.fits_nz_source.png b/src/blind_2pt_cosmosis/cosmosis_files/sim_fiducial.fits_nz_source.png new file mode 100644 index 0000000..9c776e0 Binary files /dev/null and b/src/blind_2pt_cosmosis/cosmosis_files/sim_fiducial.fits_nz_source.png differ diff --git a/src/blind_2pt_cosmosis/cosmosis_files/sim_fiducial.fits_wtheta.png b/src/blind_2pt_cosmosis/cosmosis_files/sim_fiducial.fits_wtheta.png new file mode 100644 index 0000000..94783ef Binary files /dev/null and b/src/blind_2pt_cosmosis/cosmosis_files/sim_fiducial.fits_wtheta.png differ diff --git a/src/blind_2pt_cosmosis/cosmosis_files/sim_fiducial.fits_xim.png b/src/blind_2pt_cosmosis/cosmosis_files/sim_fiducial.fits_xim.png new file mode 100644 index 0000000..bbb6578 Binary files /dev/null and b/src/blind_2pt_cosmosis/cosmosis_files/sim_fiducial.fits_xim.png differ diff --git a/src/blind_2pt_cosmosis/cosmosis_files/sim_fiducial.fits_xip.png b/src/blind_2pt_cosmosis/cosmosis_files/sim_fiducial.fits_xip.png new file mode 100644 index 0000000..80ce03d Binary files /dev/null and b/src/blind_2pt_cosmosis/cosmosis_files/sim_fiducial.fits_xip.png differ diff --git a/src/blind_2pt_cosmosis/cosmosis_files/sim_fiducial_BLINDED.fits b/src/blind_2pt_cosmosis/cosmosis_files/sim_fiducial_BLINDED.fits new file mode 100644 index 0000000..5df8eb5 Binary files /dev/null and b/src/blind_2pt_cosmosis/cosmosis_files/sim_fiducial_BLINDED.fits differ diff --git a/src/blind_2pt_cosmosis/cosmosis_files/sim_fiducial_BLINDED.fits_cov.png b/src/blind_2pt_cosmosis/cosmosis_files/sim_fiducial_BLINDED.fits_cov.png new file mode 100644 index 0000000..c053b7d Binary files /dev/null and b/src/blind_2pt_cosmosis/cosmosis_files/sim_fiducial_BLINDED.fits_cov.png differ diff --git a/src/blind_2pt_cosmosis/cosmosis_files/sim_fiducial_BLINDED.fits_gammat.png b/src/blind_2pt_cosmosis/cosmosis_files/sim_fiducial_BLINDED.fits_gammat.png new file mode 100644 index 0000000..ecfb45e Binary files /dev/null and b/src/blind_2pt_cosmosis/cosmosis_files/sim_fiducial_BLINDED.fits_gammat.png differ diff --git a/src/blind_2pt_cosmosis/cosmosis_files/sim_fiducial_BLINDED.fits_nz_lens.png b/src/blind_2pt_cosmosis/cosmosis_files/sim_fiducial_BLINDED.fits_nz_lens.png new file mode 100644 index 0000000..6491040 Binary files /dev/null and b/src/blind_2pt_cosmosis/cosmosis_files/sim_fiducial_BLINDED.fits_nz_lens.png differ diff --git a/src/blind_2pt_cosmosis/cosmosis_files/sim_fiducial_BLINDED.fits_nz_source.png b/src/blind_2pt_cosmosis/cosmosis_files/sim_fiducial_BLINDED.fits_nz_source.png new file mode 100644 index 0000000..9c776e0 Binary files /dev/null and b/src/blind_2pt_cosmosis/cosmosis_files/sim_fiducial_BLINDED.fits_nz_source.png differ diff --git a/src/blind_2pt_cosmosis/cosmosis_files/sim_fiducial_BLINDED.fits_wtheta.png b/src/blind_2pt_cosmosis/cosmosis_files/sim_fiducial_BLINDED.fits_wtheta.png new file mode 100644 index 0000000..951a4b8 Binary files /dev/null and b/src/blind_2pt_cosmosis/cosmosis_files/sim_fiducial_BLINDED.fits_wtheta.png differ diff --git a/src/blind_2pt_cosmosis/cosmosis_files/sim_fiducial_BLINDED.fits_xim.png b/src/blind_2pt_cosmosis/cosmosis_files/sim_fiducial_BLINDED.fits_xim.png new file mode 100644 index 0000000..32578f5 Binary files /dev/null and b/src/blind_2pt_cosmosis/cosmosis_files/sim_fiducial_BLINDED.fits_xim.png differ diff --git a/src/blind_2pt_cosmosis/cosmosis_files/sim_fiducial_BLINDED.fits_xip.png b/src/blind_2pt_cosmosis/cosmosis_files/sim_fiducial_BLINDED.fits_xip.png new file mode 100644 index 0000000..f6bf580 Binary files /dev/null and b/src/blind_2pt_cosmosis/cosmosis_files/sim_fiducial_BLINDED.fits_xip.png differ diff --git a/src/blind_2pt_cosmosis/cosmosis_files/type_table.txt b/src/blind_2pt_cosmosis/cosmosis_files/type_table.txt new file mode 100644 index 0000000..505ca4e --- /dev/null +++ b/src/blind_2pt_cosmosis/cosmosis_files/type_table.txt @@ -0,0 +1,16 @@ +#Type_1 Type_2 Section x y +galaxy_position_fourier galaxy_position_fourier galaxy_cl ell bin_{0}_{1} +galaxy_shear_emode_fourier galaxy_shear_emode_fourier shear_cl ell bin_{0}_{1} +galaxy_position_fourier galaxy_shear_emode_fourier galaxy_shear_cl ell bin_{0}_{1} +galaxy_shear_emode_fourier cmb_kappa_emode_fourier shear_cmbkappa_cl ell bin_{0}_{1} +cmb_shear_emode_fourier galaxy_shear_emode_fourier shear_cmbkappa_cl ell bin_{1}_{0} +galaxy_position_fourier cmb_kappa_emode_fourier galaxy_cmbkappa_cl ell bin_{0}_{1} +cmb_shear_emode_fourier galaxy_position_fourier galaxy_cmbkappa_cl ell bin_{1}_{0} +galaxy_shear_plus_real galaxy_shear_plus_real shear_xi_plus theta bin_{0}_{1} +galaxy_shear_plus_real galaxy_shear_plus_real shear_xi_plus theta bin_{0}_{1} +galaxy_shear_minus_real galaxy_shear_minus_real shear_xi_minus theta bin_{0}_{1} +galaxy_position_real galaxy_position_real galaxy_xi theta bin_{0}_{1} +galaxy_position_real galaxy_shear_plus_real galaxy_shear_xi theta bin_{0}_{1} +galaxy_position_real cmb_kappa_real galaxy_cmbkappa_xi theta bin_{0}_{1} +cmb_kappa_real cmb_kappa_real cmbkappa_xi theta bin_{0}_{1} +galaxy_shear_plus_real cmb_kappa_real shear_cmbkappa_xi theta bin_{0}_{1} \ No newline at end of file diff --git a/src/blind_2pt_cosmosis/cosmosis_files/values.ini b/src/blind_2pt_cosmosis/cosmosis_files/values.ini new file mode 100644 index 0000000..a47490d --- /dev/null +++ b/src/blind_2pt_cosmosis/cosmosis_files/values.ini @@ -0,0 +1,70 @@ +[cosmological_parameters] +omega_m = 0.1 0.3 0.9 +h0 = 0.55 0.69 0.91 +omega_b = 0.03 0.048 0.12 +n_s = 0.87 0.97 1.07 +sigma_8_input = 0.2 0.843 1.6 +A_s = 0.5e-09 2.19e-9 5.0e-09 +; omnuh2 = 0.0006 0.00083 0.01 +omnuh2 = 0.0 +massive_nu = 0 +massless_nu = 3.046 +omega_k = 0.0 +tau = 0.0697186 +;Helium mass fraction. Needed for Planck +yhe = 0.245341 + +;%include ../y3-3x2pt-methods/cosmosis/de_models/${DEMODEL}.ini +w = -1.0 +wa = 0.0 + +[shear_calibration_parameters] +m1 = -0.1 0.0 0.1 +m2 = -0.1 0.0 0.1 +m3 = -0.1 0.0 0.1 +m4 = -0.1 0.0 0.1 + +[wl_photoz_errors] +bias_1 = -0.1 0.0 0.1 +bias_2 = -0.1 0.0 0.1 +bias_3 = -0.1 0.0 0.1 +bias_4 = -0.1 0.0 0.1 + +[lens_photoz_errors] +bias_1 = -0.05 0.0 0.05 +bias_2 = -0.05 0.0 0.05 +bias_3 = -0.05 0.0 0.05 +bias_4 = -0.05 0.0 0.05 +bias_5 = -0.05 0.0 0.05 + +[bias_lens] +b1 = 0.8 1.7 3.0 +b2 = 0.8 1.7 3.0 +b3 = 0.8 1.7 3.0 +b4 = 0.8 2.0 3.0 +b5 = 0.8 2.0 3.0 + +[mag_alpha_lens] +; alpha_1 = -4. 0.903125 4. +; alpha_2 = -4. 0.68572965 4. +; alpha_3 = -4. 0.65340057 4. +; alpha_4 = -4. 1.58867862 4. +; alpha_5 = -4. 1.93754879 4. +alpha_1 = 0.903125 +alpha_2 = 0.68572965 +alpha_3 = 0.65340057 +alpha_4 = 1.58867862 +alpha_5 = 1.93754879 + +[intrinsic_alignment_parameters] +z_piv = 0.62 +A1 = -5.0 0.7 5.0 +A2 = -5.0 -1.36 5.0 +alpha1=-5.0 -1.7 5.0 +alpha2=-5.0 -2.5 5.0 +bias_ta = 0.0 1.0 2.0 +#bias_tt = 1.0 + + + + diff --git a/src/blind_2pt_cosmosis/io.py b/src/blind_2pt_cosmosis/io.py new file mode 100644 index 0000000..1cd8d8d --- /dev/null +++ b/src/blind_2pt_cosmosis/io.py @@ -0,0 +1,149 @@ +import argparse +from . import __version__ + +DEFAULT_PARAM_RANGE = {'cosmological_parameters--sigma_8_input':(0.834-3*.04,0.834+3*0.04),\ + 'cosmological_parameters--w':(-1.5,-.5)} + +# DEFAULT_PARAM_RANGE = {'cosmological_parameters--sigma_8_input':(-3*0.04, 3*0.04),\ +# 'cosmological_parameters--w':(-0.5,0.5)} + + +class DictAction(argparse.Action): + def __call__(self, parser, namespace, values, option_string=None): + setattr(namespace, self.dest, eval(values)) + +def get_stored_seed_and_tag(args): + """ + Obtains the seed and tag to be used in the output filename + + Parameters + ---------- + args : argparse.Namespace + Namespace with the arguments given by the user + + Returns + ------- + storeseed : str + Seed to be stored in the fits file + tagstr : str + Tag to be used in the output filename + """ + if args.outftag is not None: + tagstr = args.outftag # +'_'+tagstr + if args.seedinfname: + tagstr = tagstr + '_'+ args.seedstring + elif args.seedinfname: + tagstr = '_' + args.seedstring + else: + tagstr = '' + + if args.seedinfits: + storeseed = args.seedstring + else: + storeseed = 'notsaved' + return storeseed, tagstr + +def get_parser(): + """ + creates the parser to obtain the user command line input + """ + parser = argparse.ArgumentParser( + prog='blind_2pt_cosmosis', + formatter_class=argparse.RawTextHelpFormatter, + description=f''' + -------------------------------------------------------------------------------- + Blinding Module for 2pt data in Cosmosis. + Version {__version__} + + This is an adaptation of Jessie Muir's blind_2pt_usingcosmosis.py scripts. + -------------------------------------------------------------------------------- + This module will apply a blinding factor to 2pt functios stored in a fits + file, using cosmosis to compute the blinding factors. + + The workhorse function that puts everything together is do2ptblinding(). + This script can be called with command line arguments, or by calling that + function in another script. + + The script uses cosmosis, so you'll need to source its setup file before + running this. + + ............................................................ + What it does: + 1. Using a string seed, pseudo-randomly draw a shift in parameters. This + will eitherbe drawn from a flat distribution in some predefined parameter + ranges, or from a distribution defined in paramshiftmodule. Return + dictionary of shifts where keys are parameter names matching those expected + by run_cosmosis_togen_2ptdict. + -> See draw_paramshift + + 2. Using a cosmosis parameter ini file template, compute 2pt functions at + reference and shifted cosmologies by running the cosmosis twice. This + should work no matter what sampler shows up in the template ini file. + (If it uses the test sampler, make sure the output is nothing, so that + it doesn't save the cosmology and 2pt info.) Gets the 2pt functions into + dictionaries of a specific format. + -> See run_cosmosis_togen_2ptdict + + 4. Take ratio or difference of 2pt datavectors to get blinding factor, in same + dictionary format. (NOTE, MULTIPLICATIVE BLINDING IS CURRENTLY DISABLED) + -> See get_factordict + + 5. Apply blinding factor to desired datafile, saving a new fits file with + string key in filename. + -> See apply2ptblinding_tofits + ............................................................ + Script maintained by Jessica Muir (jlmuir@stanford.edu). + Module maintained by Arthur Loureiro (arthur.loureiro@fysik.su.se) + -------------------------------------------------------------------------------- + ''', fromfile_prefix_chars='@') + + + parser.add_argument("-u", "--origfits", type=str, required=True, + help='Name of unblinded fits file') + + parser.add_argument("-i", "--ini", type=str, required=False, + default='./src/blind_2pt_cosmosis/cosmosis_files/default_blinding_template.ini', + help='Ini file containing template for generating 2pt functions. \nShould'+ + 'use a binning that matches that of the origfits file' + + '\nReference a values file centered at desired reference' + + '\ncosmology') + + parser.add_argument('-s', '--seed', type=str, required=False, + default="HARD_CODED_BLINDING", + help='string used to seed parameter shift selection') + + parser.add_argument('-t', '--bftype', type=str, required=False, + default='add', + help="Blinding factor type. Can be 'add', 'mult', or" + + "'multNOCS' (mult with no cov scaling) \n[MULT OPTION IS DISABLED]. \n>> Default is 'add'") + + parser.add_argument('-o', '--outfname', type=str, required=False, + default=None, + help="Output filename; only set for testing, " + + "If outftag is set, then " + + "\noutfname = [origfits]_[outftag]_[seed].fits." + + "\n>> Default behavior is outfname = [origfits]_[seed].fits.") + + parser.add_argument('-f', '--outftag', type=str, required=False, + default="_BLINDED", + help="String to label output filename, for testing purposes. \nIf set" + + "when outfname is None (default), \nthen " + + "outfname = [origfits]_[outftag]_[seed].fits") + + parser.add_argument('-p', '--paramshifts', action=DictAction, required=False, + default=DEFAULT_PARAM_RANGE, + help="Dictionary of parameter shifts between quotes \". \nPlease use the parameter names" + + "as named in Cosmosis. \n>> Default is \"{'cosmological_parameters--sigma8_input':(0.834-3*.04,0.834"+ + "+3*0.04),\ncosmological_parameters--w':(-1.5,-.5)}\"") + + parser.add_argument('--seedinfname', action='store_true', required=False, + default=False, + help="If set, appends seed to blinded data filename. Default is False.") + + parser.add_argument('--seedinfits', action='store_true', required=False, + default=False, + help="If set, stores seed in KEYWORD entry in blinded fits file. Default True.") + + parser.add_argument("--log-level", choices=["DEBUG", "INFO", "WARNING", "ERROR", "CRITICAL"], + default="INFO", help="Specify the logging level: DEBUG, INFO, WARNING, ERROR, or CRITICAL. Default is INFO.") + return parser diff --git a/src/blind_2pt_cosmosis/param_shifts.py b/src/blind_2pt_cosmosis/param_shifts.py new file mode 100644 index 0000000..d3336cf --- /dev/null +++ b/src/blind_2pt_cosmosis/param_shifts.py @@ -0,0 +1,121 @@ +import hashlib +import logging +import numpy as np +from .io import DEFAULT_PARAM_RANGE + +logger = logging.getLogger("2pt_blinding") + +def draw_flat_param_shift(seedstring='blinded', ranges=None): + """ + Given a string seed, pseudo-randomly draws shift in parameter space. + + By default, draws random values for all the parameters with ranges defined + in the dictionary 'ranges'. Note: + - parameter names should match those used by cosmosis (see + DEFAULT_PARAM_RANGE dict for an example) + - ranges is a dictionary with parameter names as keys, the values are + tuples set up as (min,max). The parameter will be drawn from a flat + distribution between min and max. + + Make sure any parameters that are shifted here have names matching + how cosmosis uses them, so that the code in the for loop starting with + 'for parameter in pipeline.parameters' in run_cosmosis_togen_2ptdict + will work. + """ + # sets the seed: + seedind = int(int(hashlib.md5(seedstring.encode('utf-8')).hexdigest(), 16) % 1.e8) + np.random.seed(seedind) + + if ranges is None: + ranges = DEFAULT_PARAM_RANGE + + # sorting makes sure it's always the same order + params2shift = sorted(ranges.keys()) + Nparam = len(params2shift) + # arrays between 0 and 1 + shiftfrac = np.random.rand(Nparam) + mins, maxs = zip(*[ranges[k] for k in params2shift]) + dparams = np.array(mins) + (np.array(maxs) - np.array(mins)) * shiftfrac + pdict = {param: value for param, value in zip(params2shift, dparams)} + pdict['SHIFTS'] = False + + return pdict + +def apply_parameter_shifts(pipeline, pdict): + """ + Apply parameter shifts to the pipeline parameters. + """ + doshifts = len(list(pdict.keys())) + if doshifts > 0: + SHIFTS = pdict['SHIFTS'] + haveshifted = [] + # set parameters as desired + for parameter in pipeline.parameters: + key = str(parameter) + #print(">>>>>>>>>>>>>>>>>>>>>>>>>>>>>......", key) + #set the values + if doshifts > 0: + try: + if SHIFTS: + logger.debug(f"Shifting parameter {key} by {pdict[key]}") + parameter.start = parameter.start + pdict[key] + logger.debug(f">>>>>>>>>>>>>>>>>>>>>>>>>>>>>...... {parameter.start}") + else: + logger.debug(f"Shifting parameter {key} to {pdict[key]}") + parameter.start = pdict[key] + logger.debug(f">>>>>>>>>>>>>>>>>>>>>>>>>>>>>...... {parameter.start}") + haveshifted.append(key) + except KeyError: + print(" asked for shifts in:", key) + logger.debug(f"Parameter {key} not in pdict. Not shifting.") + pass + + if (doshifts > 0) and (len(haveshifted) != (len(list(pdict.keys())) -1)): + print(" asked for shifts in:",list(pdict.keys())) + print(" did shifts in:",haveshifted) + raise ValueError("WARNING: YOU ASKED FOR SHIFTS IN PARAMTERS NOT IN THE COSMOSIS PIPELINE.") + return pipeline + +def get_factordict(refdict,shiftdict,bftype='add'): + """ + Given two point dictionaries for reference and shifted cosmology, + returns dictionary of blinding factors. + + Parameters + ---------- + refdict : dict + Dictionary of reference cosmology 2pt functions. + shiftdict : dict + Dictionary of shifted cosmology 2pt functions. + bftype : str, optional + What kind of blinding factor is it? + 'add' : bf = - ref + shift + 'mult': bf = shift/ref + 'multNOCS': bf = shift/ref, but with no cosmic shear + (i.e. no E-mode or B-mode shear) + Default is 'add'. + + Returns + ------- + factordict : dict + Dictionary of blinding factors. + """ + + factordict = {} + for key in refdict: + end = key[key.rfind('_')+1:] + #print(key,end) + # don't take ratios or differences of angle/multipole info + if end in ['ell','l','theta','bins','angbins','binavg','mins','maxs']: + #print ' no change' + factordict[key] = refdict[key] + else: + if bftype=='mult' or bftype=='multNOCS': + logger.debug(' dividing!') + factordict[key] = shiftdict[key]/refdict[key] + elif bftype=='add': + logger.debug(' adding!') + factordict[key] = shiftdict[key] - refdict[key] + else: + raise ValueError('In get_factordict: blinding factor type not recognized') + return factordict \ No newline at end of file diff --git a/src/blind_2pt_cosmosis/run_cosmosis_2pt.py b/src/blind_2pt_cosmosis/run_cosmosis_2pt.py new file mode 100644 index 0000000..7986f93 --- /dev/null +++ b/src/blind_2pt_cosmosis/run_cosmosis_2pt.py @@ -0,0 +1,89 @@ +import time +import logging +from cosmosis.runtime.config import Inifile +from cosmosis.runtime.pipeline import LikelihoodPipeline +from .twopt_utils import get_twoptdict_from_pipeline_data +from .param_shifts import apply_parameter_shifts + +logger = logging.getLogger("2pt_blinding") + +def setup_pipeline(inifile, angles_file=None, nz_file=None): + """ + Initialize and set up the Cosmosis pipeline. + """ + ini = Inifile(inifile) + # Modify settings to suppress unnecessary outputs + ini = modify_settings(ini, angles_file, nz_file) + pipeline = LikelihoodPipeline(ini) + return pipeline + +def modify_settings(ini, angles_file=None, nz_file=None): + """ + Modify pipeline settings to suppress unnecessary outputs. + """ + sections_to_modify = ['test', 'output'] + for section in sections_to_modify: + if section in ini.__dict__['_sections']: + ini.__dict__['_sections'][section]['save_dir'] = '' + ini.__dict__['_sections']['pipeline']['debug'] = 'F' + ini.__dict__['_sections']['pipeline']['quiet'] = 'T' + + # a structure like what follows could be used to make sure + # angles and n(z) are being used consistently with file being blinded + # as in, make sure that they are referencing it for angles and n(z). + # However, this isn't really ideal, as how it is set up will be specific + # to the 3x2pt pipeline and may cause issues when/if the same script + # is bein used for 5x2pt or other observables. Really the template ini + # file should handle these choices, and hsould be set up in the same way + # that you would set up an ini file to do parameter estimation on your + # data file. + if angles_file is not None: #fits file to get theta ranges from + hadsections=[] + for section in ['shear_2pt_eplusb','shear_2pt_eminusb','2pt_gal','2pt_gal_shear']: + if section in ini.__dict__['_sections'].keys(): + ini.__dict__['_sections'][section]['theta_file']=angles_file + hadsections.append(section) + print("CHANGED theta_file TO ",angles_file," FOR:",hadsections," \n IF YOU MEANT TO CHANGE IT FOR OTHER 2PT FUNCTIONS, SOMETHING WENT WRONG.") + + if nz_file is not None: # fits file to get n(z) from + section = 'fits_nz' + if section in ini.__dict__['_sections'].keys(): + ini.__dict__['_sections'][section]['nz_file']=nz_file + else: + raise ValueError("You specified nz_file as "+nz_file+", but I can't find the fits_nz module settings in your ini file.") + return ini + +def run_pipeline(pipeline): + """ + Run the Cosmosis pipeline and return the data. + """ + data = pipeline.run_parameters([]) + return data + +def run_cosmosis_togen_2ptdict(inifile, pdict={}, + nz_file=None, angles_file=None): + """ + Runs cosmosis pipeline to generate 2pt functions. + """ + logger.debug("Running cosmosis pipeline to generate 2pt functions.") + pipeline = setup_pipeline(inifile, angles_file, nz_file) + logger.debug("Passed setup_pipeline.") + # need to set all of the parameters to be fixed for run_parameters([]) + # to work. Doing this will effectively run things like the test sampler + # no matter what sampler is listed in the ini file + for parameter in pipeline.parameters: + pipeline.set_fixed(parameter.section, parameter.name, parameter.start) + # if we need to do shifts: + if len(list(pdict.keys())) > 0: + logger.debug("Applying parameter shifts.") + pipeline = apply_parameter_shifts(pipeline, pdict) + else: + logger.debug("No parameter shifts to apply.") + + data = run_pipeline(pipeline) + logger.debug("Passed run_pipeline.") + + twoptdict = get_twoptdict_from_pipeline_data(data) + logger.debug(f"Passed get_twoptdict_from_pipeline_data!") + + return twoptdict diff --git a/src/blind_2pt_cosmosis/twopt_utils.py b/src/blind_2pt_cosmosis/twopt_utils.py new file mode 100644 index 0000000..6100637 --- /dev/null +++ b/src/blind_2pt_cosmosis/twopt_utils.py @@ -0,0 +1,400 @@ +import numpy as np +import os +import logging +from astropy.io import fits +import shutil +from scipy.interpolate import interp1d +from astropy.table import Table + +logger = logging.getLogger("2pt_blinding") + +# # gambiarra to use the type_table from cosmosis-standard-library +def load_type_table(): + # Taken from: + # https://github.com/joezuntz/cosmosis-standard-library/ + # blob/sacc/likelihood/2pt/twopoint_cosmosis.py + # with Joe Zuntz consent + dirname = os.path.split(__file__)[0] + table_name = os.path.join(dirname, "cosmosis_files/type_table.txt") + type_table = Table.read(table_name, format="ascii.commented_header") + table = {} + for (type1, type2, section, x, y) in type_table: + table[(type1, type2)] = (section, x, y) + return table + +type_table = load_type_table() + +class SpectrumInterp(object): + """ + This is copied from 2pt_like, for get_data_from_dict_for_2pttype + + Should not get used if bin averaging is being used; is in place + for if theory vector in datablock is very densely sampled, this + gets used to pick out values corresponding to desired angle positions + """ + def __init__(self,angle,spec,bounds_error=True): + if np.all(spec>0): + self.interp_func=interp1d(np.log(angle),np.log(spec),bounds_error=bounds_error,fill_value=-np.inf) + self.interp_type='loglog' + elif np.all(spec<0): + self.interp_func=interp1d(np.log(angle),np.log(-spec),bounds_error=bounds_error,fill_value=-np.inf) + self.interp_type='minus_loglog' + else: + self.interp_func=interp1d(np.log(angle),spec,bounds_error=bounds_error,fill_value=0.) + self.interp_type="log_ang" + + def __call__(self,angle): + if self.interp_type=='loglog': + spec=np.exp(self.interp_func(np.log(angle))) + elif self.interp_type=='minus_loglog': + spec=-np.exp(self.interp_func(np.log(angle))) + else: + assert self.interp_type=="log_ang" + spec=self.interp_func(np.log(angle)) + return spec + +def spectrum_array_from_block(block, section_name, types, xlabel='theta', bin_format = 'bin_{0}_{1}'): + """ + >> Updated 4/21/20 to work with bin averaging implemented for DES Y3 + + Initially adapting this from save_2pt.py and 2pt_like code. + + No scale cutting implemented or angular sampling; + that will be handled later. Since we + track bin numbers and angle values, which will be checked against + the unblinded fits file in get_dictdat_tomatch_fitsdat. + """ + if xlabel=='theta': + is_binavg = block[section_name,"bin_avg"] + else: + is_binavg = False #no bin averaging in fourier space + + # for cross correlations we must save bin_ji as well as bin_ij. + # but not for auto-correlations. + is_auto = (types[0] == types[1]) + if block.has_value(section_name, "nbin"): + nbin_a = block[section_name, "nbin"] + else: + nbin_a = block[section_name, "nbin_a"] + nbin_b = block[section_name, "nbin_b"] + + + #This is the ell/theta values that have been calculated by cosmosis, + # if bin averaging, these should match what's in the fits files (up to rounding) + # if interpolating, angles will be more densely sampled than values in fits files. + theory_angle = block[section_name,xlabel] + n_angle = len(theory_angle) #whatevers in the block, length of array + if is_binavg: + theory_angle_edges = block[section_name,xlabel+'_edges'] #angle bin edges + + #The fits format stores all the measurements + #as one long vector. So we build that up here from the various + #bins that we will load in. These are the different columns + value = [] + bin1 = [] + bin2 = [] + angles = [] + # these will only be used if averaging over angular bins + angle_mins = [] + angle_maxs = [] + + # n.b. don't have suffix option implemented here + + #Bin pairs. Varies depending on auto-correlation + for i in range(nbin_a): + if is_auto: + jmax = i+1 + else: + jmax = nbin_b + for j in range(jmax): + #Load and interpolate from the block + binlabel = bin_format.format(i+1,j+1) + if block.has_value(section_name, binlabel): + cl = block[section_name, binlabel] + bin1.append(np.repeat(i + 1, n_angle)) + bin2.append(np.repeat(j + 1, n_angle)) + angles.append(theory_angle) + if is_binavg: + angle_mins.append(theory_angle_edges[:-1]) + angle_maxs.append(theory_angle_edges[1:]) + value.append(cl) + + if is_auto and i!=j: #also store under flipped z bin labels + # this allows the script to work w fits files uing either convention + bin1.append(np.repeat(j + 1, n_angle)) + bin2.append(np.repeat(i + 1, n_angle)) + value.append(cl) + angles.append(theory_angle) + if is_binavg: + angle_mins.append(theory_angle_edges[:-1]) + angle_maxs.append(theory_angle_edges[1:]) + + + #Convert all the lists of vectors into long single vectors + value = np.concatenate(value) + bin1 = np.concatenate(bin1) + bin2 = np.concatenate(bin2) + bins = (bin1,bin2) + angles = np.concatenate(angles) + if is_binavg: + angle_mins = np.concatenate(angle_mins) + angle_maxs = np.concatenate(angle_maxs) + else: + angle_mins = None + angle_maxs = None + + return angles, value, bins, is_binavg, angle_mins, angle_maxs + +def get_dictkey_for_2pttype(type1, type2): + """ + Convert strings used in fits file to label spectra type in fits file to + dictionary keys expected by this script's functions. + (fits file 2pt table naming -> this script's naming) + """ + logger.debug("Requested: \ttype1: {0:s}, type2: {1:s}".format(type1, type2)) + mapping = { + 'G': 'galaxy', + 'C': 'cmb', + 'P': 'position', + 'E': 'shear_emode', + 'B': 'shear_bmode', + '+': 'shear_plus', + '-': 'shear_minus', + 'K': 'kappa', + 'R': 'real', + 'F': 'fourier' + } + + if type1 in ("GPF", "GEF", "GBF", "GPR", "G+R", "G-R", "CKR", "GPF", "GEF", "GBF"): + # Translate short codes into longer strings + newtypes = ['_'.join([mapping[t[0]], mapping[t[1]], mapping[t[2]]]) for t in [type1, type2]] + type1 = newtypes[0] + type2 = newtypes[1] + + try: + section, xlabel, ylabel = type_table[(type1, type2)] + ykey = section + xkey = f"{section}_{xlabel}" + except KeyError: + raise ValueError("Spectra type not recognized: {0:s}, {1:s} ".format(type1, type2)) + logger.debug("Returning: \txkey: {0:s}, ykey: {1:s}".format(xkey, ykey)) + return xkey, ykey + +def get_twoptdict_from_pipeline_data(data): + """ + Extract 2pt data from a cosmosis pipeline data object and return a dictionary + with keys corresponding to the 2pt spectra types. + + For auto spectra where z bin label order doesn't matter, both orders + will be stored. (e.g. the same data will be there for bins 1-2 and 2-1). + """ + outdict = {} + type_keys = [types for types in type_table if data.has_section(type_table[types][0])] + + for types in type_keys: + section, xlabel, binformat = type_table[types] + xkey, ykey = get_dictkey_for_2pttype(types[0], types[1]) + + x, y, bins, is_binavg, x_mins, x_maxs = spectrum_array_from_block(data, section, types, xlabel, binformat) + + outdict[xkey] = x + outdict[ykey] = y + outdict[ykey + '_bins'] = bins + outdict[ykey + '_binavg'] = is_binavg + + if is_binavg: + outdict[xkey + '_mins'] = x_mins + outdict[xkey + '_maxs'] = x_maxs + + return outdict + +def apply_2pt_blinding_and_save_fits(factordict, origfitsfile, outfname=None, outftag="_BLINDED", + justfname=False,bftype='add', storeseed='notsaved'): + """ + Given the dictionary of one set of blinding factors, + the name of a of a fits file containing unblinded 2pt data, + and (optional) desired output file name or tag, + multiplies 2pt data in original fits file by blinding factors and + saves results (blinded data) into a new fits file. + + If argument is passed for outfname, that will be the name of output file, + if not, will be _.fits. The script will check + that the outfname is different than the original, unblinded file. If they + match, it will revert to default behavior for naming the blinded file + so that the unblinded file isn't overwritten. + > if storeseed=True, will save whatever string is there to entry in fits file + + If justfname == True, doesn't do any file manipulation, just returns + string of output filename. (for testing) + + bftype can be 'add','mult' or 'mult-nocs'. For data vec d, blinding factor f + 'add' - do additive blinding: + d_blind = d_input + f, f = d_shift - d_ref, cov_bl = cov + 'mult' - do multiplicative blinding, scale covmat in blinded file: + d_blind = d_input*f, f = d_shift/d_ref, cov_bl = f^T*cov*f + 'multNOCS' - do multiplicative blinding, but without covariance scaling + d_blind = d_input*f, f = d_shift/d_ref, cov_bl = cov + """ + logger.info(f'apply2ptblinding for: {origfitsfile}') + logger.info(f'bftype: {bftype}') + # check whether data is already blinded and whether Nbins match + for table in fits.open(origfitsfile): #look through tables to find 2ptdata + if table.header.get('2PTDATA'): + if table.header.get('BLINDED'): #check for blinding + #if entry not there, or storing False -> not already blinded + raise ValueError('Data is already blinded!') + return + + # set up output file + if outfname == None or outfname==origfitsfile: + #make sure you can't accidentaly overwrite the original + #if output filename isn't given, add outftag onto the name of the + # unblinded file + if outftag == None or outftag =='': + outftag = 'BLINDED-{0:s}-defaulttag'.format(bftype) + outfname = origfitsfile.replace('.fits','{0}.fits'.format(outftag)) + + + if not justfname: + shutil.copyfile(origfitsfile,outfname) + + hdulist = fits.open(outfname, mode='update') #update lets us write over + + #apply blinding factors + for table in hdulist: #look all tables + if table.header.get('2PTDATA'): + factor = get_dictdat_tomatch_fitsdat(table, factordict) + + if bftype=='mult' or bftype=='multNOCS': + #print 'multiplying!' + table.data['value'] *= factor + if bftype=='mult': #store info about how to scale covmat + type1 = table.header['QUANT1'] + type2 = table.header['QUANT2'] + raise ValueError('Not currently set up to do covariance scaling needed for multiplicative blinding.') + + elif bftype=='add': + #print 'adding!' + table.data['value'] += factor + else: + raise ValueError('bftype {0:s} not recognized'.format(bftype)) + + #add new header entry to note that blinding has occurred, and what type + table.header['BLINDED'] = bftype + table.header['KEYWORD'] = storeseed + + hdulist.close() # will save new data to file if 'update' was passed when opened + logger.info(f">>>> Stored blinded data in {outfname}") + return outfname + +def get_dictdat_tomatch_fitsdat(table, dictdata): + """ + Given table of type fits.hdu.table.BinTableHDU containing 2pt data, + retrieves corresponding data from dictionary (blinding factors). + + Expects that same z and theta bin numbers correspond to the same + z and theta values (i.e. matches up bin numbers but doesn't do + any interpolation). Theta values will be checked, but z values won't. + """ + if not table.header.get('2PTDATA'): + logger.warning("Can't match dict data: this fits table doesn't contain 2pt data. Is named:",table.name) + return + type1 = table.header['QUANT1'] + type2 = table.header['QUANT2'] + + bin1 = table.data['BIN1'] #which bin is quant1 from? + bin2 = table.data['BIN2'] #which bin is quant2 from? + + # check for bin averaging + if "ANGLEMIN" in table.data.names: + fits_is_binavg = True + xfromfits_mins = table.data['ANGLEMIN'] + xfromfits_maxs = table.data['ANGLEMAX'] + else: + fits_is_binavg = False + xfromfits_mins = None + xfromfits_maxs = None + + xfromfits = table.data['ANG'] + + yfromdict = get_data_from_dict_for_2pttype(type1, type2, bin1, bin2, xfromfits, + dictdata, fits_is_binavg=fits_is_binavg, + xfits_mins=xfromfits_mins, xfits_maxs=xfromfits_maxs) + + return yfromdict + +def get_data_from_dict_for_2pttype(type1, type2, bin1fits, bin2fits, xfits, datadict, + fits_is_binavg=True, xfits_mins=None, xfits_maxs=None): + """ + Given info about 2pt data in a fits file (spectra type, z bin numbers, + and angular bin numbers, extracts 2pt data from a dictionary of + e.g. blinding factors to match, with same array format and bin ordering + as the fits file data. + """ + xkey,ykey = get_dictkey_for_2pttype(type1,type2) + + is_binavg = datadict[ykey+'_binavg'] + # this check is probably unnecessary, since the data is always bin averaged + # if dict_is_binavg != fits_is_binavg: + # raise ValueError("Theory calc and fits file aren't consistent in whether they do bin averaging vs interpolation for 2pt calculations.") + if fits_is_binavg and ((xfits_mins is None) or (xfits_maxs is None)): + raise ValueError("Fits file is bin-averaged but I couldn't find theta values for bin edges.") + + if 'theta' in xkey: #if realspace, put angle data into arcmin + xmult = 60.*180./np.pi # change to arcmin + else: + xmult = 1. #fourier space + + xfromdict = datadict[xkey]*xmult #in arcmin (is in radians in datablock) + xfromdict = xfromdict*xmult + if is_binavg: + xfromdict_mins = datadict[xkey+'_mins']*xmult + xfromdict_maxs = datadict[xkey+'_maxs']*xmult + yfromdict = datadict[ykey] + binsdict = datadict[ykey+'_bins'] + b1dict = binsdict[0] + b2dict = binsdict[1] + + #if get theory calcs in same format as fits ones + Nentries = bin1fits.size + yout = np.zeros(Nentries) + #this structure will store interp functions for b1-b2 combos + # so we don't have to keep recreating them + if not is_binavg: + Nb1fits = max(bin1fits) + Nb2fits = max(bin2fits) + interpfuncs = [[None for b2f in range(Nb2fits)]\ + for b1f in range(Nb1fits)] + + for i in range(Nentries): + b1 = bin1fits[i] + b2 = bin2fits[i] + if is_binavg: + wherebinsmatch = (b1==b1dict)*(b2==b2dict) + + roundto=4 #round for matching, since there are more decimals in fits than datablock + wherexmatch_mins = np.around(xfits_mins[i],roundto)==np.around(xfromdict_mins,roundto) + wherexmatch_maxs = np.around(xfits_maxs[i],roundto)==np.around(xfromdict_maxs,roundto) + wherexmatch = wherexmatch_mins*wherexmatch_maxs + whichinds = wherebinsmatch*wherexmatch + howmany = np.sum(whichinds) + if howmany==0: #no matches + raise ValueError("No theory calc match for data point: {0:s}, z bins = ({1:d},{2:d}), theta between [{3:0.2f},{4:0.2f}]".format(ykey,b1,b2,xfits_mins[i],xfits_maxs[i])) + elif howmany>1: #duplicate matches + raise ValueError("More than one theory calc match for data point: {0:s}, z bins = ({1:d},{2:d}), theta between [{3:0.2f},{4:0.2f}]".format(ykey,b1,b2,xfits_mins[i],xfits_maxs[i])) + + yout[i] = yfromdict[whichinds] + else: + whichinds = (b1==b1dict)*(b2==b2dict)#*(ab==angbdict) + # get x and y info for that bin combo + tempx = xfromdict[whichinds] + tempy = yfromdict[whichinds] + if interpfuncs[b1-1][b2-1] is None: #no interpfunc yet, set it up + yinterp = SpectrumInterp(tempx,tempy) + interpfuncs[b1-1][b2-1] = yinterp + else: + yinterp = interpfuncs[b1-1][b2-1] + yout[i] = yinterp(xfits[i]) + # We're returning the y data from the dictionary's array + # interpolated to match the theta values in the fits file. + return yout \ No newline at end of file diff --git a/tests/__init__.py b/tests/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/tests/test_io.py b/tests/test_io.py new file mode 100644 index 0000000..7f9a7e0 --- /dev/null +++ b/tests/test_io.py @@ -0,0 +1,51 @@ +import pytest +import argparse + +from blind_2pt_cosmosis.io import get_parser, DictAction, get_stored_seed_and_tag + +def test_dict_action(): + parser = argparse.ArgumentParser() + parser.add_argument('--test', action=DictAction) + args = parser.parse_args(['--test', "{'key': 'value'}"]) + assert args.test == {'key': 'value'} + +def test_get_parser(): + parser = get_parser() + assert isinstance(parser, argparse.ArgumentParser) + +def test_get_stored_seed_and_tag(): + # Test case 1: Default values + args = argparse.Namespace(outftag=None, seedinfname=False, seedinfits=False, seedstring='myseed') + storeseed, tagstr = get_stored_seed_and_tag(args) + assert storeseed == 'notsaved' + assert tagstr == '' + + # Test case 2: outftag provided + args = argparse.Namespace(outftag='output', seedinfname=False, seedinfits=False, seedstring='myseed') + storeseed, tagstr = get_stored_seed_and_tag(args) + assert storeseed == 'notsaved' + assert tagstr == 'output' + + # Test case 3: seedinfname and outftag provided + args = argparse.Namespace(outftag='output', seedinfname=True, seedinfits=False, seedstring='myseed') + storeseed, tagstr = get_stored_seed_and_tag(args) + assert storeseed == 'notsaved' + assert tagstr == 'output_myseed' + + # Test case 4: seedinfits provided + args = argparse.Namespace(outftag=None, seedinfname=False, seedinfits=True, seedstring='myseed') + storeseed, tagstr = get_stored_seed_and_tag(args) + assert storeseed == 'myseed' + assert tagstr == '' + + # Test case 5: seedinfname and seedinfits provided + args = argparse.Namespace(outftag=None, seedinfname=True, seedinfits=True, seedstring='myseed') + storeseed, tagstr = get_stored_seed_and_tag(args) + assert storeseed == 'myseed' + assert tagstr == '_myseed' + + # Test case 6: All options provided + args = argparse.Namespace(outftag='output', seedinfname=True, seedinfits=True, seedstring='myseed') + storeseed, tagstr = get_stored_seed_and_tag(args) + assert storeseed == 'myseed' + assert tagstr == 'output_myseed' \ No newline at end of file diff --git a/tests/test_param_shifts.py b/tests/test_param_shifts.py new file mode 100644 index 0000000..a6143d4 --- /dev/null +++ b/tests/test_param_shifts.py @@ -0,0 +1,96 @@ +import pytest +import numpy as np + +from blind_2pt_cosmosis.param_shifts import draw_flat_param_shift +from blind_2pt_cosmosis.param_shifts import DEFAULT_PARAM_RANGE +from blind_2pt_cosmosis.param_shifts import apply_parameter_shifts + +@pytest.fixture +def ranges(): + return { + 'param1': (0, 1), + 'param2': (-1, 1), + 'param3': (10, 20) + } + +def test_parameter_shifts(ranges): + pdict = draw_flat_param_shift(ranges=ranges) + assert isinstance(pdict, dict) + assert len(pdict) == len(ranges) + 1 # +1 for 'SHIFTS' key + assert not pdict['SHIFTS'] + + for param, value in pdict.items(): + if param != 'SHIFTS': + assert value >= ranges[param][0] + assert value <= ranges[param][1] + +def test_parameter_shifts_with_default_range(): + pdict = draw_flat_param_shift() + assert isinstance(pdict, dict) + assert len(pdict) == len(DEFAULT_PARAM_RANGE) + 1 # +1 for 'SHIFTS' key + assert not pdict['SHIFTS'] + + for param, value in pdict.items(): + if param != 'SHIFTS': + assert value >= DEFAULT_PARAM_RANGE[param][0] + assert value <= DEFAULT_PARAM_RANGE[param][1] + +def test_draw_flat_param_shift(): + # Set a known seed to make the random draw predictable + seedstring = 'test_seed' + np.random.seed(42) # Seed for np.random.rand + + # Test case 1: Using default parameter ranges + pdict = draw_flat_param_shift(seedstring) + assert 'SHIFTS' in pdict + assert pdict['SHIFTS'] is False + for param, (min_val, max_val) in DEFAULT_PARAM_RANGE.items(): + assert param in pdict + assert min_val <= pdict[param] <= max_val + + # Test case 2: Custom parameter ranges + custom_ranges = {'param1': (0.0, 1.0), 'param2': (-1.0, 1.0)} + pdict = draw_flat_param_shift(seedstring, custom_ranges) + assert 'SHIFTS' in pdict + assert pdict['SHIFTS'] is False + for param, (min_val, max_val) in custom_ranges.items(): + assert param in pdict + assert min_val <= pdict[param] <= max_val + + # Test case 3: Check that a different seed produces different results + pdict_1 = draw_flat_param_shift(seedstring) + pdict_2 = draw_flat_param_shift('different_seed') + assert pdict_1 != pdict_2 + + # Test case 4: Check that providing the same seed produces the same result + pdict_1 = draw_flat_param_shift(seedstring) + pdict_2 = draw_flat_param_shift(seedstring) + assert pdict_1 == pdict_2 + +# def test_apply_parameter_shifts(): +# # Create a mock pipeline with parameters +# class MockParameter: +# def __init__(self, name, start): +# self.name = name +# self.start = start + +# class MockPipeline: +# def __init__(self, parameters): +# self.parameters = parameters.name + +# # Define test parameters and expected results +# parameters = [MockParameter('param1', 1.0), MockParameter('param2', 2.0)] +# pdict = {'param1': 0.5, 'param2': 0.0, 'SHIFTS': True} +# expected_start_values = {'param1': 1.5, 'param2': 2.0} + +# # Apply parameter shifts +# pipeline = MockPipeline(parameters) +# modified_pipeline = apply_parameter_shifts(pipeline, pdict) + +# # Check that the parameter values were shifted as expected +# for parameter in modified_pipeline.parameters: +# assert parameter.name in expected_start_values +# assert parameter.start == expected_start_values[parameter.name] + +# # Ensure that the 'SHIFTS' parameter is removed from the dictionary +# assert 'SHIFTS' not in pdict diff --git a/tests/test_twopt_utils.py b/tests/test_twopt_utils.py new file mode 100644 index 0000000..875aee5 --- /dev/null +++ b/tests/test_twopt_utils.py @@ -0,0 +1,61 @@ +import pytest +from blind_2pt_cosmosis.twopt_utils import get_dictkey_for_2pttype +from blind_2pt_cosmosis.twopt_utils import get_twoptdict_from_pipeline_data +from blind_2pt_cosmosis.twopt_utils import spectrum_array_from_block + +# FIXME: MISSING TESTS FOR spectrum_array_from_block + +def test_get_dictkey_for_2pttype(): + + type1 = 'GPR' + type2 = 'G+R' + xkey, ykey = get_dictkey_for_2pttype(type1, type2) + + assert xkey == 'galaxy_shear_xi_theta' + assert ykey == 'galaxy_shear_xi' + + type1 = 'G+R' + type2 = 'CKR' + xkey, ykey = get_dictkey_for_2pttype(type1, type2) + + assert xkey == 'shear_cmbkappa_xi_theta' + assert ykey == 'shear_cmbkappa_xi' + + type1 = 'CPR' + type2 = 'GEF' + + with pytest.raises(ValueError): + get_dictkey_for_2pttype(type1, type2) + + +# # A mock datablock object for testing +# class MockDatablock: +# def has_section(self, section): +# return section in ['section1', 'section2'] + +# # A mock spectrum_array_from_block function for testing +# def mock_spectrum_array_from_block(data, section, types, xlabel, binformat): +# return [1, 2, 3], [4, 5, 6], [0.5, 1.0, 1.5], True, [0.4, 0.9, 1.4], [0.6, 1.1, 1.6] + +# def test_get_twoptdict_from_pipelinedata(): +# data = MockDatablock() +# original_spectrum_array_from_block = your_module.spectrum_array_from_block +# your_module.spectrum_array_from_block = mock_spectrum_array_from_block + +# result = get_twoptdict_from_pipelinedata(data) + +# your_module.spectrum_array_from_block = original_spectrum_array_from_block + +# assert 'galaxy_position_real_theta' in result +# assert 'galaxy_position_real_xi' in result +# assert 'galaxy_position_real_xi_bins' in result +# assert 'galaxy_position_real_xi_binavg' in result +# assert 'galaxy_position_real_theta_mins' in result +# assert 'galaxy_position_real_theta_maxs' in result + +# assert 'cmb_kappa_real_ell' in result +# assert 'cmb_kappa_real_cl' in result +# assert 'cmb_kappa_real_cl_bins' in result +# assert 'cmb_kappa_real_cl_binavg' in result +# assert 'cmb_kappa_real_ell_mins' in result +# assert 'cmb_kappa_real_ell_maxs' in result \ No newline at end of file