diff --git a/.felt/blank-input-tiles-grid-4-b9a61358.md b/.felt/blank-input-tiles-grid-4-b9a61358.md new file mode 100644 index 000000000..77c7cb6a8 --- /dev/null +++ b/.felt/blank-input-tiles-grid-4-b9a61358.md @@ -0,0 +1,9 @@ +--- +title: Blank input tiles grid_4 +depends-on: + - snakemake-crash-runtimeerror-20c2564f +created-at: 2026-07-08T17:09:48.720667191+02:00 +outcome: '12 of 37 input tiles in /n09data/hervas/skills_out/*_grid_4/images/SP_tiles are all-zero (image AND weight), identical set in all 4 shear sims: 276.280 276.281 277.282 280.280 280.282 281.281 282.281 284.279 285.279 285.280 286.279 286.280. Upstream sim-generation problem, needs regeneration by Fabian (hervas). Workaround: new tile_IDs_exclude config key (Snakefile + info.py), active in run config.yaml. Also added: log directive err_job_{tile}_{bit}.txt with tee so errors survive failure; config.yaml no longer tracked input of rule init (edits cascaded into full re-run); standard command now uses --rerun-triggers=mtime and fixed -j 4.' +--- + +Snakemake pipeline_all failed at job 512 (vignetmaker) for tile 281.281 in two sims: SExtractor (job 256) succeeds but finds 0 objects, writes header-only sexcat, vignetmaker raises ValueError empty/corrupt catalogue. Snakemake deletes outputs of failed jobs, so the shapepipe log vanished; real error was in tiles///output/run_sp_tile_ViVi_*/vignetmaker_runner_run_1/logs/process-*.log diff --git a/.felt/snakemake-crash-runtimeerror-20c2564f.md b/.felt/snakemake-crash-runtimeerror-20c2564f.md new file mode 100644 index 000000000..bd8d4d04a --- /dev/null +++ b/.felt/snakemake-crash-runtimeerror-20c2564f.md @@ -0,0 +1,9 @@ +--- +title: 'Snakemake crash: RuntimeError can''t start new thread during pipeline_all' +depends-on: + - shapepipe +created-at: 2026-07-08T09:32:38.198900721+02:00 +outcome: 'Root cause: user thread/process limit exhaustion. ulimit -u soft=1200 hard=1280 on candide. Snakemake''s scheduler Timer thread was just the victim; the real consumers are the N parallel apptainer ShapePipe jobs. apptainer_noslurm.sh strips only SLURM/PMI/PMIX/OMPI vars, so host OMP_NUM_THREADS leaks into every container job: -j $OMP_NUM_THREADS gives N jobs x N BLAS/OpenMP threads each ~ N^2 threads. Fix (agreed, NOT yet applied): in scripts/image_sims_pipeline/Snakefile add --env OMP_NUM_THREADS=1,OPENBLAS_NUM_THREADS=1,MKL_NUM_THREADS=1,NUMEXPR_NUM_THREADS=1 to APP_CMD and SPV_CMD, and run with a fixed conservative -j 4or5 instead of $OMP_NUM_THREADS. The completed job''s log file is valid, so rerunning snakemake resumes where it stopped.' +--- + +Run: snakemake -s Snakefile --configfile config.yaml -j $OMP_NUM_THREADS pipeline_all in ~/v2.0/image_sims/grids on candide (node c02). Crash in snakemake 9.19.0 job_scheduler._schedule_reevalutation -> threading.Timer().start() right after tile 288.279 job 16 completed. diff --git a/auxdir/CFIS/im_sims_202606/tile_numbers_1.txt b/auxdir/CFIS/im_sims_202606/tile_numbers_1.txt new file mode 100644 index 000000000..30673d233 --- /dev/null +++ b/auxdir/CFIS/im_sims_202606/tile_numbers_1.txt @@ -0,0 +1,40 @@ +233.293 +237.292 +238.292 +239.291 +239.292 +239.293 +240.291 +240.292 +240.293 +241.293 +242.290 +242.291 +242.292 +242.293 +243.290 +243.291 +243.292 +243.293 +244.293 +245.290 +245.291 +245.292 +246.290 +246.291 +247.289 +247.291 +247.292 +247.293 +248.288 +248.289 +248.290 +248.291 +249.289 +249.290 +249.291 +249.292 +250.289 +250.290 +250.292 +251.287 diff --git a/auxdir/CFIS/im_sims_202606/tile_numbers_2.txt b/auxdir/CFIS/im_sims_202606/tile_numbers_2.txt new file mode 100644 index 000000000..449bd74d6 --- /dev/null +++ b/auxdir/CFIS/im_sims_202606/tile_numbers_2.txt @@ -0,0 +1,43 @@ +251.288 +251.289 +251.291 +252.287 +252.288 +252.289 +252.291 +253.287 +253.288 +253.291 +254.286 +254.287 +254.288 +255.286 +255.288 +255.290 +256.286 +256.287 +257.285 +257.286 +257.287 +257.288 +258.287 +258.289 +259.285 +259.286 +259.287 +259.288 +259.289 +260.284 +260.285 +260.286 +260.287 +260.288 +260.289 +261.284 +261.285 +261.286 +261.287 +261.288 +262.285 +262.286 +262.287 diff --git a/auxdir/CFIS/im_sims_202606/tile_numbers_3.txt b/auxdir/CFIS/im_sims_202606/tile_numbers_3.txt new file mode 100644 index 000000000..d51361b45 --- /dev/null +++ b/auxdir/CFIS/im_sims_202606/tile_numbers_3.txt @@ -0,0 +1,46 @@ +262.288 +263.283 +263.284 +263.285 +263.286 +264.283 +264.284 +264.286 +264.287 +265.284 +265.285 +265.286 +265.287 +266.282 +266.283 +266.285 +266.286 +267.283 +267.284 +267.285 +267.286 +268.282 +268.283 +268.284 +268.285 +269.281 +269.282 +269.283 +269.284 +269.286 +270.281 +270.282 +270.283 +271.281 +271.282 +271.283 +271.284 +271.285 +272.281 +272.282 +272.283 +272.285 +273.280 +273.282 +273.283 +273.284 diff --git a/auxdir/CFIS/im_sims_202606/tile_numbers_4.txt b/auxdir/CFIS/im_sims_202606/tile_numbers_4.txt new file mode 100644 index 000000000..90077aea5 --- /dev/null +++ b/auxdir/CFIS/im_sims_202606/tile_numbers_4.txt @@ -0,0 +1,37 @@ +274.279 +274.280 +274.282 +274.283 +275.279 +275.280 +275.282 +275.283 +276.279 +276.280 +276.281 +276.282 +276.283 +277.279 +277.281 +277.282 +278.279 +278.281 +278.282 +279.281 +279.282 +280.280 +280.281 +280.282 +281.280 +281.281 +282.280 +282.281 +283.280 +284.279 +285.279 +285.280 +286.279 +286.280 +287.279 +288.279 +289.279 diff --git a/auxdir/CFIS/im_sims_202606/tile_numbers_5.txt b/auxdir/CFIS/im_sims_202606/tile_numbers_5.txt new file mode 100644 index 000000000..e69de29bb diff --git a/bin/shapepipe_run.py b/bin/shapepipe_run.py index dedfc5093..cacb6772b 100755 --- a/bin/shapepipe_run.py +++ b/bin/shapepipe_run.py @@ -1,3 +1,5 @@ #!/usr/bin/env python from shapepipe.shapepipe_run import main +print("MKDEBUG start") main() +print("MKDEBUG end") diff --git a/docs/source/image_sims_calibration.md b/docs/source/image_sims_calibration.md new file mode 100644 index 000000000..e1b550318 --- /dev/null +++ b/docs/source/image_sims_calibration.md @@ -0,0 +1,365 @@ +# Image Simulations Calibration Pipeline + +## Overview + +The image simulations calibration pipeline derives multiplicative (m) and additive (c) shear bias corrections from synthetic image simulations. It uses ShapePipe to measure galaxy shapes in five pre-sheared grids, then computes bias by comparing measurements across shear directions. + +**Key outputs:** +- `m_bias_results.yaml` / `m_bias_results.txt` — final m₁, m₂, c₁, c₂ bias estimates with bootstrap errors +- `mbias_cumulative.yaml` — convergence history (m/c as function of n_tiles) +- `mbias_convergence.png` — m and c vs tile count (2 panels) +- `mbias_errors.png` — error shrinkage vs tile count (2 panels) + +--- + +## Quick Start + +### Prerequisites +- container: `/n17data/mkilbing/sp_validation_im_sims.sif` +- Pipeline scripts: `/path/to/shapepipe/scripts/image_sims_pipeline/` + +### Setup (one time) + +```bash +# Copy template config to run directory and rename to config.yaml +cp /path/to/shapepipe/scripts/image_sims_pipeline/config.yaml.template /config.yaml + +# Edit config.yaml with your tile IDs, grid number, etc. +``` + +### Run (fully automatic, snakemake) + +To get help on how to run snakemake, and the different targets and options, type +```bash +python ~/shapepipe/scripts/image_sims_pipeline/info.py -h +``` + +To monitor the status of the run and completed number of tiles, type +```bash +python ~/shapepipe/scripts/image_sims_pipeline/info.py -m [-v] +``` + +First, for convenience define the short cut +```bash +alias sn_ims="snakemake -s ~/shapepipe/scripts/image_sims_pipeline/Snakefile --configfile config.yaml -j $OMP_NUM_THREADS' +``` + +To go step by step (one target at a time), do + +```bash +# Initialise subdirs (one for each sheared verions) +sn_ims init_all + +# Optional: scan input weight images for tile coverage. Tiles with a low +fraction of non-zero weight pixels (<50%, can be set by the user) has a small +number of detected galaxies. In some cases this number can be zero which would +lead to a failure of the ShapePipe module vignet_maker_runner. + +This call creates a tile-by-tile coverage report and adds the below-threshold +tiles to an exclusion list in config.yaml. + +# Run shapepipe +sn_ims pipeline_all +``` + +```bash + +# Full pipeline to calibrated catalogues +snakemake -s ~/astro/repositories/github/shapepipe/scripts/image_sims_pipeline/Snakefile \ + --configfile config.yaml -j 5 calibrate_all + +# M-bias computation with convergence tracking +snakemake -s ~/astro/repositories/github/shapepipe/scripts/image_sims_pipeline/Snakefile \ + --configfile config.yaml -j 1 mbias + +# Monitor progress +~/astro/repositories/github/shapepipe/scripts/image_sims_pipeline/info.py -m -v + +# Incremental m-bias as tiles finish +~/astro/repositories/github/shapepipe/scripts/image_sims_pipeline/monitor_mbias.py -v +``` + +### Run (non-automatic) + +Activate Docker container with +```bash +app_gen --bind /n09data /n17data/mkilbing/shapepipe_im_sims.sif +``` + +Get config file +```bash +cp /path/to/shapepipe/scripts/image_sims_pipeline/config.yaml.template ./config.yaml +``` + +Edit config file. + +Initialise subdir: +```bash +init_run_v2.0.py [-s ] +``` +with e.g. `subdir=1p2z_grid_1`. + +Start pipeline runs + +```bash +cd +run_job_sp_canfar_v2.0.bash -t image_sims -j 1 +``` + +This will run the first job (`-j 1`) for all tiles found in `tile_numbers.txt`. + +To get help on the options, use the flag `-h`. + +To test it one a single tile for testing, add the option `-e ID`, e.g. `-e 272.281`. + +To run all pipeline jobs, either call `run_job_sp_canfar_v2.0.bash` for each subsequent job, +e.g. `-j 1`, `-j 2`, `-j 4`, or group them by (bit-)adding the job numbers. The entire pipeline +will be run with `-j 4095`. + + +## Configuration + +### Snakemake Config: `config.yaml` + +```yaml +# Run directory +base: /n17data/mkilbing/astro/Runs/shapepipe/CFIS/v2.0/image_sims + +# Tile IDs (all tiles in the survey, or subset for testing) +tile_IDs: ["254.286", "254.287"] + +# catalogue subversion (mask config: mask_v1.X.{sample}_im_sim.yaml) +sample: 9 + +# Simulation type (grid or others) +type: grid + +# Grid number (e.g., grid_2) +num: 2 + +# SMP batch size per ShapePipe job (-1 = use container default) +n_smp: -1 +``` + +### M-bias Calibration Settings + +M-bias computation parameters are hard-coded in the Snakefile `mbias` rule: + +```python +shear_amplitude: 0.02 # Input shear amplitude (absolute value) +match_radius_deg: 0.0002 # Position matching radius for cross-matching +w_col: w_des # Weight column name +n_bootstrap: 500 # Bootstrap resamples for error estimation +catalog_name: shape_catalog_cut_ngmix.fits +``` + +To modify these parameters for a different run, edit the `mbias` rule in the Snakefile. + +--- + +## Pipeline Stages + +The pipeline runs in order, with each stage dependent on the previous: + +### 1. `init_all` — Initialize run directories + +Creates per-sim directories, parameter files (`params.py`), and mask configs. + +**Outputs:** `{grids}/{sim}/params.py`, `{grids}/{sim}/config_mask.yaml` + +### 2. `pipeline_all` — Run ShapePipe on all tiles + +Executes ShapePipe's full job sequence (bits 1→2→4→...→2048) for each tile in each simulation. + +**Note:** `-j 5` means 5 concurrent Snakemake jobs (tile×sim pairs run sequentially within each job). + +**Outputs:** ShapePipe catalogues, logs at `{grids}/{sim}/logs/log_job_{tile}_{bit}.txt` + +**Time:** ~hours per grid depending on tile count and cluster load. + +### 3. `merge_all` — Merge ShapePipe tile outputs + +Combines per-tile catalogs into a single HDF5 file per simulation using `create_final_cat.py`. + +**Outputs:** `{grids}/{sim}/final_cat_{sim}.hdf5` (HDF5 with tile counts as attributes) + +### 4. `extract_all` — Extract comprehensive catalogues + +Reads merged HDF5 and extracts shape information, PSF quantities, and pre-calibration columns using `extract_info.py`. + +**Outputs:** `{grids}/{sim}/shape_catalog_comprehensive_ngmix.hdf5` (with `n_tiles` attribute in header) + +### 5. `calibrate_all` — Apply cuts and calibrate + +Applies selection masks (flags, magnitude, signal-to-noise) and computes shear calibration (m/c per object) using `calibrate_comprehensive_cat.py`. + +**Outputs:** `{grids}/{sim}/shape_catalog_cut_ngmix.fits` (final catalogue, ready for m-bias) + +### 6. `mbias` — Compute m-bias with convergence tracking + +Computes m₁, m₂, c₁, c₂ from the five shear pairs. With `--cumulative`, tracks convergence as tile counts grow. + +**Outputs:** +- `results/m_bias_results.yaml` — final results (YAML) +- `results/m_bias_results.txt` — final results (human-readable text) +- `results/mbias_cumulative.yaml` — convergence history +- `results/mbias_convergence.png` — m/c vs n_tiles (2 panels) +- `results/mbias_errors.png` — error convergence (2 panels) + +--- + +## Monitoring Convergence + +### Live Status: `info.py` + +Monitor pipeline progress with a status table: + +```bash +cd /path/to/run/image_sims +~/astro/repositories/github/shapepipe/scripts/image_sims_pipeline/info.py -m -v +``` + +**Status table (example output):** + +| Simulation | #Jobs | Job Bits | Merge | Extract | Calibrate | +|-----------|-------|----------|-------|---------|-----------| +| 1m2z_grid_2 | 12/12 | ✓✓✓✓✓✓✓✓✓✓✓✓ | 2 | ✓ | ✓ | +| 1p2z_grid_2 | 12/12 | ✓✓✓✓✓✓✓✓✓✓✓✓ | 2 | ✓ | ✓ | +| 1z2m_grid_2 | 12/12 | ✓✓✓✓✓✓✓✓✓✓✓✓ | 2 | ✓ | ✓ | +| 1z2p_grid_2 | 12/12 | ✓✓✓✓✓✓✓✓✓✓✓✓ | 2 | ✓ | ✓ | +| 1z2z_grid_2 | 12/12 | ✓✓✓✓✓✓✓✓✓✓✓✓ | 2 | ✓ | ✓ | + +**Column meanings:** +- **#Jobs** — ShapePipe job completion count (max 12 bits) +- **Job Bits** — per-bit status (✓ = all tiles done, · = incomplete) +- **Merge** — tile count in final HDF5 catalogue +- **Extract** — comprehensive catalog status +- **Calibrate** — calibrated catalogue status + +### Incremental M-bias: `monitor_mbias.py` + +Recompute m-bias each time new tiles finish all ShapePipe jobs. Automatically detects completed tiles and updates convergence tracking. + +```bash +cd /path/to/run/image_sims +~/astro/repositories/github/shapepipe/scripts/image_sims_pipeline/monitor_mbias.py -v +``` + +Run this repeatedly as tiles complete to watch m-bias converge. + +--- + +## Results Interpretation + +### M-Bias Results Files + +**YAML format** (`m_bias_results.yaml`): +```yaml +m1: -0.7909173336321427 +m1_err: 0.27407603929732716 +c1: -0.0030413689177988708 +c1_err: 0.005899515490890933 +m2: -1.350016719584532 +m2_err: 0.20884008300375373 +c2: -0.004950048084875509 +c2_err: 0.004208875641566414 +``` + +**Text format** (`m_bias_results.txt`): +``` +Multiplicative and additive shear bias from image simulations +============================================================ + +m1 = -0.790917 ± 0.274076 +c1 = -0.003041 ± 0.005900 + +m2 = -1.350017 ± 0.208840 +c2 = -0.004950 ± 0.004209 + +Errors computed via bootstrap resampling (n=500 resamples) +``` + +### Convergence History (`mbias_cumulative.yaml`) + +Tracks m/c evolution as tile count grows: + +```yaml +'2': + c1: -0.0030413689177988708 + c1_err: 0.005899515490890933 + m1: -0.7909173336321427 + m1_err: 0.27407603929732716 + # ... (c2, m2, errors) +'4': + c1: -0.0025... + # ... (more tiles) +``` + +### Understanding Error Bars + +Errors are computed via **bootstrap resampling**: + +1. Draw N=500 random resamples (with replacement) from measured galaxy ellipticities +2. Recompute m and c for each resample +3. Error = standard deviation of the bootstrap distribution + +This captures: +- **Photometric noise** — measurement uncertainties per galaxy +- **Cosmic variance** — shape correlations across the sky +- **Calibration uncertainties** — from shear responsivity scatter + +**Error shrinkage:** Error ∝ 1/√(n_tiles). More tiles → tighter constraints. + +### Convergence Plots + +**`mbias_convergence.png`** — m and c with error bars vs n_tiles +- Left panel: multiplicative bias (m₁, m₂) with error bars +- Right panel: additive bias (c₁, c₂) with error bars +- Shows systematic trends and statistical uncertainties + +**`mbias_errors.png`** — error shrinkage vs n_tiles +- Left panel: m error bars (m₁, m₂) only +- Right panel: c error bars (c₁, c₂) only +- Shows constraint tightening as data accumulates +- Useful for deciding when m-bias is "converged" (errors small enough for science) + +--- + +## File Structure + +``` +/n17data/mkilbing/astro/Runs/shapepipe/CFIS/v2.0/image_sims/ +├── config.yaml # Run configuration +├── .snakemake/ # Snakemake metadata +├── logs/ # Snakemake logs +├── grids/ +│ ├── {sim}_grid_2/ +│ │ ├── logs/ # ShapePipe job logs +│ │ ├── tiles/ # Per-tile ShapePipe outputs +│ │ ├── final_cat_{sim}.hdf5 # Merged catalogue +│ │ ├── shape_catalog_comprehensive_ngmix.hdf5 +│ │ └── shape_catalog_cut_ngmix.fits +│ └── results/ +│ ├── m_bias_results.yaml +│ ├── m_bias_results.txt +│ ├── mbias_cumulative.yaml +│ ├── mbias_convergence.png +│ └── mbias_errors.png +└── monitoring/ # Incremental m-bias workspace + └── {sim}/ + ├── final_cat_{sim}.hdf5 + ├── shape_catalog_comprehensive_ngmix.hdf5 + └── shape_catalog_cut_ngmix.fits +``` + +--- + +## References + +- **Pipeline code:** `~/astro/repositories/github/shapepipe/scripts/image_sims_pipeline/` +- **ShapePipe:** `~/astro/repositories/github/shapepipe/` +- **SP Validation:** `~/astro/repositories/github/sp_validation/` + +--- + +*Last updated: 2026-06-26* diff --git a/example/cfis/config_tile_Ng_template.ini b/example/cfis/config_tile_Ng_template.ini index 5bfbc2c79..6214765fc 100644 --- a/example/cfis/config_tile_Ng_template.ini +++ b/example/cfis/config_tile_Ng_template.ini @@ -83,5 +83,14 @@ PIXEL_SCALE = 0.186 # noisy for stars and flagged as incorrect — see #767). CENTROID_SOURCE = wcs +# SEED_FROM_POSITION (optional): FALSE (default) seeds one RNG per tile. +# TRUE seeds a per-object RNG from the object's sky position instead, so that +# metacal's fixnoise counter-noise and the fit guesses are identical for the +# same object across Pujol image-simulation shear branches and cancel in the +# branch difference — shrinking the m-bias error. This is FOR IMAGE +# SIMULATIONS ONLY (Pujol noise cancellation, ngmix#796); leave it off (or +# omit) for real data, where it has no benefit. +# SEED_FROM_POSITION = FALSE + ID_OBJ_MIN = X ID_OBJ_MAX = X diff --git a/example/cfis_image_sims/config_exp_Gie_symlink.ini b/example/cfis_image_sims/config_exp_Gie_symlink.ini new file mode 100644 index 000000000..a6a3aed92 --- /dev/null +++ b/example/cfis_image_sims/config_exp_Gie_symlink.ini @@ -0,0 +1,90 @@ +# ShapePipe configuration file for: get images + + +## Default ShapePipe options +[DEFAULT] + +# verbose mode (optional), default: True, print messages on terminal +VERBOSE = False + +# Name of run (optional) default: shapepipe_run +RUN_NAME = run_sp_exp_Gie + +# Add date and time to RUN_NAME, optional, default: False +RUN_DATETIME = True + + +## ShapePipe execution options +[EXECUTION] + +# Module name, single string or comma-separated list of valid module runner names +MODULE = get_images_runner + +# Parallel processing mode, SMP or MPI +MODE = SMP + + +## ShapePipe file handling options +[FILE] + +# Log file master name, optional, default: shapepipe +LOG_NAME = log_sp + +# Runner log file name, optional, default: shapepipe_runs +RUN_LOG_NAME = log_run_sp + +# Input directory, containing input files, single string or list of names +INPUT_DIR = $SP_RUN + +# Output directory +OUTPUT_DIR = $SP_RUN/output + + +## ShapePipe job handling options +[JOB] + +# Batch size of parallel processing (optional), default is 1, i.e. run all jobs in serial +SMP_BATCH_SIZE = 1 + +# Timeout value (optional), default is None, i.e. no timeout limit applied +TIMEOUT = 96:00:00 + + +## Module options + +# Get exposures +[GET_IMAGES_RUNNER] + +INPUT_MODULE = last:find_exposures_runner + +FILE_PATTERN = exp_numbers + +FILE_EXT = .txt + +# NUMBERING_SCHEME (optional) string with numbering pattern for input files +NUMBERING_SCHEME = -000-000 + + +# Paths + +# Output path (optional, default is [FILE]:OUTPUT_DIR +# OUTPUT_PATH = input_images + +# Input path where original images are stored. Can be local path or vos url. +# Single string or list of strings +INPUT_PATH = $SP_DIR/input_exp/,$SP_DIR/input_exp/,$SP_DIR/input_exp/ + +# Input file pattern including tile number as dummy template +INPUT_FILE_PATTERN = simu_image-0000000,simu_weight-0000000,simu_flag-0000000 + +# Input file extensions +INPUT_FILE_EXT = .fits, .fits,.fits + +# Input numbering scheme, python regexp +INPUT_NUMBERING = \d{7} + +# Output file pattern without number +OUTPUT_FILE_PATTERN = image-,weight-,flag- + +# Method to retrieve images, one in 'vos', 'symlink' +RETRIEVE = symlink diff --git a/example/cfis_image_sims/config_exp_Ma_onthefly.ini b/example/cfis_image_sims/config_exp_Ma_onthefly.ini new file mode 100644 index 000000000..df0307a19 --- /dev/null +++ b/example/cfis_image_sims/config_exp_Ma_onthefly.ini @@ -0,0 +1,76 @@ +# ShapePipe configuration file for masking of exposures + + +## Default ShapePipe options +[DEFAULT] + +# verbose mode (optional), default: True, print messages on terminal +VERBOSE = True + +# Name of run (optional) default: shapepipe_run +RUN_NAME = run_sp_exp_Ma + +# Add date and time to RUN_NAME, optional, default: False +; RUN_DATETIME = False + + +## ShapePipe execution options +[EXECUTION] + +# Module name, single string or comma-separated list of valid module runner names +MODULE = mask_runner + +# Parallel processing mode, SMP or MPI +MODE = SMP + + +## ShapePipe file handling options +[FILE] + +# Log file master name, optional, default: shapepipe +LOG_NAME = log_sp + +# Runner log file name, optional, default: shapepipe_runs +RUN_LOG_NAME = log_run_sp + +# Input directory, containing input files, single string or list of names +INPUT_DIR = . + +# Output directory +OUTPUT_DIR = $SP_RUN/output + + +## ShapePipe job handling options +[JOB] + +# Batch size of parallel processing (optional), default is 1, i.e. run all jobs in serial +SMP_BATCH_SIZE = 1 + +# Timeout value (optional), default is None, i.e. no timeout limit applied +TIMEOUT = 96:00:00 + + +## Module options + +### Mask exposures +[MASK_RUNNER] + +# Parent module +INPUT_DIR = last:split_exp_runner + +# Update numbering convention, accounting for HDU number of +# single-exposure single-HDU files +NUMBERING_SCHEME = -0000000-0 + +# Path of mask config file +MASK_CONFIG_PATH = $SP_CONFIG/config_onthefly.mask_simu + +# External mask file flag, use if True, otherwise ignore +USE_EXT_FLAG = True + +# External star catalogue flag, use external cat if True, +# obtain from online catalogue if False +USE_EXT_STAR = False + +# File name suffix for the output flag files (optional) +PREFIX = pipeline diff --git a/example/cfis_image_sims/config_exp_Sp.ini b/example/cfis_image_sims/config_exp_Sp.ini new file mode 100644 index 000000000..ca010b643 --- /dev/null +++ b/example/cfis_image_sims/config_exp_Sp.ini @@ -0,0 +1,74 @@ +# ShapePipe configuration file for single-exposures, +# split images + + +## Default ShapePipe options +[DEFAULT] + +# verbose mode (optional), default: True, print messages on terminal +VERBOSE = True + +# Name of run (optional) default: shapepipe_run +RUN_NAME = run_sp_exp_Sp + +# Add date and time to RUN_NAME, optional, default: True +RUN_DATETIME = True + + +## ShapePipe execution options +[EXECUTION] + +# Module name, single string or comma-separated list of valid module runner names +MODULE = split_exp_runner + +# Run mode, SMP or MPI +MODE = SMP + + +## ShapePipe file handling options +[FILE] + +# Log file master name, optional, default: shapepipe +LOG_NAME = log_sp + +# Runner log file name, optional, default: shapepipe_runs +RUN_LOG_NAME = log_run_sp + +# Input directory, containing input files, single string or list of names with length matching FILE_PATTERN +INPUT_DIR = . + +# Output directory +OUTPUT_DIR = $SP_RUN/output + + +## ShapePipe job handling options +[JOB] + +# Batch size of parallel processing (optional), default is 1, i.e. run all jobs in serial +SMP_BATCH_SIZE = 1 + +# Timeout value (optional), default is None, i.e. no timeout limit applied +TIMEOUT = 96:00:00 + + +## Module options + +[SPLIT_EXP_RUNNER] + +INPUT_DIR = last:get_images_runner + +FILE_PATTERN = image, weight, flag + +# Matches compressed single-exposure files +FILE_EXT = .fits, .fits, .fits + +NUMBERING_SCHEME = -0000000 + +# OUTPUT_SUFFIX, actually file name prefixes. +# Expected keyword "flag" will lead to a behavior where the data are saved as int. +# The code also expects the image data to use the "image" suffix +# (default value in the pipeline). +OUTPUT_SUFFIX = image, weight, flag + +# Number of HDUs/CCDs of mosaic +N_HDU = 40 diff --git a/example/cfis_image_sims/config_exp_psfex.ini b/example/cfis_image_sims/config_exp_psfex.ini new file mode 100644 index 000000000..f44bb38d7 --- /dev/null +++ b/example/cfis_image_sims/config_exp_psfex.ini @@ -0,0 +1,54 @@ +# ShapePipe configuration file for fake PSF postage stamps (image simulations) + + +## Default ShapePipe options +[DEFAULT] + +VERBOSE = True + +RUN_NAME = run_sp_tile_fpsf + +RUN_DATETIME = True + + +## ShapePipe execution options +[EXECUTION] + +MODULE = fake_psf_runner + +MODE = SMP + + +## ShapePipe file handling options +[FILE] + +LOG_NAME = log_sp + +RUN_LOG_NAME = log_run_sp + +INPUT_DIR = $SP_RUN/output + +OUTPUT_DIR = $SP_RUN/output + + +## ShapePipe job handling options +[JOB] + +SMP_BATCH_SIZE = 1 + +TIMEOUT = 96:00:00 + + +## Module options +[FAKE_PSF_RUNNER] + +INPUT_MODULE = sextractor_runner + +FILE_PATTERN = sexcat + +FILE_EXT = .fits + +NUMBERING_SCHEME = -000-000 + +# Path to the pickled PSF dictionary +PSF_DICT_PATH = /home/hervas/fhervas/workdir_skills/input/psf_files/Full_psf_dict.pickle diff --git a/example/cfis_image_sims/config_onthefly.mask_simu b/example/cfis_image_sims/config_onthefly.mask_simu new file mode 100644 index 000000000..1a63cc2e5 --- /dev/null +++ b/example/cfis_image_sims/config_onthefly.mask_simu @@ -0,0 +1,86 @@ +# Mask module configuration file for single-exposure images + +## Paths to executables +[PROGRAM_PATH] + +WW_PATH = weightwatcher +WW_CONFIG_FILE = $SP_CONFIG/mask_default/default.ww + +# Indicate cds client executable if no external star catalogue is available +# (e.g. no internet access on run nodes) +CDSCLIENT_PATH = findgsc2.2 + + +## Border mask +[BORDER_PARAMETERS] + +BORDER_MAKE = True + +BORDER_WIDTH = 50 +BORDER_FLAG_VALUE = 4 + + +## Halo mask +[HALO_PARAMETERS] + +HALO_MAKE = False + +HALO_MASKMODEL_PATH = $SP_CONFIG/mask_default/halo_mask.reg +HALO_MAG_LIM = 13. +HALO_SCALE_FACTOR = 0.05 +HALO_MAG_PIVOT = 13.8 +HALO_FLAG_VALUE = 2 +HALO_REG_FILE = halo.reg + + +## Diffraction spike mask +[SPIKE_PARAMETERS] + +SPIKE_MAKE = False + +SPIKE_MASKMODEL_PATH = $SP_CONFIG/mask_default/MEGAPRIME_star_i_13.8.reg +SPIKE_MAG_LIM = 18. +SPIKE_SCALE_FACTOR = 0.3 +SPIKE_MAG_PIVOT = 13.8 +SPIKE_FLAG_VALUE = 128 +SPIKE_REG_FILE = spike.reg + + +## Messier mask +[MESSIER_PARAMETERS] + +MESSIER_MAKE = False + +MESSIER_CAT_PATH = $SP_CONFIG/mask_default/Messier_catalog_updated.fits +MESSIER_SIZE_PLUS = 0. +MESSIER_FLAG_VALUE = 16 + + +## NGC mask +[NGC_PARAMETERS] + +NGC_MAKE = False + +NGC_CAT_PATH = $SP_CONFIG/mask_default/ngc_cat.fits +NGC_SIZE_PLUS = 0. +NGC_FLAG_VALUE = 32 + + + +## Missing data parameters +[MD_PARAMETERS] + +MD_MAKE = False + +MD_THRESH_FLAG = 0.3 +MD_THRESH_REMOVE = 0.75 +MD_REMOVE = False + + +## Other parameters +[OTHER] + +TEMP_DIRECTORY = .temp + +KEEP_REG_FILE = False +KEEP_INDIVIDUAL_MASK = False diff --git a/example/cfis_image_sims/config_tile_Fe.ini b/example/cfis_image_sims/config_tile_Fe.ini new file mode 100644 index 000000000..9c3cb8cd5 --- /dev/null +++ b/example/cfis_image_sims/config_tile_Fe.ini @@ -0,0 +1,68 @@ +# ShapePipe configuration file for: find exposures + + +## Default ShapePipe options +[DEFAULT] + +# verbose mode (optional), default: True, print messages on terminal +VERBOSE = False + +# Name of run (optional) default: shapepipe_run +RUN_NAME = run_sp_tile_Fe + +# Add date and time to RUN_NAME, optional, default: False +RUN_DATETIME = True + + +## ShapePipe execution options +[EXECUTION] + +# Module name, single string or comma-separated list of valid module runner names +MODULE = find_exposures_runner + +# Parallel processing mode, SMP or MPI +MODE = SMP + + +## ShapePipe file handling options +[FILE] + +# Log file master name, optional, default: shapepipe +LOG_NAME = log_sp + +# Runner log file name, optional, default: shapepipe_runs +RUN_LOG_NAME = log_run_sp + +# Input directory, containing input files, single string or list of names +INPUT_DIR = $SP_RUN + +# Output directory +OUTPUT_DIR = $SP_RUN/output + + +## ShapePipe job handling options +[JOB] + +# Batch size of parallel processing (optional), default is 1, i.e. run all jobs in serial +SMP_BATCH_SIZE = 1 + +# Timeout value (optional), default is None, i.e. no timeout limit applied +TIMEOUT = 96:00:00 + + +## Module options + +[FIND_EXPOSURES_RUNNER] + +INPUT_MODULE = last:get_images_runner + +FILE_PATTERN = CFIS_simu_image + +FILE_EXT = .fits + +# NUMBERING_SCHEME (optional) string with numbering pattern for input files +NUMBERING_SCHEME = -000-000 + +COLNUM = 2 + +EXP_PREFIX = simu_image- diff --git a/example/cfis_image_sims/config_tile_Git_symlink.ini b/example/cfis_image_sims/config_tile_Git_symlink.ini new file mode 100644 index 000000000..87c4a8816 --- /dev/null +++ b/example/cfis_image_sims/config_tile_Git_symlink.ini @@ -0,0 +1,87 @@ +# ShapePipe configuration file for: get tile images + + +## Default ShapePipe options +[DEFAULT] + +# verbose mode (optional), default: True, print messages on terminal +VERBOSE = False + +# Name of run (optional) default: shapepipe_run +RUN_NAME = run_sp_tile_Git + +# Add date and time to RUN_NAME, optional, default: False +RUN_DATETIME = True + + +## ShapePipe execution options +[EXECUTION] + +# Module name, single string or comma-separated list of valid module runner names +MODULE = get_images_runner + +# Parallel processing mode, SMP or MPI +MODE = SMP + + +## ShapePipe file handling options +[FILE] + +# Log file master name, optional, default: shapepipe +LOG_NAME = log_sp + +# Runner log file name, optional, default: shapepipe_runs +RUN_LOG_NAME = log_run_sp + +# Input directory, containing input files, single string or list of names +INPUT_DIR = $SP_RUN + +# Output directory +OUTPUT_DIR = $SP_RUN/output + + +## ShapePipe job handling options +[JOB] + +# Batch size of parallel processing (optional), default is 1, i.e. run all jobs in serial +SMP_BATCH_SIZE = 1 + +# Timeout value (optional), default is None, i.e. no timeout limit applied +TIMEOUT = 96:00:00 + + +## Module options + +# Get tiles +[GET_IMAGES_RUNNER] + +FILE_PATTERN = tile_numbers + +FILE_EXT = .txt + +# NUMBERING_SCHEME (optional) string with numbering pattern for input files +NUMBERING_SCHEME = + +# Paths + +# Input path where original images are stored. Can be local path or vos url. +# Single string or list of strings +INPUT_PATH = $SP_DIR/input_tiles, $SP_DIR/input_tiles + +# Input file pattern including tile number as dummy template +INPUT_FILE_PATTERN = CFIS_simu_image-000-000, CFIS_simu_weight-000-000 + +# Input file extensions +INPUT_FILE_EXT = .fits, .fits + +# Input numbering scheme, python regexp +INPUT_NUMBERING = \d{3}-\d{3} + +# Output file pattern without number +OUTPUT_FILE_PATTERN = CFIS_simu_image-, CFIS_simu_weight- + +# Copy/download method, one in 'vos', 'symlink' +RETRIEVE = symlink + +# Copy command options, optional +RETRIEVE_OPTIONS = -L diff --git a/example/cfis_image_sims/config_tile_Ma_onthefly.ini b/example/cfis_image_sims/config_tile_Ma_onthefly.ini new file mode 100644 index 000000000..0f49eccea --- /dev/null +++ b/example/cfis_image_sims/config_tile_Ma_onthefly.ini @@ -0,0 +1,82 @@ +# ShapePipe configuration file for masking of tiles + + +## Default ShapePipe options +[DEFAULT] + +# verbose mode (optional), default: True, print messages on terminal +VERBOSE = True + +# Name of run (optional) default: shapepipe_run +RUN_NAME = run_sp_tile_Ma + +# Add date and time to RUN_NAME, optional, default: False +; RUN_DATETIME = False + + +## ShapePipe execution options +[EXECUTION] + +# Module name, single string or comma-separated list of valid module runner names +MODULE = mask_runner + +# Parallel processing mode, SMP or MPI +MODE = SMP + + +## ShapePipe file handling options +[FILE] + +# Log file master name, optional, default: shapepipe +LOG_NAME = log_sp + +# Runner log file name, optional, default: shapepipe_runs +RUN_LOG_NAME = log_run_sp + +# Input directory, containing input files, single string or list of names +INPUT_DIR = $SP_RUN/output + +# Output directory +OUTPUT_DIR = $SP_RUN/output + + +## ShapePipe job handling options +[JOB] + +# Batch size of parallel processing (optional), default is 1, i.e. run all jobs in serial +SMP_BATCH_SIZE = 1 + +# Timeout value (optional), default is None, i.e. no timeout limit applied +TIMEOUT = 96:00:00 + + +## Module options + +### Mask tiles +[MASK_RUNNER] + +# Input directory, containing input files, single string or list of names +INPUT_DIR = run_sp_Git:get_images_runner, last:uncompress_fits_runner + +# NUMBERING_SCHEME (optional) string with numbering pattern for input files +NUMBERING_SCHEME = -000-000 + +# Input file pattern(s), list of strings with length matching number of expected input file types +# Cannot contain wild cards +FILE_PATTERN = CFIS_simu_image, CFIS_simu_weight + +# FILE_EXT (optional) list of string extensions to identify input files +FILE_EXT = .fits, .fits + +# Path of mask config file +MASK_CONFIG_PATH = $SP_CONFIG/config_tile_onthefly.mask_simu + +# External mask file flag, use if True, otherwise ignore +USE_EXT_FLAG = False + +# External star catalogue flag, use external cat if True, +# obtain from online catalogue if False +USE_EXT_STAR = False + +# File name suffix for the output flag files (optional) +PREFIX = pipeline diff --git a/example/cfis_image_sims/config_tile_Mc_psfex.ini b/example/cfis_image_sims/config_tile_Mc_psfex.ini new file mode 100644 index 000000000..3b186484e --- /dev/null +++ b/example/cfis_image_sims/config_tile_Mc_psfex.ini @@ -0,0 +1,79 @@ +# ShapePipe post-run configuration file: create final catalogs for psfex + + +## Default ShapePipe options +[DEFAULT] + +# verbose mode (optional), default: True, print messages on terminal +VERBOSE = True + +# Name of run (optional) default: shapepipe_run +RUN_NAME = run_sp_tile_Mc_psfex + +# Add date and time to RUN_NAME, optional, default: True +; RUN_DATETIME = False + + +## ShapePipe execution options +[EXECUTION] + +# Module name, single string or comma-separated list of valid module runner names +MODULE = make_cat_runner + +# Parallel processing mode, SMP or MPI +MODE = SMP + + +## ShapePipe file handling options +[FILE] + +# Log file master name, optional, default: shapepipe +LOG_NAME = log_sp + +# Runner log file name, optional, default: shapepipe_runs +RUN_LOG_NAME = log_run_sp + +# Input directory, containing input files, single string or list of names with length matching FILE_PATTERN +INPUT_DIR = . + +# Output directory +OUTPUT_DIR = ./output + + +## ShapePipe job handling options +[JOB] + +# Batch size of parallel processing (optional), default is 1, i.e. run all jobs in serial +SMP_BATCH_SIZE = 1 + +# Timeout value (optional), default is None, i.e. no timeout limit applied +TIMEOUT = 96:00:00 + + +## Module options + +[MAKE_CAT_RUNNER] + +# Input directory, containing input files, single string or list of names with length matching FILE_PATTERN +INPUT_DIR = run_sp_tile_Sx:sextractor_runner, last:fake_psf_runner, last:ngmix_runner + +# Input file pattern(s), list of strings with length matching number of expected input file types +# Cannot contain wild cards +FILE_PATTERN = sexcat, galaxy_psf, ngmix + +# FILE_EXT (optional) list of string extensions to identify input files +FILE_EXT = .fits, .sqlite, .fits + +# Numbering convention, string that exemplifies a numbering pattern. +# Matches input single exposures (with 'p' removed) +# Needs to be given in this section, will be updated in module +# sections below +NUMBERING_SCHEME = -000-000 + +# Star/galaxy spread-model classification. The image-sims pipeline has no +# spread_model_runner stage (make_cat inputs are sextractor + fake_psf + +# ngmix only), so there is nothing to classify — disable it. The newer +# make_cat_runner reads this key unconditionally (NoOptionError if absent). +SM_DO_CLASSIFICATION = False + +SHAPE_MEASUREMENT_TYPE = ngmix diff --git a/example/cfis_image_sims/config_tile_Mh_exp.ini b/example/cfis_image_sims/config_tile_Mh_exp.ini new file mode 100644 index 000000000..94270a7e6 --- /dev/null +++ b/example/cfis_image_sims/config_tile_Mh_exp.ini @@ -0,0 +1,77 @@ +# ShapePipe configuration file for merging per-exposure WCS headers +# at the tile level. Input is the exp_numbers file produced by +# find_exposures_runner; EXP_BASE_DIR tells the runner where to find +# the per-exposure split_exp_runner header .npy files. + + +## Default ShapePipe options +[DEFAULT] + +# verbose mode (optional), default: True, print messages on terminal +VERBOSE = True + +# Name of run (optional) default: shapepipe_run +RUN_NAME = run_sp_tile_Mh_exp + +# Add date and time to RUN_NAME, optional, default: True +RUN_DATETIME = True + + +## ShapePipe execution options +[EXECUTION] + +# Module name, single string or comma-separated list of valid module runner names +MODULE = merge_headers_runner + +# Run mode, SMP or MPI +MODE = SMP + + +## ShapePipe file handling options +[FILE] + +# Log file master name, optional, default: shapepipe +LOG_NAME = log_sp + +# Runner log file name, optional, default: shapepipe_runs +RUN_LOG_NAME = log_run_sp + +# Input directory, containing input files, single string or list of names with length matching FILE_PATTERN +INPUT_DIR = . + +# Output directory +OUTPUT_DIR = $SP_RUN/output + + +## ShapePipe job handling options +[JOB] + +# Batch size of parallel processing (optional), default is 1, i.e. run all jobs in serial +SMP_BATCH_SIZE = 1 + +# Timeout value (optional), default is None, i.e. no timeout limit applied +TIMEOUT = 96:00:00 + + +## Module options + +[MERGE_HEADERS_RUNNER] + +# Input: exp_numbers txt file from find_exposures_runner +INPUT_DIR = last:find_exposures_runner + +FILE_PATTERN = exp_numbers + +FILE_EXT = .txt + +# Tile numbering scheme (RA-Dec, e.g. -301-279) +NUMBERING_SCHEME = -000-000 + +# Root directory containing all per-exposure work directories. +# The runner will walk this tree to collect headers-.npy files. +EXP_BASE_DIR = $SP_EXP + +# If True, log a warning and continue when split_exp_runner output is missing +# for some exposures (e.g. those that failed job 16). If False (default), +# raise an error and stop. +WARN_MISSING_EXP = False diff --git a/example/cfis_image_sims/config_tile_Ng_batch_psfex_sx.ini b/example/cfis_image_sims/config_tile_Ng_batch_psfex_sx.ini new file mode 100644 index 000000000..d29f0ab8f --- /dev/null +++ b/example/cfis_image_sims/config_tile_Ng_batch_psfex_sx.ini @@ -0,0 +1,85 @@ +# ShapePipe configuration file for tiles: ngmix + + +## Default ShapePipe options +[DEFAULT] + +# verbose mode (optional), default: True, print messages on terminal +VERBOSE = True + +# Name of run (optional) default: shapepipe_run +RUN_NAME = run_sp_tile_Ng + +# Add date and time to RUN_NAME, optional, default: False +RUN_DATETIME = True + + +## ShapePipe execution options +[EXECUTION] + +# Module name, single string or comma-separated list of valid module runner names +MODULE = ngmix_runner + +# Parallel processing mode, SMP or MPI +MODE = SMP + + +## ShapePipe file handling options +[FILE] + +# Log file master name, optional, default: shapepipe +LOG_NAME = log_sp + +# Runner log file name, optional, default: shapepipe_runs +RUN_LOG_NAME = log_run_sp + +# Input directory, containing input files, single string or list of names +INPUT_DIR = . + +# Output directory +OUTPUT_DIR = $SP_RUN/output + + +## ShapePipe job handling options +[JOB] + +# Batch size of parallel processing (optional), default is 1, i.e. run all jobs in serial +SMP_BATCH_SIZE = 1 + +# Timeout value (optional), default is None, i.e. no timeout limit applied +TIMEOUT = 96:00:00 + + +## Module options + +# Model-fitting shapes with ngmix +[NGMIX_RUNNER] + +INPUT_DIR = run_sp_tile_Sx:sextractor_runner,run_sp_tile_fpsf:fake_psf_runner,last:vignetmaker_runner_run_2,run_sp_tile_Mh_exp:merge_headers_runner + +FILE_PATTERN = sexcat, image_vignet, galaxy_psf, weight_vignet, flag_vignet, log_exp_headers + +FILE_EXT = .fits, .sqlite, .sqlite, .sqlite, .sqlite, .sqlite + +# NUMBERING_SCHEME (optional) string with numbering pattern for input files +NUMBERING_SCHEME = -000-000 + +# Magnitude zero-point +MAG_ZP = 30.0 + +# No background subtraction for image sims (background not simulated) +BKG_SUB = False + +# Pixel scale in arcsec +PIXEL_SCALE = 0.186 + +SAVE_BATCH = 1000 + +# Position-seeded per-object fixnoise RNG (#796/#803): the SKiLLS shear +# branches share their sky-noise realization, so seeding the metacal fixnoise +# from sky position makes the added noise identical across branches too and +# the Pujol estimator cancels it. +SEED_FROM_POSITION = True + +ID_OBJ_MIN = -1 +ID_OBJ_MAX = -1 diff --git a/example/cfis_image_sims/config_tile_PiViVi_canfar_sx.ini b/example/cfis_image_sims/config_tile_PiViVi_canfar_sx.ini new file mode 100644 index 000000000..025097f81 --- /dev/null +++ b/example/cfis_image_sims/config_tile_PiViVi_canfar_sx.ini @@ -0,0 +1,127 @@ +# ShapePipe configuration file for tile, from detection up to shape measurement. +# Fake PSF model. + + +## Default ShapePipe options +[DEFAULT] + +# verbose mode (optional), default: True, print messages on terminal +VERBOSE = True + +# Name of run (optional) default: shapepipe_run +RUN_NAME = run_sp_tile_ViVi + +# Add date and time to RUN_NAME, optional, default: False +; RUN_DATETIME = False + + +## ShapePipe execution options +[EXECUTION] + +# Module name, single string or comma-separated list of valid module runner names + +MODULE = vignetmaker_runner, vignetmaker_runner + +# Parallel processing mode, SMP or MPI +MODE = SMP + + +## ShapePipe file handling options +[FILE] + +# Log file master name, optional, default: shapepipe +LOG_NAME = log_sp + +# Runner log file name, optional, default: shapepipe_runs +RUN_LOG_NAME = log_run_sp + +# Input directory, containing input files, single string or list of names +INPUT_DIR = . + +# Output directory +OUTPUT_DIR = $SP_RUN/output + + +## ShapePipe job handling options +[JOB] + +# Batch size of parallel processing (optional), default is 1, i.e. run all jobs in serial +SMP_BATCH_SIZE = 1 + +# Timeout value (optional), default is None, i.e. no timeout limit applied +TIMEOUT = 96:00:00 + + +## Module options + +# Create vignets for tiles weights +[VIGNETMAKER_RUNNER_RUN_1] + +INPUT_DIR = run_sp_tile_Sx:sextractor_runner, last:uncompress_fits_runner + +FILE_PATTERN = sexcat, CFIS_simu_weight + +FILE_EXT = .fits, .fits + +# NUMBERING_SCHEME (optional) string with numbering pattern for input files +NUMBERING_SCHEME = -000-000 + +MASKING = False +MASK_VALUE = 0 + +# Run mode for psfex interpolation: +# CLASSIC: 'classical' run, interpolate to object positions +# MULTI-EPOCH: interpolate for multi-epoch images +# VALIDATION: validation for single-epoch images +MODE = CLASSIC + +# Coordinate frame type, one in PIX (pixel frame), SPHE (spherical coordinates) +COORD = PIX +POSITION_PARAMS = XWIN_IMAGE,YWIN_IMAGE + +# Vignet size in pixels +STAMP_SIZE = 51 + +# Output file name prefix, file name is _vignet.fits +PREFIX = weight + + +[VIGNETMAKER_RUNNER_RUN_2] + +# Create multi-epoch vignets for tiles corresponding to +# positions on single-exposures + +INPUT_DIR = run_sp_tile_Sx:sextractor_runner, run_sp_tile_Mh_exp:merge_headers_runner, last:find_exposures_runner + +FILE_PATTERN = sexcat, log_exp_headers, exp_numbers + +FILE_EXT = .fits, .sqlite, .txt + +# NUMBERING_SCHEME (optional) string with numbering pattern for input files +NUMBERING_SCHEME = -000-000 + +MASKING = False +MASK_VALUE = 0 + +# Run mode for psfex interpolation: +# CLASSIC: 'classical' run, interpolate to object positions +# MULTI-EPOCH: interpolate for multi-epoch images +# VALIDATION: validation for single-epoch images +MODE = MULTI-EPOCH + +# Coordinate frame type, one in PIX (pixel frame), SPHE (spherical coordinates) +COORD = SPHE +POSITION_PARAMS = XWIN_WORLD,YWIN_WORLD + +# Vignet size in pixels +STAMP_SIZE = 51 + +# Output file name prefix, file name is vignet.fits +PREFIX = + +# Additional parameters for path and file pattern corresponding to single-exposure +# run outputs. ME_IMAGE_EXP_DIR/ME_IMAGE_EXP_RUNNERS replace ME_IMAGE_DIR for +# the v2.0 per-exposure pipeline; output dirs are discovered by scanning $SP_EXP. +ME_IMAGE_EXP_DIR = $SP_EXP +ME_IMAGE_EXP_RUNNERS = split_exp_runner, split_exp_runner, split_exp_runner +ME_IMAGE_PATTERN = flag, image, weight diff --git a/example/cfis_image_sims/config_tile_Sx_nomask.ini b/example/cfis_image_sims/config_tile_Sx_nomask.ini new file mode 100644 index 000000000..810614c5d --- /dev/null +++ b/example/cfis_image_sims/config_tile_Sx_nomask.ini @@ -0,0 +1,114 @@ +# ShapePipe configuration file for tile detection + + +## Default ShapePipe options +[DEFAULT] + +# verbose mode (optional), default: True, print messages on terminal +VERBOSE = True + +# Name of run (optional) default: shapepipe_run +RUN_NAME = run_sp_tile_Sx + +# Add date and time to RUN_NAME, optional, default: True +; RUN_DATETIME = False + + +## ShapePipe execution options +[EXECUTION] + +# Module name, single string or comma-separated list of valid module runner names +MODULE = sextractor_runner + + +# Run mode, SMP or MPI +MODE = SMP + + +## ShapePipe file handling options +[FILE] + +# Log file master name, optional, default: shapepipe +LOG_NAME = log_sp + +# Runner log file name, optional, default: shapepipe_runs +RUN_LOG_NAME = log_run_sp + +# Input directory, containing input files, single string or list of names with length matching FILE_PATTERN +INPUT_DIR = $SP_RUN/output + +# Output directory +OUTPUT_DIR = $SP_RUN/output + + +## ShapePipe job handling options +[JOB] + +# Batch size of parallel processing (optional), default is 1, i.e. run all jobs in serial +SMP_BATCH_SIZE = 1 + +# Timeout value (optional), default is None, i.e. no timeout limit applied +TIMEOUT = 96:00:00 + + +## Module options + +[SEXTRACTOR_RUNNER] + +INPUT_DIR = run_sp_tile_Git:get_images_runner, last:uncompress_fits_runner, run_sp_tile_Mh_exp:merge_headers_runner + +FILE_PATTERN = CFIS_simu_image, CFIS_simu_weight, log_exp_headers + +FILE_EXT = .fits, .fits, .sqlite + +# NUMBERING_SCHEME (optional) string with numbering pattern for input files +NUMBERING_SCHEME = -000-000 + +# SExtractor executable path +EXEC_PATH = source-extractor + +# SExtractor configuration files +DOT_SEX_FILE = $SP_CONFIG/default_tile.sex +DOT_PARAM_FILE = $SP_CONFIG/default_noimaflags.param +DOT_CONV_FILE = $SP_CONFIG/default.conv + +# Use input weight image if True +WEIGHT_IMAGE = True + +# Use input flag image if True +FLAG_IMAGE = False + +# Use input PSF file if True +PSF_FILE = False + +# Use distinct image for detection (SExtractor in +# dual-image mode) if True +DETECTION_IMAGE = False + +# Distinct weight image for detection (SExtractor +# in dual-image mode) +DETECTION_WEIGHT = False + +ZP_FROM_HEADER = False + +BKG_FROM_HEADER = False + +# Type of image check (optional), default not used, can be a list of +# BACKGROUND, BACKGROUND_RMS, INIBACKGROUND, +# MINIBACK_RMS, -BACKGROUND, #FILTERED, +# OBJECTS, -OBJECTS, SEGMENTATION, APERTURES +CHECKIMAGE = BACKGROUND + +# File name suffix for the output sextractor files (optional) +SUFFIX = sexcat + +## Post-processing + +# Necessary for tiles, to enable multi-exposure processing +MAKE_POST_PROCESS = True + +# World coordinate keywords, SExtractor output. Format: KEY_X,KEY_Y +WORLD_POSITION = XWIN_WORLD,YWIN_WORLD + +# Number of pixels in x,y of a CCD. Format: Nx,Ny +CCD_SIZE = 33,2080,1,4612 diff --git a/example/cfis_image_sims/config_tile_onthefly.mask_simu b/example/cfis_image_sims/config_tile_onthefly.mask_simu new file mode 100644 index 000000000..f42d05d39 --- /dev/null +++ b/example/cfis_image_sims/config_tile_onthefly.mask_simu @@ -0,0 +1,90 @@ +# Mask module config file for tiles + +## Paths to executables +[PROGRAM_PATH] + +WW_PATH = weightwatcher +WW_CONFIG_FILE = $SP_CONFIG/mask_default/default.ww + +# Indicate cds client executable if no external star catalogue is available +# (e.g. no internet access on run nodes) +CDSCLIENT_PATH = findgsc2.2 + +## Border parameters +[BORDER_PARAMETERS] + +BORDER_MAKE = False + +BORDER_WIDTH = 1 +BORDER_FLAG_VALUE = 4 + + +## Halo parameters +[HALO_PARAMETERS] + +HALO_MAKE = False + +HALO_MASKMODEL_PATH = $SP_CONFIG/mask_default/halo_mask.reg +HALO_MAG_LIM = 13. +HALO_SCALE_FACTOR = 0.05 +HALO_MAG_PIVOT = 13.8 +HALO_FLAG_VALUE = 2 +HALO_REG_FILE = halo.reg + + +## Diffraction pike parameters +[SPIKE_PARAMETERS] + +SPIKE_MAKE = False + +SPIKE_MASKMODEL_PATH = $SP_CONFIG/mask_default/MEGAPRIME_star_i_13.8.reg +SPIKE_MAG_LIM = 18. +SPIKE_SCALE_FACTOR = 0.3 +SPIKE_MAG_PIVOT = 13.8 +SPIKE_FLAG_VALUE = 128 +SPIKE_REG_FILE = spike.reg + + +## Messier parameters +[MESSIER_PARAMETERS] + +MESSIER_MAKE = False + +MESSIER_CAT_PATH = $SP_CONFIG/mask_default/Messier_catalog_updated.fits +MESSIER_PIXEL_SCALE = 0.187 +MESSIER_SIZE_PLUS = 0. +MESSIER_FLAG_VALUE = 16 + +## NGC mask +[NGC_PARAMETERS] + +NGC_MAKE = False + +NGC_CAT_PATH = $SP_CONFIG/mask_default/ngc_cat.fits +NGC_SIZE_PLUS = 0. +NGC_FLAG_VALUE = 32 + + +## External flag +[EXTERNAL_FLAG] + +EF_MAKE = False + + +## Missing data parameters +[MD_PARAMETERS] + +MD_MAKE = False + +MD_THRESH_FLAG = 0.3 +MD_THRESH_REMOVE = 0.75 +MD_REMOVE = False + + +## Other parameters +[OTHER] + +KEEP_REG_FILE = False +KEEP_INDIVIDUAL_MASK = False + +TEMP_DIRECTORY = .temp_tiles diff --git a/example/cfis_image_sims/default.conv b/example/cfis_image_sims/default.conv new file mode 120000 index 000000000..bd71df850 --- /dev/null +++ b/example/cfis_image_sims/default.conv @@ -0,0 +1 @@ +../cfis/default.conv \ No newline at end of file diff --git a/example/cfis_image_sims/default.param b/example/cfis_image_sims/default.param new file mode 120000 index 000000000..49e000314 --- /dev/null +++ b/example/cfis_image_sims/default.param @@ -0,0 +1 @@ +../cfis/default.param \ No newline at end of file diff --git a/example/cfis_image_sims/default_noimaflags.param b/example/cfis_image_sims/default_noimaflags.param new file mode 120000 index 000000000..75451801e --- /dev/null +++ b/example/cfis_image_sims/default_noimaflags.param @@ -0,0 +1 @@ +../cfis/default_noimaflags.param \ No newline at end of file diff --git a/example/cfis_image_sims/default_tile.sex b/example/cfis_image_sims/default_tile.sex new file mode 120000 index 000000000..8770da87e --- /dev/null +++ b/example/cfis_image_sims/default_tile.sex @@ -0,0 +1 @@ +../cfis/default_tile.sex \ No newline at end of file diff --git a/example/cfis_image_sims/final_cat.param b/example/cfis_image_sims/final_cat.param new file mode 100644 index 000000000..5f96d959e --- /dev/null +++ b/example/cfis_image_sims/final_cat.param @@ -0,0 +1,111 @@ +# Final-catalogue column selection for the image-simulation pipeline. +# +# ShapePipe-v2 (post-#761) ngmix column grammar: ellipticity in named scalar +# components NGMIX_G{1,2}_*, PSF size split into original (T_PSF_ORIG) and +# reconvolved (T_PSF_RECONV) measurements. create_final_cat.py selects exactly +# these columns from the per-tile make_cat output into final_cat_{sim}.hdf5; +# every entry must exist in that catalogue. +# +# This is the image-sims sibling of example/cfis/final_cat.param. The two +# differ only in columns the simulation pipeline does not produce: +# - IMAFLAGS_ISO : the sims skip the imaging-flag masking stage +# - NGMIX_MOM_FAIL : the v2 sims path emits NGMIX_MCAL_TYPES_FAIL instead, +# which the downstream calibration does not read + +# coordinates +XWIN_WORLD +YWIN_WORLD + +# tile ID, for plot of tile-dependent additive bias. +TILE_ID + +# SExtractor number +NUMBER + +# flags +FLAGS +NGMIX_MCAL_FLAGS +NGMIX_MCAL_TYPES_FAIL + +# PSF ellipticity (original image PSF) +NGMIX_G1_PSF_ORIG_NOSHEAR +NGMIX_G2_PSF_ORIG_NOSHEAR + +# Number of epochs (exposures) +N_EPOCH +NGMIX_N_EPOCH + +## Shape measurement outputs +## Ngmix: model fitting + +# galaxy ellipticity +NGMIX_G1_1M +NGMIX_G2_1M +NGMIX_G1_1P +NGMIX_G2_1P +NGMIX_G1_2M +NGMIX_G2_2M +NGMIX_G1_2P +NGMIX_G2_2P +NGMIX_G1_NOSHEAR +NGMIX_G2_NOSHEAR +NGMIX_G1_ERR_NOSHEAR +NGMIX_G2_ERR_NOSHEAR + +# flags +NGMIX_FLAGS_1M +NGMIX_FLAGS_1P +NGMIX_FLAGS_2M +NGMIX_FLAGS_2P +NGMIX_FLAGS_NOSHEAR + +# size and error +NGMIX_T_1M +NGMIX_T_1P +NGMIX_T_2M +NGMIX_T_2P +NGMIX_T_NOSHEAR +NGMIX_T_ERR_1M +NGMIX_T_ERR_1P +NGMIX_T_ERR_2M +NGMIX_T_ERR_2P +NGMIX_T_ERR_NOSHEAR + +# reconvolved PSF size (used by the metacal size cut / Tpsf) +NGMIX_T_PSF_RECONV_1M +NGMIX_T_PSF_RECONV_1P +NGMIX_T_PSF_RECONV_2M +NGMIX_T_PSF_RECONV_2P +NGMIX_T_PSF_RECONV_NOSHEAR + +# flux and error +NGMIX_FLUX_1M +NGMIX_FLUX_1P +NGMIX_FLUX_2M +NGMIX_FLUX_2P +NGMIX_FLUX_NOSHEAR +NGMIX_FLUX_ERR_1M +NGMIX_FLUX_ERR_1P +NGMIX_FLUX_ERR_2M +NGMIX_FLUX_ERR_2P +NGMIX_FLUX_ERR_NOSHEAR + +# magnitudes +MAG_AUTO +MAGERR_AUTO +MAG_WIN +MAGERR_WIN +FLUX_AUTO +FLUXERR_AUTO +FLUX_APER +FLUXERR_APER +FLUX_RADIUS + +# SNR from SExtractor +SNR_WIN + +FWHM_IMAGE +FWHM_WORLD + +# PSF size measured on original image +NGMIX_T_PSF_ORIG_NOSHEAR diff --git a/example/cfis_image_sims/input_dirs.conf b/example/cfis_image_sims/input_dirs.conf new file mode 100644 index 000000000..cd3661983 --- /dev/null +++ b/example/cfis_image_sims/input_dirs.conf @@ -0,0 +1,12 @@ +# Configuration template for init_run_v2.0.sh +# Specifies input directories for image simulations +# Copy and edit this file, then pass with: -c /path/to/config + +# Base directory for input data +input_dir_base="/n09data/hervas/skills_out" + +# Tile images subdirectory +input_dir_tiles="${input_dir_base}/${subdir}/images/SP_tiles" + +# Exposure images subdirectory +input_dir_exp="${input_dir_base}/${subdir}/images/SP_exp" diff --git a/scripts/image_sims_pipeline/Snakefile b/scripts/image_sims_pipeline/Snakefile new file mode 100644 index 000000000..ab04ff6e8 --- /dev/null +++ b/scripts/image_sims_pipeline/Snakefile @@ -0,0 +1,443 @@ +"""Snakemake workflow for image-simulation pipeline. + +Run from the run directory: + snakemake -s Snakefile --configfile config.yaml -j 1 -- coverage # once per grid, runs init if needed + snakemake -s Snakefile --configfile config.yaml -j 4 --rerun-triggers=mtime -- pipeline_all + snakemake -s Snakefile --configfile config.yaml -j 1 --dry-run -- pipeline_all + +ShapePipe stages run in one-shot mode by default: one job per (sim, tile) +covering the whole bit sequence (-j 4095); n_sims x n_tiles independent +jobs. Use --config per_bit=1 (implied by job=N) for per-bit jobs. + +Cluster execution (SLURM, snakemake >= 8 with snakemake-executor-plugin-slurm): +run from a login-node tmux; each job becomes its own sbatch submission: + snakemake -s Snakefile --configfile config.yaml --executor slurm -j 80 \ + --default-resources slurm_partition= cpus_per_task=2 \ + mem_mb=8000 runtime=1440 --rerun-triggers=mtime -- pipeline_all + +When running locally, use a small fixed -j (e.g. 4), not -j +$OMP_NUM_THREADS: each job spawns MPI processes, and ulimit -u on candide +is 1200 ("can't start new thread" crash). + +--rerun-triggers=mtime avoids re-running completed jobs after edits to this +file or to config.yaml (default triggers include rule-code changes). +""" + +import os + +# ---------- configuration ---------- +configfile: "config.yaml" + +BASE = config["base"] +SAMPLE = str(config["sample"]) +NUM = config["num"] +TYPE = config["sims_type"] + +# Load tile IDs from file (always a file path) +# If tile_IDs not specified, construct path with num-specific file +if "tile_IDs" in config: + _tile_ids_file = config["tile_IDs"].replace("{num}", str(NUM)) +else: + # Try num-specific tile file first + _tile_ids_file = f"{os.path.expanduser('~')}/shapepipe/auxdir/CFIS/im_sims_202606/tile_numbers_{NUM}.txt" + if not os.path.exists(_tile_ids_file): + # Fall back to default + _tile_ids_file = f"{os.path.expanduser('~')}/shapepipe/auxdir/CFIS/im_sims_202606/tile_numbers.txt" + +with open(_tile_ids_file) as f: + TILE_IDs = [line.strip() for line in f if line.strip()] + +# Apply exclusive tile list if provided (for testing/limiting tiles) +_exclusive = config.get("tile_IDs_exclusive") +if _exclusive and isinstance(_exclusive, list) and _exclusive: # Non-empty list + TILE_IDs = _exclusive + +# Remove excluded tiles if provided (e.g. blank input images) +_exclude = config.get("tile_IDs_exclude") +if _exclude and isinstance(_exclude, list): + _exclude_set = set(str(t) for t in _exclude) + TILE_IDs = [t for t in TILE_IDs if t not in _exclude_set] + +# Remove low-coverage tiles found by check_tile_coverage.py, if the scan +# has been run ("snakemake coverage"). Read at DAG-build time and kept out +# of rule inputs so a rescan does not invalidate completed jobs. +COVERAGE_FILE = f"{BASE}/tile_coverage_grid_{config['num']}.yaml" +if os.path.isfile(COVERAGE_FILE): + import yaml as _yaml + with open(COVERAGE_FILE) as _f: + _cov_exclude = set(str(t) for t in _yaml.safe_load(_f).get("exclude", [])) + _n_before = len(TILE_IDs) + TILE_IDs = [t for t in TILE_IDs if t not in _cov_exclude] + if len(TILE_IDs) < _n_before: + print( + f"Coverage scan: excluding {_n_before - len(TILE_IDs)} " + f"low-coverage tiles ({len(TILE_IDs)} remain), see {COVERAGE_FILE}" + ) +else: + print( + f"No coverage scan found ({COVERAGE_FILE}); " + f"run 'snakemake coverage' first to exclude (near-)empty input tiles" + ) +FORCE = "--force" if config.get("force", 0) else "" +N_SMP = config.get("n_smp", -1) +SHEAR_AMP = config.get("shear_amplitude", 0.02) +MATCH_RAD = config.get("match_radius_deg", 0.0002) +PAIR_MATCH = "true" if config.get("pair_match", True) else "false" +# Weight schemes for m-bias; "none" means uniform weights, first is primary +_w_cols = config.get("w_cols") or ( + [config["w_col"]] if config.get("w_col") else ["none", "w_iv"] +) +W_COLS = "[" + ", ".join(f"'{w}'" for w in _w_cols) + "]" +N_BOOT = config.get("n_bootstrap", 500) + +FULL_JOB_SEQUENCE = [1, 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024, 2048] +# Filter by bitmask if specified via --config "job=N" (N may combine multiple bits) +if config.get("job"): + target_job = int(config["job"]) + JOB_SEQUENCE = [j for j in FULL_JOB_SEQUENCE if j & target_job] +else: + JOB_SEQUENCE = list(FULL_JOB_SEQUENCE) +JOB = sum(JOB_SEQUENCE) + +PSF = "psfex" + +STR_TYPE = f"_{TYPE}" if TYPE == "grid" else "" +SIMS = [ + f"1{dg1}2{dg2}{STR_TYPE}_{NUM}" + for dg1, dg2 in zip(["m", "p", "z", "z"], + ["z", "z", "m", "p"]) +] + +HOME = os.path.expanduser("~") +PATH_SH = f"{HOME}/astro/repositories/github/shapepipe/scripts/sh" +PATH_PY = f"{HOME}/astro/repositories/github/shapepipe/scripts/python" +PATH_IMSIM = f"{HOME}/astro/repositories/github/shapepipe/scripts/image_sims_pipeline" +PATH_NB = f"{HOME}/astro/repositories/github/sp_validation/notebooks" +PATH_SPV = f"{HOME}/astro/repositories/github/sp_validation" + +SPV_SIF = "/n17data/mkilbing/sp_validation_im_sims.sif" +# The same container holds shapepipe and sp_validation code +SP_SIF = SPV_SIF +SPV_SRC = f"{HOME}/astro/repositories/github/sp_validation/src" +SP_SRC = f"{HOME}/astro/repositories/github/shapepipe/src" + +# Cap threading libraries inside containers; the host OMP_NUM_THREADS would +# otherwise leak in and N parallel jobs x N threads exhausts ulimit -u +_THREAD_CAPS = ( + "OMP_NUM_THREADS=1,OPENBLAS_NUM_THREADS=1," + "MKL_NUM_THREADS=1,NUMEXPR_NUM_THREADS=1" +) + +# ShapePipe pipeline stages (uses shapepipe_im_sims.sif — correct MPI build) +_APPTAINER = f"{PATH_SH}/apptainer_noslurm.sh" +APP_CMD = ( + f"{_APPTAINER} --bind /n17data,/n09data,/home" + f" --env PYTHONPATH={SP_SRC}" + f" --env PYTHONWARNINGS=ignore::SyntaxWarning" + f" --env {_THREAD_CAPS}" + f" {SP_SIF}" +) + +# sp_validation stages (uses sp_validation_im_sims.sif with local src override) +SPV_CMD = ( + f"{_APPTAINER} --bind /n17data,/n09data,/home" + f" --env PYTHONPATH={SPV_SRC}" + f" --env {_THREAD_CAPS}" + f" {SPV_SIF}" +) +MBIAS_OUT = f"{BASE}/results/m_bias_results.yaml" + +# ---------- helper functions ---------- +def validate_output_files(output_files): + """Validate that output files exist and are non-empty. + + Raises an error if any output file is missing or zero-sized, + preventing empty/corrupt files from being treated as valid outputs. + """ + if isinstance(output_files, str): + output_files = [output_files] + + for f in output_files: + if not os.path.isfile(f): + raise FileNotFoundError(f"Output file missing: {f}") + size = os.path.getsize(f) + if size == 0: + raise ValueError(f"Output file is empty (0 bytes): {f}") + +# ---------- targets ---------- +rule all: + input: + expand(f"{BASE}/{{sim}}/shape_catalog_cut_ngmix.fits", sim=SIMS) + +# Individual job bit targets (run without config filtering) +rule job_1_only: + input: + expand(f"{BASE}/{{sim}}/logs/log_job_{{tile}}_1.txt", sim=SIMS, tile=TILE_IDs) + +rule job_2_only: + input: + expand(f"{BASE}/{{sim}}/logs/log_job_{{tile}}_2.txt", sim=SIMS, tile=TILE_IDs) + +rule job_4_only: + input: + expand(f"{BASE}/{{sim}}/logs/log_job_{{tile}}_4.txt", sim=SIMS, tile=TILE_IDs) + +rule job_8_only: + input: + expand(f"{BASE}/{{sim}}/logs/log_job_{{tile}}_8.txt", sim=SIMS, tile=TILE_IDs) + +rule job_16_only: + input: + expand(f"{BASE}/{{sim}}/logs/log_job_{{tile}}_16.txt", sim=SIMS, tile=TILE_IDs) + +rule job_32_only: + input: + expand(f"{BASE}/{{sim}}/logs/log_job_{{tile}}_32.txt", sim=SIMS, tile=TILE_IDs) + +rule job_64_only: + input: + expand(f"{BASE}/{{sim}}/logs/log_job_{{tile}}_64.txt", sim=SIMS, tile=TILE_IDs) + +rule job_128_only: + input: + expand(f"{BASE}/{{sim}}/logs/log_job_{{tile}}_128.txt", sim=SIMS, tile=TILE_IDs) + +rule job_256_only: + input: + expand(f"{BASE}/{{sim}}/logs/log_job_{{tile}}_256.txt", sim=SIMS, tile=TILE_IDs) + +rule job_512_only: + input: + expand(f"{BASE}/{{sim}}/logs/log_job_{{tile}}_512.txt", sim=SIMS, tile=TILE_IDs) + +rule job_1024_only: + input: + expand(f"{BASE}/{{sim}}/logs/log_job_{{tile}}_1024.txt", sim=SIMS, tile=TILE_IDs) + +rule job_2048_only: + input: + expand(f"{BASE}/{{sim}}/logs/log_job_{{tile}}_2048.txt", sim=SIMS, tile=TILE_IDs) + +# ---------- convenience targets (run in order) ---------- +# 1. init_all +# 2. pipeline_all +# 3. merge_all +# 4. extract_all +# 5. calibrate_all +# 6. mbias +# 7. diagnostics + +# Pixel thinning for the coverage scan: read every Nth pixel per axis +COVERAGE_SUBSAMPLE = 20 +# Weight-map location is read from this variant's get-images ini (from init) +SIM_REF = f"1p2z{STR_TYPE}_{NUM}" + +rule coverage: + """Scan input weight maps; write low-coverage tile exclude list. + + Run this once per grid, after init_all and BEFORE pipeline_all: the + exclude list is read at DAG-build time, so it only takes effect on + the next invocation. Rerun with -f after the input sims change + (e.g. tile regeneration). + """ + input: + f"{BASE}/{SIM_REF}/params.py" + output: + COVERAGE_FILE + shell: + f"{SPV_CMD} python {PATH_IMSIM}/check_tile_coverage.py" + f" -c {BASE}/config.yaml -o {{output}}" + f" --subsample {COVERAGE_SUBSAMPLE} -v" + +rule generate_tile_list: + """Generate tile_numbers.txt from config's tile_IDs.""" + output: + f"{BASE}/tile_numbers.txt" + run: + os.makedirs(f"{BASE}", exist_ok=True) + with open(output[0], "w") as f: + for tile in TILE_IDs: + f.write(f"{tile}\n") + print(f"Generated {output[0]} with {len(TILE_IDs)} tiles") + +rule init_all: + input: + expand(f"{BASE}/{{sim}}/params.py", sim=SIMS) + +rule pipeline_all: + input: + expand( + f"{BASE}/{{sim}}/logs/log_job_{{tile}}_{JOB_SEQUENCE[-1]}.txt", + sim=SIMS, tile=TILE_IDs + ) + +rule merge_all: + input: + expand(f"{BASE}/{{sim}}/final_cat_{{sim}}.hdf5", sim=SIMS) + +rule extract_all: + input: + expand(f"{BASE}/{{sim}}/shape_catalog_comprehensive_ngmix.hdf5", sim=SIMS) + +rule calibrate_all: + input: + expand(f"{BASE}/{{sim}}/shape_catalog_cut_ngmix.fits", sim=SIMS) + +# ---------- rules ---------- +# config.yaml is deliberately not a tracked input: editing it (e.g. to +# exclude tiles) must not invalidate init and cascade into a full re-run. +# After config changes that require re-init, run: snakemake -f init_all +rule init: + output: + params = f"{BASE}/{{sim}}/params.py", + mask = f"{BASE}/{{sim}}/config_mask.yaml", + tiles = f"{BASE}/{{sim}}/tile_numbers.txt", + cfis = directory(f"{BASE}/{{sim}}/cfis"), + shell: + f"python3 {PATH_PY}/init_run_v2.0.py" + f" -c {BASE}/config.yaml -s {{wildcards.sim}} -q" + +# Two execution modes for the ShapePipe stages: +# - one-shot (default): one job per (sim, tile) running the whole (filtered) +# bit sequence via a single -j {JOB} call. Ideal for cluster execution: +# n_sims x n_tiles independent jobs. Resume-safe: the run script skips +# bits whose existing run is complete. +# - per-bit (--config per_bit=1, or implied by --config job=N): one job per +# (sim, tile, bit) with inter-bit dependencies; for surgical reruns. +PER_BIT = bool(config.get("job")) or bool(config.get("per_bit", 0)) + +if PER_BIT: + def _shapepipe_job_input(wildcards): + bit = int(wildcards.bit) + tile = wildcards.tile + bit_idx = JOB_SEQUENCE.index(bit) + + if bit_idx > 0: + # Previous bit in the filtered sequence + return f"{BASE}/{wildcards.sim}/logs/log_job_{tile}_{JOB_SEQUENCE[bit_idx - 1]}.txt" + + # First bit in the filtered sequence: find predecessor in the full sequence + full_idx = FULL_JOB_SEQUENCE.index(bit) + if full_idx == 0: + return [ + f"{BASE}/{wildcards.sim}/params.py", + f"{BASE}/{wildcards.sim}/config_mask.yaml", + ] + return f"{BASE}/{wildcards.sim}/logs/log_job_{tile}_{FULL_JOB_SEQUENCE[full_idx - 1]}.txt" + + rule shapepipe_job: + input: _shapepipe_job_input + output: f"{BASE}/{{sim}}/logs/log_job_{{tile}}_{{bit}}.txt" + # log file keeps stdout/stderr; unlike output, it survives job failure + log: f"{BASE}/{{sim}}/logs/err_job_{{tile}}_{{bit}}.txt" + wildcard_constraints: + bit = "|".join(str(b) for b in JOB_SEQUENCE), + tile = "|".join(t.replace(".", r"\.") for t in TILE_IDs), + shell: + f"( cd {BASE}/{{wildcards.sim}} &&" + f" {APP_CMD} {PATH_SH}/run_job_sp_canfar_v2.0.bash" + f" {FORCE} -N {N_SMP} -e {{wildcards.tile}} -t image_sims -j {{wildcards.bit}} &&" + f" python3 -c \"import os; f='{{output}}'; " + f"os.path.getsize(f) > 0 or exit(1)\"" + f" ) 2>&1 | tee {{log}}" +else: + rule shapepipe_tile: + input: + f"{BASE}/{{sim}}/params.py", + f"{BASE}/{{sim}}/config_mask.yaml", + # Last bit's log marks the whole chain done; per-bit logs are still + # written along the way (info.py monitoring keeps working) + output: + f"{BASE}/{{sim}}/logs/log_job_{{tile}}_{JOB_SEQUENCE[-1]}.txt" + log: + f"{BASE}/{{sim}}/logs/err_job_{{tile}}_all.txt" + wildcard_constraints: + tile = "|".join(t.replace(".", r"\.") for t in TILE_IDs), + shell: + f"( cd {BASE}/{{wildcards.sim}} &&" + f" {APP_CMD} {PATH_SH}/run_job_sp_canfar_v2.0.bash" + f" {FORCE} -N {N_SMP} -e {{wildcards.tile}} -t image_sims -j {JOB} &&" + f" python3 -c \"import os; f='{{output}}'; " + f"os.path.getsize(f) > 0 or exit(1)\"" + f" ) 2>&1 | tee {{log}}" + +rule merge: + input: + expand( + f"{BASE}/{{{{sim}}}}/logs/log_job_{{tile}}_{JOB_SEQUENCE[-1]}.txt", + tile=TILE_IDs + ) + output: + f"{BASE}/{{sim}}/final_cat_{{sim}}.hdf5" + shell: + f"cd {BASE}/{{wildcards.sim}} &&" + f" {SPV_CMD} python {PATH_PY}/create_final_cat.py" + f" -I -m final_cat_{{wildcards.sim}}.hdf5" + f" -i .. -p cfis/final_cat.param -P {{wildcards.sim}}" + f" -o n_tiles_final.txt -v" + +rule extract: + input: + cat = f"{BASE}/{{sim}}/final_cat_{{sim}}.hdf5", + params = f"{BASE}/{{sim}}/params.py", + output: + f"{BASE}/{{sim}}/shape_catalog_comprehensive_ngmix.hdf5" + shell: + f"cd {BASE}/{{wildcards.sim}} &&" + f" {SPV_CMD} python {PATH_NB}/extract_info.py" + +rule calibrate: + input: + cat = f"{BASE}/{{sim}}/shape_catalog_comprehensive_ngmix.hdf5", + mask = f"{BASE}/{{sim}}/config_mask.yaml", + output: + f"{BASE}/{{sim}}/shape_catalog_cut_ngmix.fits" + shell: + f"cd {BASE}/{{wildcards.sim}} &&" + f" {SPV_CMD} python {PATH_SPV}/notebooks/calibrate_comprehensive_cat.py -s calibrate" + +rule mbias: + input: + expand(f"{BASE}/{{sim}}/shape_catalog_cut_ngmix.fits", sim=SIMS) + # mbias_cumulative.yaml is deliberately not declared as an output: + # it accumulates results across runs, and Snakemake deletes declared + # outputs before rerunning a rule + output: + MBIAS_OUT, + shell: + f"mkdir -p {BASE}/results && cat > /tmp/mbias_cfg.yaml << 'EOF'\n" + f"base: {BASE}\n" + f"num: {NUM}\n" + f"catalog_name: shape_catalog_cut_ngmix.fits\n" + f"shear_amplitude: {SHEAR_AMP}\n" + f"match_radius_deg: {MATCH_RAD}\n" + f"pair_match: {PAIR_MATCH}\n" + f"w_cols: {W_COLS}\n" + f"n_bootstrap: {N_BOOT}\n" + f"output_path: {BASE}/results/m_bias_results.yaml\n" + f"results_dir: {BASE}/results\n" + f"EOF\n" + f"{SPV_CMD} python {PATH_SPV}/scripts/compute_m_bias_image_sims.py" + f" -c /tmp/mbias_cfg.yaml -v --cumulative" + +rule diagnostics: + input: + expand(f"{BASE}/{{sim}}/shape_catalog_cut_ngmix.fits", sim=SIMS) + output: + f"{BASE}/results/footprint.png", + shell: + f"mkdir -p {BASE}/results && cat > /tmp/diag_cfg.yaml << 'EOF'\n" + f"base: {BASE}\n" + f"num: {NUM}\n" + f"catalog_name: shape_catalog_cut_ngmix.fits\n" + f"shear_amplitude: {SHEAR_AMP}\n" + f"match_radius_deg: {MATCH_RAD}\n" + f"pair_match: {PAIR_MATCH}\n" + f"w_cols: {W_COLS}\n" + f"n_bootstrap: {N_BOOT}\n" + f"output_path: {BASE}/results/m_bias_results.yaml\n" + f"results_dir: {BASE}/results\n" + f"diagnostics_dir: {BASE}/results\n" + f"EOF\n" + f"{SPV_CMD} python {PATH_SPV}/scripts/diagnostics_image_sims.py" + f" -c /tmp/diag_cfg.yaml -v" diff --git a/scripts/image_sims_pipeline/check_tile_coverage.py b/scripts/image_sims_pipeline/check_tile_coverage.py new file mode 100644 index 000000000..b36684d14 --- /dev/null +++ b/scripts/image_sims_pipeline/check_tile_coverage.py @@ -0,0 +1,166 @@ +#!/usr/bin/env python3 +"""Check input-tile coverage for image simulations. + +Scans the weight maps of the input simulation tiles and computes the +fraction of nonzero pixels (subsampled). Tiles below the coverage +threshold are written to an exclude list, which the Snakefile reads at +DAG-build time to drop them from the run. + +The weight-map location is read from the ShapePipe get-images config +{base}/{sim}/cfis/config_tile_Git_symlink.ini (created by init), for the +first sheared variant (1p2z); coverage is identical across variants. + +(Near-)empty tiles occur at the edge of the simulated footprint; they +crash the pipeline downstream (vignetmaker on 0 objects, ngmix on empty +PSF vignettes) and their edge objects bias the m-bias estimate. + +Usage: + check_tile_coverage.py -c config.yaml [-o output.yaml] [-t THRESH] + [-s SUBSAMPLE] [-v] +""" + +import argparse +import configparser +import os +import re +import sys + +import numpy as np +import yaml +from astropy.io import fits + + +def parse_args(): + p = argparse.ArgumentParser(description=__doc__) + p.add_argument("-c", "--config", required=True, help="run config YAML file") + p.add_argument( + "-o", "--output", + help="output YAML path (default {base}/tile_coverage_grid_{num}.yaml)", + ) + p.add_argument( + "-t", "--threshold", type=float, + help="minimum coverage fraction (default: min_tile_coverage from" + " config, or 0.5)", + ) + p.add_argument( + "-s", "--subsample", type=int, default=20, + help="read every Nth pixel per axis (default: %(default)s)", + ) + p.add_argument("-v", "--verbose", action="store_true", help="verbose output") + return p.parse_args() + + +def get_tile_ids(config): + """Read the tile ID list the same way the Snakefile does.""" + num = config["num"] + path = config["tile_IDs"].replace("{num}", str(num)) + with open(path) as f: + tile_ids = [line.strip() for line in f if line.strip()] + + exclusive = config.get("tile_IDs_exclusive") + if exclusive and isinstance(exclusive, list): + tile_ids = exclusive + + return tile_ids + + +def get_weight_template(config): + """Weight-map path template from the get-images ShapePipe config. + + Reads {base}/{sim}/cfis/config_tile_Git_symlink.ini of the first + sheared variant (1p2z) and returns (directory, pattern, extension, + numbering regex) of the weight entry. + """ + num = config["num"] + sims_type = config.get("sims_type", "grid") + str_type = f"_{sims_type}" if sims_type == "grid" else "" + sim_dir = os.path.join(config["base"], f"1p2z{str_type}_{num}") + ini_path = os.path.join(sim_dir, "cfis", "config_tile_Git_symlink.ini") + + cp = configparser.ConfigParser(interpolation=None) + if not cp.read(ini_path): + raise FileNotFoundError( + f"{ini_path} not found; run 'snakemake init_all' first" + ) + sec = cp["GET_IMAGES_RUNNER"] + paths = [p.strip() for p in sec["INPUT_PATH"].split(",")] + patterns = [p.strip() for p in sec["INPUT_FILE_PATTERN"].split(",")] + exts = [e.strip() for e in sec["INPUT_FILE_EXT"].split(",")] + numbering = sec["INPUT_NUMBERING"].split(",")[0].strip() + + idx = next( + (i for i, p in enumerate(patterns) if "weight" in p.lower()), None + ) + if idx is None: + raise ValueError(f"no weight entry in INPUT_FILE_PATTERN of {ini_path}") + + weight_dir = paths[idx].replace("$SP_DIR", sim_dir) + return weight_dir, patterns[idx], exts[idx], numbering + + +def coverage_fraction(path, subsample): + """Fraction of nonzero pixels in a weight map, subsampled.""" + with fits.open(path) as hdul: + data = hdul[0].data[::subsample, ::subsample] + return float(np.count_nonzero(data) / data.size) + + +def main(): + args = parse_args() + + with open(args.config) as f: + config = yaml.safe_load(f) + + threshold = args.threshold + if threshold is None: + threshold = config.get("min_tile_coverage", 0.5) + + num = config["num"] + tile_ids = get_tile_ids(config) + weight_dir, pattern, ext, numbering = get_weight_template(config) + out_path = args.output or os.path.join( + config["base"], f"tile_coverage_grid_{num}.yaml" + ) + + print(f"Scanning {len(tile_ids)} tiles in {weight_dir}") + print(f"Coverage threshold: {threshold}") + + coverage = {} + exclude = [] + missing = [] + for tile in tile_ids: + fname = re.sub(numbering, tile.replace(".", "-"), pattern) + ext + path = os.path.join(weight_dir, fname) + if not os.path.isfile(path): + missing.append(tile) + exclude.append(tile) + print(f" {tile}: WEIGHT FILE MISSING -> exclude") + continue + frac = coverage_fraction(path, args.subsample) + coverage[tile] = round(frac, 4) + if frac < threshold: + exclude.append(tile) + if args.verbose or frac < threshold: + flag = " <-- exclude" if frac < threshold else "" + print(f" {tile}: coverage {frac:.3f}{flag}") + + result = { + "num": num, + "min_coverage": threshold, + "coverage": coverage, + "missing": missing, + "exclude": sorted(exclude), + } + with open(out_path, "w") as f: + yaml.dump(result, f, default_flow_style=False, sort_keys=False) + + n_keep = len(tile_ids) - len(exclude) + print(f"\n{len(exclude)} of {len(tile_ids)} tiles below threshold" + f" ({n_keep} remain)") + print(f"Written to {out_path}") + + return 0 + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/scripts/image_sims_pipeline/config.yaml.template b/scripts/image_sims_pipeline/config.yaml.template new file mode 100644 index 000000000..8a68667d1 --- /dev/null +++ b/scripts/image_sims_pipeline/config.yaml.template @@ -0,0 +1,67 @@ +# Basic configuration + +## Type: data or image_sims +type: image_sims + +## Base run directory (-d flag) +base: /n17data/mkilbing/astro/Runs/shapepipe/CFIS/v2.0/image_sims/grids + +## Tile IDs path (-T flag) +tile_IDs: "/home/mkilbing/shapepipe/auxdir/CFIS/im_sims_202606/tile_numbers_{num}.txt" + +## Config directory for ShapePipe (-C flag); if omitted, type-specific default is used +config_dir: /home/mkilbing/shapepipe/example/cfis_image_sims + +## Optional: exclusive list of tiles to process (for testing) +## If provided, only these tiles will be used instead of all tiles from tile_IDs +## Leave empty/commented out to use all tiles from tile_IDs +# tile_IDs_exclusive: ["251.288", "251.289", "252.287"] + +## Optional: tiles to remove from the list (e.g. blank/corrupt input images) +## Low-coverage tiles are excluded automatically via "snakemake coverage"; +## use this for manual exclusions on top of that +# tile_IDs_exclude: ["281.281"] + +## Minimum input-tile coverage (nonzero weight fraction); tiles below are +## excluded after running "snakemake coverage" (writes tile_coverage_grid_{num}.yaml) +min_tile_coverage: 0.5 + +## Input simulation base dir (-i flag) +input_sims_base_dir: /n09data/hervas/skills_out/ + +## Input tiles and exposures directories (supports {subdir} and {input_sims_base_dir} placeholders) +input_tiles: "{input_sims_base_dir}/{subdir}/images/SP_tiles" +input_exp: "{input_sims_base_dir}/{subdir}/images/SP_exp" + +## Galaxy sample (-S flag) (for cuts and selection); catalogue subversion +sample: 9 + +## parameter file for extraction (-P flag) +params_py: /home/mkilbing/sp_validation/notebooks/params.py + +## SMP batch size per ShapePipe job (-1 = use config file default = 1) +n_smp: -1 + +## Simulation type +sims_type: grid + +## Grid number +num: 2 + + +## Calibration + +## Input shear amplitude +shear_amplitude: 0.025 + +## Position matching (for grids) +match_radius_deg: 0.0002 + +## Pair-match objects between sheared catalogues; if false, use all objects +pair_match: true + +## Weight schemes for m-bias; "none" = uniform weights; first entry is primary +w_cols: ['none', 'w_iv'] + +## Bootstrap error estimation +n_bootstrap: 500 diff --git a/scripts/image_sims_pipeline/info.py b/scripts/image_sims_pipeline/info.py new file mode 100644 index 000000000..77a4ff669 --- /dev/null +++ b/scripts/image_sims_pipeline/info.py @@ -0,0 +1,329 @@ +"""Info script for image-simulation Snakemake pipeline.""" + +import os +import argparse +import yaml +import subprocess + + +SCRIPT_DIR = os.path.dirname(os.path.abspath(__file__)) +SNAKEFILE = os.path.join(SCRIPT_DIR, "Snakefile") + + +def load_config(configfile): + with open(configfile) as f: + return yaml.safe_load(f) + + +def load_tile_ids(cfg): + """Load tile IDs from config, handling both lists and file paths.""" + num = cfg.get("num", "") + tile_ids = cfg.get("tile_IDs", None) + + if tile_ids is None: + # Auto-construct from num + tile_ids = os.path.expanduser( + f"~/shapepipe/auxdir/CFIS/im_sims_202606/tile_numbers_{num}.txt" + ) + + if isinstance(tile_ids, str): + # Expand {num} placeholder + tile_ids = tile_ids.replace("{num}", str(num)) + if os.path.isfile(tile_ids): + with open(tile_ids) as f: + tile_ids = [line.strip() for line in f if line.strip()] + else: + tile_ids = [tile_ids] + + if not isinstance(tile_ids, list): + tile_ids = [tile_ids] + + # Same filtering as in the Snakefile + exclusive = cfg.get("tile_IDs_exclusive") + if exclusive and isinstance(exclusive, list): + tile_ids = exclusive + exclude = cfg.get("tile_IDs_exclude") + if exclude and isinstance(exclude, list): + exclude_set = set(str(t) for t in exclude) + tile_ids = [t for t in tile_ids if t not in exclude_set] + + # Low-coverage tiles from the coverage scan, if present + coverage_file = os.path.join(cfg["base"], f"tile_coverage_grid_{num}.yaml") + if os.path.isfile(coverage_file): + with open(coverage_file) as f: + cov_exclude = set( + str(t) for t in (yaml.safe_load(f) or {}).get("exclude", []) + ) + tile_ids = [t for t in tile_ids if t not in cov_exclude] + + return tile_ids + + +def get_hdf5_tile_count(hdf5_path): + """Query tile count from HDF5 file using create_final_cat.py -l. + + Returns the count as a string, or "." if file doesn't exist or query fails. + """ + if not os.path.isfile(hdf5_path): + return "." + + try: + HOME = os.path.expanduser("~") + script = f"{HOME}/astro/repositories/github/shapepipe/scripts/python/create_final_cat.py" + result = subprocess.run( + [script, "-I", "-l", "-m", hdf5_path, "-v"], + capture_output=True, + text=True, + timeout=10, + ) + + if result.returncode != 0: + return "?" + + for line in result.stdout.split('\n'): + if line.startswith("Total:"): + parts = line.split() + if len(parts) >= 2: + return parts[1] + + return "?" + except Exception: + return "?" + + +def get_extract_tile_count(sim_path): + """Count unique tile IDs in shape_catalog_comprehensive_ngmix.hdf5. + + Returns the count as a string, or "." if file doesn't exist or count not found. + """ + try: + import h5py + import numpy as np + except ImportError: + return "?" + + hdf5_path = os.path.join(sim_path, "shape_catalog_comprehensive_ngmix.hdf5") + if not os.path.isfile(hdf5_path): + return "." + + try: + with h5py.File(hdf5_path, 'r') as hf: + tile_ids = hf['data']['TILE_ID'][:] + return str(len(np.unique(tile_ids))) + except Exception: + return "?" + + +def get_bias_tile_count(results_dir): + """Return max n_tiles from mbias_cumulative.yaml, or '.' if absent.""" + path = os.path.join(results_dir, "mbias_cumulative.yaml") + if not os.path.isfile(path): + return "." + try: + with open(path) as f: + data = yaml.safe_load(f) or {} + if not data: + return "." + return str(max(int(k) for k in data.keys())) + except Exception: + return "?" + + +def monitor(cfg, verbose=0): + base = cfg["base"] + num = cfg["num"] + sim_type = cfg["sims_type"] + tile_ids = load_tile_ids(cfg) + n_tiles = len(tile_ids) + job_seq = [1, 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024, 2048] + job_labels = ["1","2","4","8","16","32","64","1h","2h","5h","1k","2k"] + job_names = ["Gt","Uz","Fe","Ge","Sp","Ma","Ps","Mh","Sx","Vi","Ng","Mc"] + job_widths = [max(len(h), 2) for h in job_labels] + + str_type = f"_{sim_type}" if sim_type == "grid" else "" + sims = [ + f"1{d1}2{d2}{str_type}_{num}" + for d1, d2 in zip(["m","p","z","z"], ["z","z","m","p"]) + ] + grids = base + + col_w = max(len(s) for s in sims) + 2 + downstream = [ + ("merge", "final_cat_{sim}.hdf5"), + ("extr", "shape_catalog_comprehensive_ngmix.hdf5"), + ("calib", "shape_catalog_cut_ngmix.fits"), + ("diag", "results/footprint.png"), + ("bias", "results/mbias_cumulative.yaml"), + ] + ds_w = max(len(s) for s, _ in downstream) + + tile_label = f"({n_tiles} tile{'s' if n_tiles > 1 else ''})" + mid_head = " ".join(f"{h:>{w}}" for h, w in zip(job_labels, job_widths)) + name_head = " ".join(f"{n:>{w}}" for n, w in zip(job_names, job_widths)) + bits_w = len(mid_head) + header1 = f" {'sim':<{col_w}} {'n':<{bits_w}} " + \ + " ".join(f"{s:<{ds_w}}" for s, _ in downstream) + header2 = f" {'':^{col_w}} {mid_head}" + header3 = f" {'':^{col_w}} {name_head}" + print(f" {tile_label}") + print(header1) + print(header2) + print(header3) + print("-" * len(header1)) + + for sim in sims: + logs_dir = os.path.join(grids, sim, "logs") + + # For each job bit, count how many tiles have the log file + tile_counts = {} + for j in job_seq: + tile_counts[j] = sum( + 1 for t in tile_ids + if os.path.isfile(os.path.join(logs_dir, f"log_job_{t}_{j}.txt")) + ) + + bits_str = " ".join( + f"{'o' if tile_counts[j] == n_tiles else '.':>{w}}" + for j, w in zip(job_seq, job_widths) + ) + + ds_status = [] + n_tiles_extr = None + for idx, (name, pattern) in enumerate(downstream): + sim_path = os.path.join(grids, sim) + if name in ("diag", "bias"): + # global outputs, not per-sim + path = os.path.join(grids, pattern.format(sim=sim)) + else: + path = os.path.join(grids, sim, pattern.format(sim=sim)) + if name == "merge": + ds_status.append(get_hdf5_tile_count(path)) + elif name == "extr": + n_tiles_extr = get_extract_tile_count(sim_path) + ds_status.append(n_tiles_extr) + elif name == "calib": + if os.path.isfile(path): + ds_status.append(n_tiles_extr if n_tiles_extr is not None else "o") + else: + ds_status.append(".") + elif name == "bias": + ds_status.append(get_bias_tile_count(os.path.join(grids, "results"))) + else: + ds_status.append("o" if os.path.isfile(path) else ".") + + if verbose >= 1: + counts_str = " ".join( + f"{tile_counts[j]:>{w}}" + for j, w in zip(job_seq, job_widths) + ) + row = f" {sim:<{col_w}} {counts_str} " + \ + " ".join(f"{s:>{ds_w}}" for s in ds_status) + else: + row = f" {sim:<{col_w}} {bits_str} " + \ + " ".join(f"{s:>{ds_w}}" for s in ds_status) + print(row) + + +def parse_args(): + parser = argparse.ArgumentParser( + description="Info script for image-simulation Snakemake pipeline." + ) + parser.add_argument( + "-c", "--config", default="config.yaml", + help="path to config file (default: %(default)s)" + ) + parser.add_argument( + "-m", "--monitor", action="store_true", + help="monitor pipeline progress" + ) + parser.add_argument( + "-v", "--verbose", action="count", default=0, + help="increase verbosity (can be repeated)" + ) + return parser.parse_args() + + +def main(): + args = parse_args() + + # config is only loaded after argument parsing, so -h works without it + cfg = load_config(args.config) + base_cmd = f"snakemake -s {SNAKEFILE} --configfile {args.config}" + + if args.monitor: + monitor(cfg, verbose=args.verbose) + return + + base = cfg["base"] + tile_ids = load_tile_ids(cfg) + sample = cfg["sample"] + num = cfg["num"] + sim_type = cfg["sims_type"] + job_seq = [1, 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024, 2048] + + str_type = f"_{sim_type}" if sim_type == "grid" else "" + sims = [ + f"1{d1}2{d2}{str_type}_{num}" + for d1, d2 in zip(["m","p","z","z"], ["z","z","m","p"]) + ] + grids = base + + print("=" * 70) + print(" Image simulations pipeline") + print("=" * 70) + print() + print(f" Config : {args.config}") + print(f" Base : {base}") + print(f" Tile IDs : {tile_ids}") + print(f" Sample : {sample} (mask config mask_v1.X.{sample}_im_sim.yaml)") + print(f" Run num : {num}") + print(f" Sims : {', '.join(sims)}") + print(f" Job seq : {job_seq}") + print() + + print("# Basic command") + print() + print(f" {base_cmd} -j 4 --rerun-triggers=mtime -- ") + print() + + print("# Target order") + print() + stages = [ + ("init_all", "initialise run directories"), + ("coverage", "scan input tiles, exclude (near-)empty ones (once per grid)"), + ("pipeline_all", f"run full job sequence {job_seq} in order"), + ("merge_all", "merge ShapePipe output into final_cat_.hdf5"), + ("extract_all", "extract comprehensive shape catalogue (HDF5)"), + ("calibrate_all", "apply masks and calibrate (= all)"), + ("diagnostics", "compute and plot diagnostics"), + ("mbias", "compute multiplicative biases") + ] + for target, desc in stages: + print(f" {target:<16} {desc}") + print() + + print("# Useful options") + print() + print(f" -p print shell commands") + print(f" --config 'force=1' pass --force to run_job_sp_canfar_v2.0.bash") + print(f" --forcerun pipeline_all force Snakemake to rerun all pipeline jobs") + print(f" --config 'job=J' target a specific job bit J") + print(f" --rerun-incomplete if errors occur due to previous interrupted run") + print(f" pipeline --config 'job=J' --rerun-incomplete check and fix runs") + print() + + print("# Monitor progress") + print() + print(f" python info_image_sims.py -m") + print() + + print("# Run a single shapepipe job for one sim") + print() + print(f" {base_cmd} -j 1 \\") + print(f" {grids}//logs/log_job_J.txt") + print() + + + +if __name__ == "__main__": + main() diff --git a/scripts/image_sims_pipeline/monitor_mbias.py b/scripts/image_sims_pipeline/monitor_mbias.py new file mode 100644 index 000000000..848562ec8 --- /dev/null +++ b/scripts/image_sims_pipeline/monitor_mbias.py @@ -0,0 +1,264 @@ +#!/usr/bin/env python +"""Monitor m-bias convergence as pipeline_all tiles complete. + +For tiles that have finished all ShapePipe jobs (bit 2048) in every sim, +runs merge/extract/calibrate in a monitoring subdirectory, then computes +m1, m2, c1, c2 and reports them alongside the tile count. + +Run repeatedly while pipeline_all is in progress: merge is incremental +(adds new tiles, skips tiles already in the HDF5), so each run picks up +newly finished tiles automatically. + +Usage (from run directory): + python ~/astro/repositories/github/shapepipe/scripts/image_sims_pipeline/monitor_mbias.py [-c config.yaml] [-v] +""" + +import sys +import os +import subprocess +import shutil +import yaml +import argparse + + +CONFIGFILE = "config_image_sims.yaml" +JOB_LAST = 2048 + +HOME = os.path.expanduser("~") +PATH_SH = f"{HOME}/astro/repositories/github/shapepipe/scripts/sh" +PATH_PY = f"{HOME}/astro/repositories/github/shapepipe/scripts/python" +PATH_NB = f"{HOME}/astro/repositories/github/sp_validation/notebooks" +PATH_SPV = f"{HOME}/astro/repositories/github/sp_validation" +SPV_SRC = f"{PATH_SPV}/src" +SPV_SIF = "/n17data/mkilbing/sp_validation_im_sims.sif" +APPTAINER = f"{PATH_SH}/apptainer_noslurm.sh" + +SPV_PREFIX = [ + APPTAINER, + "--bind", "/n17data,/n09data,/home", + "--env", f"PYTHONPATH={SPV_SRC}", + SPV_SIF, + "python", +] + + +def parse_args(): + p = argparse.ArgumentParser(description=__doc__) + p.add_argument( + "-c", "--config", default=CONFIGFILE, + help=f"Snakemake config YAML (default: {CONFIGFILE})" + ) + p.add_argument("-v", "--verbose", action="store_true") + return p.parse_args() + + +def load_config(path): + with open(path) as f: + return yaml.safe_load(f) + + +def sims_from_config(cfg): + num = cfg["num"] + str_type = f"_{cfg['type']}" if cfg["type"] == "grid" else "" + return [ + f"1{d1}2{d2}{str_type}_{num}" + for d1, d2 in zip(["m", "p", "z", "z", "z"], + ["z", "z", "m", "p", "z"]) + ] + + +def done_tiles(grids, sims, tile_ids): + """Return tiles that have log_job_{tile}_2048.txt in every sim.""" + result = [] + for tile in tile_ids: + if all( + os.path.isfile(f"{grids}/{sim}/logs/log_job_{tile}_{JOB_LAST}.txt") + for sim in sims + ): + result.append(tile) + return result + + +def setup_sim_dir(mon_sim_dir, real_sim_dir, verbose): + """Prepare monitoring sim dir with symlinks and copied params.py.""" + os.makedirs(mon_sim_dir, exist_ok=True) + + for name in ("cfis", "config_mask.yaml"): + link = os.path.join(mon_sim_dir, name) + if not os.path.exists(link): + target = os.path.join(real_sim_dir, name) + if os.path.exists(target): + os.symlink(target, link) + elif verbose: + print(f" Warning: {target} not found, skipping symlink") + + params_src = os.path.join(real_sim_dir, "params.py") + params_dst = os.path.join(mon_sim_dir, "params.py") + if os.path.exists(params_src): + shutil.copy2(params_src, params_dst) + elif verbose: + print(f" Warning: {params_src} not found") + + +def run_spv(script_args, cwd, label, verbose): + """Run a script inside the SPV container; return True on success.""" + cmd = SPV_PREFIX + script_args + if verbose: + print(f" cmd: {' '.join(cmd)}") + print(f" cwd: {cwd}") + result = subprocess.run( + cmd, cwd=cwd, + capture_output=(not verbose), + text=True, + ) + if result.returncode != 0: + print(f" {label} FAILED (exit {result.returncode})") + if not verbose and result.stderr: + print(result.stderr[-2000:]) + return False + return True + + +def run_merge(mon_sim_dir, grids, sim, verbose): + return run_spv( + [ + f"{PATH_PY}/create_final_cat.py", + "-I", + "-m", f"final_cat_{sim}.hdf5", + "-i", grids, + "-P", sim, + "-p", f"{grids}/{sim}/cfis/final_cat.param", + "-o", "n_tiles_monitor.txt", + ] + (["-v"] if verbose else []), + cwd=mon_sim_dir, + label="merge", + verbose=verbose, + ) + + +def run_extract(mon_sim_dir, verbose): + return run_spv( + [f"{PATH_NB}/extract_info.py"], + cwd=mon_sim_dir, + label="extract", + verbose=verbose, + ) + + +def run_calibrate(mon_sim_dir, verbose): + return run_spv( + [f"{PATH_NB}/calibrate_comprehensive_cat.py", "-s", "calibrate"], + cwd=mon_sim_dir, + label="calibrate", + verbose=verbose, + ) + + +def run_mbias(mon_dir, num, verbose): + cfg_path = os.path.join(mon_dir, "image_sims_m_bias_monitor.yaml") + out_path = os.path.join(mon_dir, "m_bias_monitor.yaml") + mbias_cfg = { + "base": mon_dir, + "num": num, + "catalog_name": "shape_catalog_cut_ngmix.fits", + "shear_amplitude": 0.02, + "match_radius_deg": 0.0002, + "w_col": "w_des", + "n_bootstrap": 500, + "output_path": out_path, + } + with open(cfg_path, "w") as f: + yaml.dump(mbias_cfg, f, default_flow_style=False) + + ok = run_spv( + [f"{PATH_SPV}/scripts/compute_m_bias_image_sims.py", "-c", cfg_path] + + (["-v"] if verbose else []), + cwd=mon_dir, + label="mbias", + verbose=verbose, + ) + return ok, out_path + + +def n_tiles_from_file(path): + """Read tile count written by create_final_cat.py -o.""" + try: + with open(path) as f: + return int(f.read().strip()) + except Exception: + return None + + +def main(): + args = parse_args() + cfg = load_config(args.config) + + base = cfg["base"] + num = cfg["num"] + tile_ids = cfg.get("tile_IDs", [cfg.get("tile_ID", "254.286")]) + grids = f"{base}/{cfg['type']}s" + sims = sims_from_config(cfg) + mon_dir = os.path.join(grids, "monitoring") + + # -- Detect done tiles ----------------------------------------------- + done = done_tiles(grids, sims, tile_ids) + n_done = len(done) + n_total = len(tile_ids) + + print(f"Tiles finished in all sims: {n_done} / {n_total}") + if done: + print(f" {done}") + + if n_done == 0: + print("No tiles complete yet. Nothing to compute.") + return 0 + + # -- Set up monitoring directory ------------------------------------ + print(f"\nMonitoring dir: {mon_dir}") + os.makedirs(mon_dir, exist_ok=True) + + # -- Per-sim: merge → extract → calibrate -------------------------- + for sim in sims: + real_sim_dir = os.path.join(grids, sim) + mon_sim_dir = os.path.join(mon_dir, sim) + print(f"\n{sim}") + + setup_sim_dir(mon_sim_dir, real_sim_dir, args.verbose) + + print(" merge ...", end="", flush=True) + if not run_merge(mon_sim_dir, grids, sim, args.verbose): + return 1 + n_found = n_tiles_from_file(os.path.join(mon_sim_dir, "n_tiles_monitor.txt")) + print(f" ok ({n_found} tiles in HDF5)") + + print(" extract ...", end="", flush=True) + if not run_extract(mon_sim_dir, args.verbose): + return 1 + print(" ok") + + print(" calibrate...", end="", flush=True) + if not run_calibrate(mon_sim_dir, args.verbose): + return 1 + print(" ok") + + # -- m-bias -------------------------------------------------------- + print("\nComputing m-bias...") + ok, results_path = run_mbias(mon_dir, num, args.verbose) + if not ok: + return 1 + + if os.path.isfile(results_path): + with open(results_path) as f: + res = yaml.safe_load(f) + print() + print(f" n_tiles = {n_done} / {n_total}") + print(f" m1 = {res['m1']:+.4f} ± {res['m1_err']:.4f}") + print(f" c1 = {res['c1']:+.4f} ± {res['c1_err']:.4f}") + print(f" m2 = {res['m2']:+.4f} ± {res['m2_err']:.4f}") + print(f" c2 = {res['c2']:+.4f} ± {res['c2_err']:.4f}") + + return 0 + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/scripts/python/create_final_cat.py b/scripts/python/create_final_cat.py index a0b10e8c9..75e3a67c0 100755 --- a/scripts/python/create_final_cat.py +++ b/scripts/python/create_final_cat.py @@ -61,7 +61,7 @@ def params_default(): Set default parameter values. """ - _params = { + _params = { "input_root_dir": ".", "merged_cat_path": "final_cat.hdf5", "param_path": None, @@ -71,29 +71,33 @@ def params_default(): "output_summary": "n_tiles_final.txt", "ID": None, "single_op": None, + "image_sims": False, } - _short_options = { - "input_root_dir": "-i", + _short_options = { + "input_root_dir": "-i", "merged_cat_path": "-m", "param_path": "-p", "patch": "-P", "list_only": "-l", "output_summary": "-o", "single_op": "-s", - } - _types = { - "hdu_num": "int", + "image_sims": "-I", + } + _types = { + "hdu_num": "int", "list_only": "bool", - } - _help_strings = { + "image_sims": "bool", + } + _help_strings = { "input_root_dir": "input root_dir for tile catalogues, default={}", "merged_cat_path": "merged catalogue path (hdf5 file), default={}", "param_path": "parameter file path, if not given use all columns, default={}", - "patch": "patch number, default={}", + "patch": "patch number (data) or grid subdir (image_sims), default={}", "list_only": "print list of patches and IDs only, default={}", "output_summary": "output file for numbre of tiles, default={}", "ID": "ID for single-ID operation, default={}", "single_op": "single ID operation, allowed are 'check', 'add', 'remove'; default={}", + "image_sims": "image simulations mode (different dir layout and run prefix), default={}", } return _params, _short_options, _types, _help_strings @@ -231,17 +235,34 @@ def check_ID(merged_cat_path, ID, verbose=False): def print_list(params): + verbose = params.get("verbose", False) n_tiles = 0 + + if not os.path.exists(params["merged_cat_path"]): + print(f"File {params['merged_cat_path']} not found") + return + with h5py.File(params["merged_cat_path"], "r") as hdf5_file: - for patch in hdf5_file["patches"]: - #print(patch) - for id in hdf5_file[f"patches/{patch}"]: - n_tiles += 1 + if "patches" not in hdf5_file: + print("Warning: no 'patches' group in output file (0 tiles added?)") + else: + for patch in hdf5_file["patches"]: + for id in hdf5_file[f"patches/{patch}"]: + n_tiles += 1 + if verbose: + print(f" {patch}/{id}") + + if verbose: + print(f"Total: {n_tiles} tiles") with open(params["output_summary"], "w") as f_out: print(n_tiles, file=f_out) + # Write n_tiles to HDF5 file header + with h5py.File(params["merged_cat_path"], "a") as hdf5_file: + hdf5_file.attrs["n_tiles"] = n_tiles + def get_patch_group(hdf5_file, patch, verbose=False): """Get Patch group. @@ -339,20 +360,50 @@ def copy_data(param_list, extracted_data, dtype): return structured_data +def collect_tile_ids_image_sims(patch_path): + """Collect tile IDs from image-sims layout: tiles/// + + Parameters + ---------- + patch_path : str + path to the grid subdir (e.g. .../1p2z_grid_1) + + Returns + ------- + list of (tile_id, tile_path) tuples + """ + id_pattern = re.compile(r"^\d+\.\d+$") + tiles_root = os.path.join(patch_path, "tiles") + result = [] + if not os.path.isdir(tiles_root): + return result + for prefix in os.listdir(tiles_root): + prefix_path = os.path.join(tiles_root, prefix) + if not os.path.isdir(prefix_path): + continue + for tile_id in os.listdir(prefix_path): + if id_pattern.match(tile_id): + result.append((tile_id, os.path.join(prefix_path, tile_id))) + return result + + def process(params): - - patch = rf"P{params['patch']}" - patch_pattern = re.compile(patch) - - # Define the nested path pattern for locating `.fits` files - file_pattern = "output/run_sp_Mc_*/make_cat_runner/output/*.fits" + + if params["image_sims"]: + patch_name = params["patch"] + run_prefix = "run_sp_tile_Mc_*" + else: + patch_name = rf"P{params['patch']}" + run_prefix = "run_sp_Mc_*" + + patch_pattern = re.compile(patch_name) # Regex pattern for tile IDs id_pattern = re.compile(r"^\d+\.\d+$") n_added = 0 IDs_added = [] - + # Open the HDF5 file (create it if it doesn't exist) if params["verbose"]: print(f"Initializing file {params['merged_cat_path']}") @@ -364,21 +415,30 @@ def process(params): # Skip non-matching entries if not patch_pattern.fullmatch(patch): continue - + # Full path to patch patch_path = os.path.join(params["input_root_dir"], patch) if not os.path.isdir(patch_path): if params["verbose"]: print(f"Path {patch_path} not found, skipping") continue - + # Get hdf5 group for this patch patch_group = get_patch_group(hdf5_file, patch, params["verbose"]) - # Get paths to all tile IDs and loop - tile_runs_path = os.path.join(patch_path, "tile_runs") - subdirs = os.listdir(tile_runs_path) - for id in tqdm.tqdm(subdirs, total=len(subdirs)): + # Collect (tile_id, tile_path) pairs depending on layout + if params["image_sims"]: + tile_items = collect_tile_ids_image_sims(patch_path) + else: + tile_runs_path = os.path.join(patch_path, "tile_runs") + tile_items = [ + (tid, os.path.join(tile_runs_path, tid)) + for tid in os.listdir(tile_runs_path) + if id_pattern.match(tid) + and os.path.isdir(os.path.join(tile_runs_path, tid)) + ] + + for id, id_path in tqdm.tqdm(tile_items, total=len(tile_items)): # Skip if the patch/ID data already exists if id in patch_group: @@ -386,44 +446,40 @@ def process(params): print(f"Skipping {id} (already processed)") continue - if id_pattern.match(id) and os.path.isdir(os.path.join(tile_runs_path, id)): - id_path = os.path.join(tile_runs_path, id) - - base_pattern = os.path.join(id_path, "output", "run_sp_Mc_*") - all_matches = [d for d in glob.glob(base_pattern) if os.path.isdir(d)] - if not all_matches: - if params["verbose"]: - print(f"Final cat for {id} not found, continuing") - continue - newest_dir = max(all_matches, key=os.path.getmtime) - - id_dash = re.sub("\.", "-", id) - fits_file = f"{newest_dir}/make_cat_runner/output/final_cat-{id_dash}.fits" - - # Exclude unsuccessful run without output FITS file - if not os.path.exists(fits_file): - if params["verbose"]: - print(f"Run without output file found for {id}, skipping") - continue - - if True: - extracted_data, dtype = read_data(fits_file, params) - - structured_data = copy_data(params["param_list"], extracted_data, dtype) - - # Create a new dataset - try: - patch_group.create_dataset( - str(id), - data=structured_data, - dtype=dtype, - ) - except: - print(f"Error for {id}: Could not create dataset in group {patch}") - raise - - n_added += 1 - IDs_added.append(id) + base_pattern = os.path.join(id_path, "output", run_prefix) + all_matches = [d for d in glob.glob(base_pattern) if os.path.isdir(d)] + if not all_matches: + if params["verbose"]: + print(f"Final cat for {id} not found, continuing") + continue + newest_dir = max(all_matches, key=os.path.getmtime) + + id_dash = re.sub(r"\.", "-", id) + fits_file = f"{newest_dir}/make_cat_runner/output/final_cat-{id_dash}.fits" + + # Exclude unsuccessful run without output FITS file + if not os.path.exists(fits_file): + if params["verbose"]: + print(f"Run without output file found for {id}, skipping") + continue + + extracted_data, dtype = read_data(fits_file, params) + + structured_data = copy_data(params["param_list"], extracted_data, dtype) + + # Create a new dataset + try: + patch_group.create_dataset( + str(id), + data=structured_data, + dtype=dtype, + ) + except: + print(f"Error for {id}: Could not create dataset in group {patch}") + raise + + n_added += 1 + IDs_added.append(id) if params["verbose"]: print(f"{n_added} tiles added ({' '.join(IDs_added)})") diff --git a/scripts/python/init_run_v2.0.py b/scripts/python/init_run_v2.0.py new file mode 100755 index 000000000..81b301ba5 --- /dev/null +++ b/scripts/python/init_run_v2.0.py @@ -0,0 +1,449 @@ +#!/usr/bin/env python3 + +""" +Initialise the directory structure for a v2.0 ShapePipe run. +Creates tile and exposure staging directories, config symlinks, +and a tile number list. + +Config file is loaded first (if provided), then overridden by command-line arguments. +""" + +import argparse +import sys +from pathlib import Path +import os + +try: + import yaml +except ImportError: + print("ERROR: PyYAML not installed. Install with: pip install pyyaml", file=sys.stderr) + sys.exit(1) + + +class InitRun: + def __init__(self): + self.sp_root = Path.home() / "shapepipe" + self.quiet = False + self.defaults = { + "type": "data", + "subdir": "", + "base_dir": ".", + "config_dir": "", + "params_src": str(Path.home() / "sp_validation" / "notebooks" / "params.py"), + "tiles_src": "", + "tile_ids_list": None, + "sample": "", + "input_sims_base_dir": "", + "input_tiles": "", + "input_exp": "", + "sims_type": "", + "num": "", + } + self.config = {} + self.args = {} + + def load_config(self, config_file): + """Load YAML config file.""" + config_path = Path(config_file) + if not config_path.exists(): + print(f"ERROR: Config file not found: {config_file}", file=sys.stderr) + sys.exit(2) + + with open(config_path) as f: + self.config = yaml.safe_load(f) or {} + + # Map config keys to script parameters + if "base" in self.config: + self.config["base_dir"] = self.config.pop("base") + if "type" in self.config: + # Keep type from config (data or image_sims) + self.config["type"] = self.config["type"] + if "tile_IDs" in self.config: + tile_ids = self.config.pop("tile_IDs") + # tile_IDs can be a file path (string) or a list of IDs + if isinstance(tile_ids, str): + self.config["tiles_src"] = tile_ids + elif isinstance(tile_ids, list): + # Will be written to file in run directory + self.config["tile_ids_list"] = tile_ids + if "params_py" in self.config: + self.config["params_src"] = self.config.pop("params_py") + # Keep input_tiles and input_exp if present (they may contain placeholders) + # They will be substituted before creating symlinks + + def merge_config(self): + """Merge: defaults <- config <- command-line args.""" + params = dict(self.defaults) + params.update(self.config) + params.update(self.args) + return params + + def substitute_placeholders(self, value, params): + """Substitute {key} placeholders in path strings using params dict.""" + if not isinstance(value, str): + return value + result = value + for key, val in params.items(): + if val is not None and val != "": + # Convert to string and strip trailing slashes + clean_val = str(val).rstrip('/') + result = result.replace(f"{{{key}}}", clean_val) + # Clean up any double slashes that may result + while '//' in result: + result = result.replace('//', '/') + return result + + def construct_default_subdir(self, params): + """Construct default subdir from sims_type and num if not explicitly set.""" + if not params["subdir"] and params.get("sims_type") and params.get("num"): + params["subdir"] = f"1p2z_{params['sims_type']}_{params['num']}" + + def set_type_defaults(self, params): + """Set type-specific defaults for config_dir and tiles_src.""" + if params["type"] == "data": + config_dir = self.sp_root / "example" / "cfis" + if not params["tiles_src"]: + params["tiles_src"] = str( + self.sp_root / "auxdir" / "CFIS" / "tiles_202604" / "tiles_r.txt" + ) + elif params["type"] == "image_sims": + config_dir = self.sp_root / "example" / "cfis_image_sims" + if not params["tiles_src"]: + # Try to use num-specific tile file first + if params.get("num"): + num_tile_file = self.sp_root / "auxdir" / "CFIS" / "im_sims_202606" / f"tile_numbers_{params['num']}.txt" + if num_tile_file.exists(): + params["tiles_src"] = str(num_tile_file) + else: + params["tiles_src"] = str( + self.sp_root / "auxdir" / "CFIS" / "im_sims_202606" / "numbers.txt" + ) + else: + params["tiles_src"] = str( + self.sp_root / "auxdir" / "CFIS" / "im_sims_202606" / "numbers.txt" + ) + else: + print(f"ERROR: Invalid type '{params['type']}'. Must be 'data' or 'image_sims'", file=sys.stderr) + sys.exit(3) + + # Use config_dir from config, or fall back to type-specific default + if not params["config_dir"]: + params["config_dir"] = str(config_dir) + return params + + def create_symlink(self, link_name, target, force=False): + """Create symlink, return True if created.""" + link_path = Path(link_name) + + if link_path.exists() or link_path.is_symlink(): + if force: + link_path.unlink() + else: + if link_path.is_symlink() and not self.quiet: + print(f"{link_name} symlink already exists, skipping") + elif not self.quiet: + print(f"WARNING: {link_name} exists as a regular file, not creating symlink") + return False + + link_path.symlink_to(target) + if not self.quiet: + print(f"Created symlink: {link_name} -> {target}") + return True + + def write_tile_ids(self, tile_ids_list): + """Write tile IDs list to tile_numbers.txt.""" + tile_file = Path("tile_numbers.txt") + with open(tile_file, "w") as f: + for tile_id in tile_ids_list: + f.write(f"{tile_id}\n") + if not self.quiet: + print(f"Created tile_numbers.txt with {len(tile_ids_list)} tiles") + + def copy_params(self, params, force=False): + """Copy and adapt params.py.""" + params_dest = Path("params.py") + + if params_dest.exists() and not force: + if not self.quiet: + print("params.py already exists, skipping") + return + + params_src = Path(params["params_src"]) + if not params_src.exists(): + print(f"ERROR: Source params.py not found: {params_src}", file=sys.stderr) + sys.exit(4) + + content = params_src.read_text() + # Remove IMAFLAGS_ISO + content = content.replace("'IMAFLAGS_ISO', ", "").replace(", 'IMAFLAGS_ISO'", "") + content = content.replace('"IMAFLAGS_ISO", ', "").replace(', "IMAFLAGS_ISO"', "") + + if params["type"] == "image_sims": + # Adapt for image_sims + lines = content.split("\n") + adapted = [] + for line in lines: + if line.startswith("name = "): + line = f"name = '{params['subdir']}'" + elif line.startswith("star_cat_path = "): + line = "star_cat_path = None" + elif line.startswith("output_format = "): + line = "output_format = '.hdf5'" + adapted.append(line) + content = "\n".join(adapted) + if not self.quiet: + print("Copied and adapted: params.py <- " + params["params_src"]) + print(f" name = '{params['subdir']}', star_cat_path = None, output_format = '.hdf5'") + else: + if not self.quiet: + print(f"Copied: params.py <- {params['params_src']}") + + params_dest.write_text(content) + + def run(self, config_file=None, quiet=False, **kwargs): + """Run initialization.""" + self.quiet = quiet + + # Load config first + if config_file: + self.load_config(config_file) + + # Apply command-line overrides (exclude quiet) + self.args = {k: v for k, v in kwargs.items() if v is not None and k != "quiet"} + + # Merge all sources + params = self.merge_config() + + # Construct default subdir from sims_type and num if needed + self.construct_default_subdir(params) + + # Validate that we have a subdir + if not params["subdir"]: + print("ERROR: subdir must be provided (-s) or derivable from sims_type and num in config", file=sys.stderr) + sys.exit(4) + + # Set type-specific defaults (including num-based tile_IDs for image_sims) + params = self.set_type_defaults(params) + + # Expand placeholders in tiles_src (e.g. {num}) + if params.get("tiles_src"): + params["tiles_src"] = self.substitute_placeholders(params["tiles_src"], params) + + # Use pwd if base_dir still empty + if not params["base_dir"]: + params["base_dir"] = str(Path.cwd()) + + base_dir = Path(params["base_dir"]) / params["subdir"] + + if not self.quiet: + print(f"Initialising ShapePipe run directory: {params['base_dir']}") + print() + + # Create base directory + base_dir.mkdir(parents=True, exist_ok=True) + os.chdir(base_dir) + + # Image sims specific: input symlinks + if params["type"] == "image_sims": + # input_tiles and input_exp symlinks (with placeholder substitution) + if params.get("input_tiles"): + input_tiles_path = self.substitute_placeholders(params["input_tiles"], params) + self.create_symlink("input_tiles", input_tiles_path, params["force"]) + if params.get("input_exp"): + input_exp_path = self.substitute_placeholders(params["input_exp"], params) + self.create_symlink("input_exp", input_exp_path, params["force"]) + + # Image sims specific: mask config symlink + if params["type"] == "image_sims" and params.get("sample"): + mask_src = ( + Path.home() + / "sp_validation" + / "config" + / "calibration" + / f"mask_v1.X.{params['sample']}_im_sim.yaml" + ) + self.create_symlink("config_mask.yaml", str(mask_src), params["force"]) + + # Create directories + if not self.quiet: + print("Creating tiles/ directory...") + Path("tiles").mkdir(exist_ok=True) + if not self.quiet: + print(" tiles/ created (subdirs added at download time)") + + if not self.quiet: + print("Creating exp/ directory...") + Path("exp").mkdir(exist_ok=True) + if not self.quiet: + print(" exp/ created (subdirs added at download time)") + + if not self.quiet: + print("Creating log and debug directories...") + Path("logs").mkdir(exist_ok=True) + Path("debug").mkdir(exist_ok=True) + + # Config symlink + self.create_symlink("cfis", params["config_dir"], params["force"]) + + # Tile numbers: either symlink to file or write list + if params.get("tile_ids_list"): + if not self.quiet: + print("Creating tile_numbers.txt from tile_IDs list...") + self.write_tile_ids(params["tile_ids_list"]) + else: + if not self.quiet: + print("Creating tile_numbers.txt symlink...") + self.create_symlink("tile_numbers.txt", params["tiles_src"], params["force"]) + n_tiles = sum(1 for _ in open("tile_numbers.txt")) + if not self.quiet: + print(f" {n_tiles} tiles") + + # params.py + if not self.quiet: + print("Creating params.py...") + self.copy_params(params, params["force"]) + + # Print summary + print(f"Done: {params['subdir']}") + + if not self.quiet: + print() + print("Directory structure:") + print(f" {base_dir}") + print(" ├── tiles/") + print(" ├── exp/") + print(" ├── logs/") + print(" ├── debug/") + print(f" ├── cfis -> {params['config_dir']}") + + # Show tile_numbers.txt info + if params.get("tile_ids_list"): + print(f" ├── tile_numbers.txt ({len(params['tile_ids_list'])} tiles)") + else: + print(f" ├── tile_numbers.txt -> {params['tiles_src']}") + + if params["type"] == "image_sims": + if params.get("sample"): + mask_src = ( + Path.home() + / "sp_validation" + / "config" + / "calibration" + / f"mask_v1.X.{params['sample']}_im_sim.yaml" + ) + print(f" ├── config_mask.yaml -> {mask_src}") + if params.get("input_tiles"): + input_tiles_path = self.substitute_placeholders(params["input_tiles"], params) + print(f" ├── input_tiles -> {input_tiles_path}") + if params.get("input_exp"): + input_exp_path = self.substitute_placeholders(params["input_exp"], params) + print(f" ├── input_exp -> {input_exp_path}") + + print(f" └── params.py <- {params['params_src']}") + + +def main(): + parser = argparse.ArgumentParser( + description="Initialise ShapePipe v2.0 run directory", + formatter_class=argparse.RawDescriptionHelpFormatter, + epilog=""" +Examples: + # Read all from config file (subdir auto-generated from sims_type + num) + init_run_v2.0.py -c config.yaml -t image_sims + + # Override config values with command-line options + init_run_v2.0.py -c config.yaml -t image_sims -d /custom/base -S 6 + + # Explicit subdir (overrides auto-generation) + init_run_v2.0.py -c config.yaml -t image_sims -s my_subdir + + # No config file, use defaults and flags + init_run_v2.0.py -t image_sims -d /tmp/test -s grid_1 + """, + ) + + parser.add_argument( + "-c", "--config", + help="Config file (YAML); entries can be overridden by following flags" + ) + parser.add_argument( + "-t", "--type", + choices=["data", "image_sims"], + help="Input type" + ) + parser.add_argument( + "-d", "--dir", + dest="base_dir", + help="Base run directory" + ) + parser.add_argument( + "-C", "--config-dir", + dest="config_dir", + help="Config directory (overrides type-specific default)" + ) + parser.add_argument( + "-s", "--subdir", + help="Subdirectory for simulations (auto-generated from sims_type+num if not provided)" + ) + parser.add_argument( + "-T", "--tiles", + dest="tiles_src", + help="Tile numbers file" + ) + parser.add_argument( + "-i", "--input-sims", + dest="input_sims_base_dir", + help="Input simulations base directory (image_sims only)" + ) + parser.add_argument( + "--input-tiles", + dest="input_tiles", + help="Input tiles directory path (supports {subdir} and {input_sims_base_dir} placeholders)" + ) + parser.add_argument( + "--input-exp", + dest="input_exp", + help="Input exposures directory path (supports {subdir} and {input_sims_base_dir} placeholders)" + ) + parser.add_argument( + "-S", "--sample", + help="Sample version for mask config (e.g., 6)" + ) + parser.add_argument( + "-P", "--params", + dest="params_src", + help="Source params.py to copy" + ) + parser.add_argument( + "--force", + action="store_true", + help="Recreate existing symlinks and parameter files" + ) + parser.add_argument( + "-q", "--quiet", + action="store_true", + help="Suppress all output except final 'Done' message" + ) + + args = parser.parse_args() + + init = InitRun() + init.run( + config_file=args.config, + type=args.type, + base_dir=args.base_dir, + config_dir=args.config_dir, + subdir=args.subdir, + tiles_src=args.tiles_src, + input_sims_base_dir=args.input_sims_base_dir, + input_tiles=args.input_tiles, + input_exp=args.input_exp, + sample=args.sample, + params_src=args.params_src, + force=args.force, + quiet=args.quiet, + ) + + +if __name__ == "__main__": + main() diff --git a/scripts/sh/apptainer_noslurm.sh b/scripts/sh/apptainer_noslurm.sh new file mode 100755 index 000000000..4d28dfb64 --- /dev/null +++ b/scripts/sh/apptainer_noslurm.sh @@ -0,0 +1,10 @@ +#!/bin/bash +# Run apptainer exec after stripping all SLURM/PMI/PMIX/OMPI environment +# variables that cause OpenMPI to crash when running inside a SLURM job. +# Usage: apptainer_noslurm.sh [apptainer_args...] -- [cmd_args...] + +for v in $(env | grep -Eo "^(SLURM|PMI|PMIX|OMPI)[^=]*"); do + unset "$v" +done + +exec apptainer exec "$@" diff --git a/scripts/sh/init_run_v2.0.sh b/scripts/sh/init_run_v2.0.sh deleted file mode 100755 index 74ab260c0..000000000 --- a/scripts/sh/init_run_v2.0.sh +++ /dev/null @@ -1,103 +0,0 @@ -#!/usr/bin/env bash - -# Name: init_run_v2.0.sh -# Description: Initialise the directory structure for a v2.0 ShapePipe run. -# Creates tile and exposure staging directories, config symlinks, -# and a tile number list from the 202604 tile catalogue. -# Author: Martin Kilbinger - -# Version -version="2.0" - -# Default base run directory (permanent storage) -#base_dir="$HOME/cosmostat/v2/v${version}" -base_dir=`pwd` - -# ShapePipe repository root (for config symlink and tile list) -sp_root="$HOME/shapepipe" - -# Tile list source (full filenames, will be stripped to NNN.MMM) -tiles_src="$sp_root/auxdir/CFIS/tiles_202604/tiles_r.txt" - -# Config directory (will be symlinked as $base_dir/cfis) -config_dir="$sp_root/example/cfis" - -## Help string -usage="Usage: $(basename "$0") [OPTIONS] -\nOptions:\n - -h\t\tthis message\n - -d, --dir DIR\tbase run directory, default='$base_dir'\n -" - -## Parse command line -while [ $# -gt 0 ]; do - case "$1" in - -h) - echo -ne $usage - exit 0 - ;; - -d|--dir) - base_dir="$2" - shift - ;; - *) - echo "Unknown option: $1" - echo -ne $usage - exit 1 - ;; - esac - shift -done - -echo "Initialising ShapePipe v${version} run directory: $base_dir" -echo "" - -# --- Base directory --- -mkdir -p "$base_dir" -cd "$base_dir" - -echo "Creating tiles/ directory..." -mkdir -p tiles -echo " exp/ created (subdirs e.g. tiles/301/301.279 added at download time)" - -# --- Exposure staging directory --- -echo "Creating exp/ directory..." -mkdir -p exp -echo " exp/ created (subdirs e.g. exp/28/ added at download time)" - -# --- Output and working directories --- -echo "Creating log and debug directories..." -mkdir -p logs -mkdir -p debug - -# --- Config symlink --- -if [ -L cfis ]; then - echo "cfis symlink already exists, skipping" -elif [ -d cfis ]; then - echo "WARNING: cfis/ exists as a directory, not creating symlink" -else - ln -s "$config_dir" cfis - echo "Created symlink: cfis -> $config_dir" -fi - -# --- Tile number list --- -echo "Creating tile_numbers.txt symlink..." -if [ -L tile_numbers.txt ]; then - echo " tile_numbers.txt symlink already exists, skipping" -elif [ -f tile_numbers.txt ]; then - echo " WARNING: tile_numbers.txt exists as a regular file, not creating symlink" -else - ln -s "$tiles_src" tile_numbers.txt - echo " Created symlink: tile_numbers.txt -> $tiles_src" -fi -n_tiles=$(wc -l < tile_numbers.txt) -echo " $n_tiles tiles" - -echo "" -echo "Done. Directory structure:" -echo " $base_dir/" -echo " ├── tiles/" -echo " ├── exp/" -echo " ├── logs/" -echo " ├── cfis -> ${config_dir}" -echo " └── tile_numbers.txt -> ${tiles_src}" diff --git a/scripts/sh/job_sp_canfar_v2.0.bash b/scripts/sh/job_sp_canfar_v2.0.bash index 1c1ce89ba..592cae70d 100755 --- a/scripts/sh/job_sp_canfar_v2.0.bash +++ b/scripts/sh/job_sp_canfar_v2.0.bash @@ -9,13 +9,17 @@ # machine. # Author: Martin Kilbinger +# Repository root, derived from this script's location (scripts/sh/ -> repo +# root) so the runner is relocatable: no hardcoded ~/shapepipe clone layout. +SP_ROOT="$(cd "$(dirname "${BASH_SOURCE[0]}")/../.." && pwd)" + # Shared job-list description -source $HOME/shapepipe/scripts/sh/job_list_help.bash +source "$SP_ROOT/scripts/sh/job_list_help.bash" # Command line arguments ## Default values job=255 -config_dir=$HOME/shapepipe/example/cfis +config_dir="$SP_ROOT/example/cfis" psf='mccd' retrieve='vos' star_cat_for_mask='onthefly' @@ -134,9 +138,8 @@ if [ "$star_cat_for_mask" != "onthefly" ] && [ "$star_cat_for_mask" != "save" ]; exit 4 fi -#if [ "$retrieve" != "vos" ] && [ "$retrieve" != "symlink" ]; then -if [ "$retrieve" != "vos" ]; then - echo "method to retrieve images (option -r) needs to be 'vos' for v2.0" +if [ "$retrieve" != "vos" ] && [ "$retrieve" != "symlink" ]; then + echo "method to retrieve images (option -r) needs to be 'vos' (data) or 'symlink' (image sims)" exit 5 fi @@ -291,7 +294,7 @@ function command_cfg_shapepipe() { batch_flag="--batch_size $_n_smp" fi - local cmd="shapepipe_run.py -c $config $batch_flag" + local cmd="shapepipe_run -c $config $batch_flag" command "$cmd" "$str" } @@ -347,12 +350,12 @@ if [[ $do_job != 0 ]]; then fi -## Retrieve exposure images (online, vos) +## Retrieve exposure images (online if retrieve=vos, symlinked if image sims) (( do_job = $job & 8 )) if [[ $do_job != 0 ]]; then command_cfg_shapepipe \ - "config_exp_Gie_vos.ini" \ + "config_exp_Gie_$retrieve.ini" \ "Run shapepipe (get exposure images)" \ $n_smp \ $exclusive @@ -495,11 +498,13 @@ fi (( do_job = $job & 2048 )) if [[ $do_job != 0 ]]; then - suff_sm="_nosm" - - ### Merge all relevant information into final catalogue + ### Merge all relevant information into final catalogue. + ### Config name matches run_job's job-2048 "Mc_${psf}" completeness check and + ### the config present in example/cfis_image_sims. (The data path uses a + ### differently-named make_cat config; harmonizing the two dirs is a separate + ### #766 reconciliation item, not needed for the image-sims rerun.) command_cfg_shapepipe \ - "config_tile_make_cat_$psf${suff_sm}.ini" \ + "config_tile_Mc_$psf.ini" \ "Run shapepipe (tile: create final cat $psf)" \ $n_smp \ $exclusive diff --git a/scripts/sh/run_job_sp_canfar_v2.0.bash b/scripts/sh/run_job_sp_canfar_v2.0.bash index 9da08f598..da207bbcc 100755 --- a/scripts/sh/run_job_sp_canfar_v2.0.bash +++ b/scripts/sh/run_job_sp_canfar_v2.0.bash @@ -3,8 +3,13 @@ # run_job_sp_canfar_v2.0 # Description: Initialise tile/exposure run directory and launch ShapePipe job. +# Repository root, derived from this script's location (scripts/sh/ -> repo root) +# so the runner is relocatable: no hardcoded ~/shapepipe clone layout. All +# config dirs and sibling scripts are resolved under $SP_ROOT. +SP_ROOT="$(cd "$(dirname "${BASH_SOURCE[0]}")/../.." && pwd)" + # Shared job-list description -source $HOME/shapepipe/scripts/sh/job_list_help.bash +source "$SP_ROOT/scripts/sh/job_list_help.bash" # ShapePipe version version="2.0" @@ -16,15 +21,22 @@ ID=-1 psf='psfex' tile_det='uc' tile_mask=0 -N_SMP=1 +N_SMP=-1 dry_run=0 dir=`pwd` debug_out="" + +# Input type: data or image_sims +type="data" + #scratch="/scratch/$USER/shapepipe/v${version}" scratch="" test_only=0 -check=0 +check_only=0 +run_and_check=1 force=0 +retry=0 +quiet=0 VERBOSE=1 pat="-- " @@ -40,14 +52,18 @@ ${JOB_LIST_HELP} -e, --exclusive ID\timage ID\n \t\t\tPSF model, one in ['psfex'|'mccd'], default='$psf'\n --tile_det DET\t\ttile detection mode, one in ['sx'|'uc'], default='$tile_det'\n --tile_mask MASK\ttile masking, default='$tile_mask'\n + -t, --type TYPE input type, allowed are 'data', 'image_sims', default='$type'\n -N, --N_SMP N_SMP\tnumber of SMP jobs, default from original config files\n -d, --directory DIR\trun directory, default is pwd ($dir)\n -S, --scratch DIR\tprocessing scratch directory, default=none\n -n, --dry_run\t\tDRY RUN, no actual processing; default is $dry_run\n --debug_out PATH\tdebug output file PATH, default=none\n --test\t\t\ttest mode, no processing\n - --check\t\tcheck download completeness only (job 8), no processing\n + --check_only\t\tcheck completeness only, no processing\n --force\t\tremove existing module output dir(s) before running\n + --retry\t\tskip jobs whose existing run is complete; remove and rerun\n + \t\t\tonly those whose existing run is incomplete/failed\n + -q, --quiet\t\tsuppress all output except Complete/Missing/Incomplete/WARNING lines\n " ## Help if no arguments @@ -75,6 +91,10 @@ while [ $# -gt 0 ]; do psf="$2" shift ;; + -t|--type) + type="$2" + shift + ;; --tile_det) tile_det="$2" shift @@ -106,12 +126,18 @@ while [ $# -gt 0 ]; do --test) test_only=1 ;; - --check) - check=1 + --check_only) + check_only=1 ;; --force) force=1 ;; + --retry) + retry=1 + ;; + -q|--quiet) + quiet=1 + ;; esac shift done @@ -124,10 +150,15 @@ function message() { my_debug_out=$2 my_exit=$3 - echo $msg + if [ "$quiet" == "0" ] || [[ "$msg" =~ (Complete|[Mm]issing|Incomplete|WARNING) ]]; then + echo $msg + fi if [ -n "$my_debug_out" ]; then echo ${pat}$msg >> $my_debug_out fi + if [ -n "$log_file" ]; then + echo ${pat}$msg >> $log_file + fi if [ "$my_exit" != "-1" ]; then if [ -n "$my_debug_out" ]; then @@ -135,12 +166,14 @@ function message() { else echo "${pat}exiting with code $my_exit" fi + if [ -n "$log_file" ]; then + echo "${pat}exiting with code $my_exit" >> $log_file + fi exit $my_exit fi } - # Initialise exposure work directory: create dirs, exp_numbers file, config symlink. # The exp_numbers-000-000.txt file is created only once (skipped if already exists). # Args: $1 = exp_id, $2 = exp_work_dir @@ -158,11 +191,35 @@ function init_exp_work_dir() { fi if [ ! -e "$exp_work_dir/cfis" ]; then - ln -sf ~/shapepipe/example/cfis "$exp_work_dir/cfis" + ln -sf $config_dir "$exp_work_dir/cfis" fi } +# Remove existing run_sp__* output directories. +# Args: $1 = base dir containing the "output" subdir (work_dir or exp_work_dir) +# $2 = "tile" or "exp" +# $3 = space-separated run_prefixes +# $4 = label for the removal message (e.g. "Force-removing") +function remove_run_dirs() { + local base_dir=$1 + local kind=$2 + local run_prefixes=$3 + local label=$4 + local run_prefix + for run_prefix in $run_prefixes; do + local dirs_to_remove + dirs_to_remove=$(ls -d "$base_dir/output/run_sp_${kind}_${run_prefix}"* 2>/dev/null) + if [ -n "$dirs_to_remove" ]; then + for d in $dirs_to_remove; do + message "${label} $d" "$debug_out" -1 + command "rm -rf $d" $dry_run + done + fi + done +} + + # Run a per-exposure job (e.g. job 8, 16). # Args: $1 = job number # $2 = space-separated list of run_sp_exp output dir prefixes (e.g. "Gie") @@ -185,7 +242,7 @@ function run_exp_job() { message "Exposure numbers file exp_numbers-${IDra}-${IDdec}.txt not found in $work_dir/output" "$debug_out" 10 fi - if [ "$check" == "1" ]; then + if [ "$check_only" == "1" ]; then message "Check mode: skipping job $exp_job" "$debug_out" -1 fi @@ -200,30 +257,27 @@ function run_exp_job() { (( n_total++ )) - # exp_id e.g. "2182795p": ab = first 2 chars, abcdefg = all but last char + # exp_id e.g. "2182795p" (data) or "208659" (image_sims) + # Strip trailing letter if present (data format); keep full id if numeric only. local exp_prefix="${exp_id:0:2}" - local exp_base="${exp_id%?}" + local exp_base + if [[ "${exp_id: -1}" =~ [a-zA-Z] ]]; then + exp_base="${exp_id%?}" + else + exp_base="$exp_id" + fi local exp_id_disp="${exp_prefix}/${exp_base}" - local exp_work_dir="$HOME/v${version}/exp/$exp_prefix/$exp_base" + local exp_work_dir="$dir/exp/$exp_prefix/$exp_base" local exp_log_file="$exp_work_dir/job_sp_canfar_v2.0.log" # Create exp_numbers-000-000.txt and cfis link if not existent init_exp_work_dir "$exp_id" "$exp_work_dir" - # force: remove all existing run directories for each prefix before running - if [ "$force" == "1" ]; then - local run_prefix - for run_prefix in $run_prefixes; do - local dirs_to_remove - dirs_to_remove=$(ls -d "$exp_work_dir/output/run_sp_exp_${run_prefix}"* 2>/dev/null) - if [ -n "$dirs_to_remove" ]; then - for d in $dirs_to_remove; do - message "Force-removing $d" "$debug_out" -1 - command "rm -rf $d" $dry_run - done - fi - done - fi + # Acquire a per-exposure lock so parallel tiles that share this exposure + # don't race: one runs, the other waits and then finds it already complete. + local lock_fd + exec {lock_fd}>"${exp_work_dir}/.sp_lock" + flock -x "$lock_fd" # Check completeness of existing run output (main prefix) local run_dir=$(ls -dt "$exp_work_dir/output/run_sp_exp_${main_prefix}"* 2>/dev/null | head -1) @@ -286,20 +340,30 @@ function run_exp_job() { fi done - if [ "$is_complete" == "1" ]; then + # force: always remove and rerun. + # retry: remove and rerun only if the existing run is incomplete/failed. + if [ "$force" == "1" ]; then + remove_run_dirs "$exp_work_dir" "exp" "$run_prefixes" "Force-removing" + elif [ "$retry" == "1" ] && [ "$is_complete" == "0" ]; then + remove_run_dirs "$exp_work_dir" "exp" "$run_prefixes" "Retry: removing incomplete run" + fi + + if [ "$force" != "1" ] && [ "$is_complete" == "1" ]; then message "Complete $exp_id_disp: run_sp_exp_${main_prefix} ( $check_desc)" "$debug_out" -1 (( n_complete++ )) + exec {lock_fd}>&- continue fi # Report incomplete/missing in check mode; in run mode handle and proceed - if [ "$check" == "1" ]; then + if [ "$check_only" == "1" ]; then if [ -n "$run_dir" ]; then message " Benign incomplete: $exp_id_disp ($check_desc)" "$debug_out" -1 else message " missing: $exp_id_disp" "$debug_out" -1 fi (( n_incomplete++ )) + exec {lock_fd}>&- continue fi @@ -311,9 +375,10 @@ function run_exp_job() { [ -n "$debug_out" ] && debug_flag="--debug_out $debug_out" echo "$(basename "$0") -j $exp_job -e $exp_id" > "$exp_log_file" - echo "pwd=`pwd`" - command "job_sp_canfar_v2.0.bash -p $psf --tile_det $tile_det --tile_mask $tile_mask -j $exp_job --n_smp $N_SMP --nsh_jobs $N_SMP $debug_flag" $dry_run 2>&1 | tee -a "$exp_log_file" - echo "Done with job_sp_canfar_v2.0.bash" + [ "$quiet" == "0" ] && echo "pwd=`pwd`" + command "$SP_ROOT/scripts/sh/job_sp_canfar_v2.0.bash -c $config_dir -p $psf -r $retrieve --tile_det $tile_det --tile_mask $tile_mask -j $exp_job --n_smp $N_SMP --nsh_jobs $N_SMP $debug_flag" $dry_run 2>&1 | tee -a "$exp_log_file" + [ "$quiet" == "0" ] && echo "Done with job_sp_canfar_v2.0.bash" + exec {lock_fd}>&- done < "$exp_numbers_file" @@ -343,20 +408,12 @@ function run_tile_job() { local complete_checks=$3 local main_prefix="${run_prefixes%% *}" - # force: remove all existing run directories for each prefix before running - if [ "$force" == "1" ]; then - local run_prefix - for run_prefix in $run_prefixes; do - local dirs_to_remove - dirs_to_remove=$(ls -d "$work_dir/output/run_sp_tile_${run_prefix}"* 2>/dev/null) - if [ -n "$dirs_to_remove" ]; then - for d in $dirs_to_remove; do - message "Force-removing $d" "$debug_out" -1 - command "rm -rf $d" $dry_run - done - fi - done - fi + # Tile jobs must run in the tile work dir: SP_RUN/SP_CONFIG derive from cwd, + # and run_exp_job leaves cwd in the last exposure dir. When run_job is invoked + # once for the whole pipeline (-j 4095, as the sp_validation workflow does), + # the tile jobs after the exposure loop would otherwise inherit that stale + # exposure cwd and write into it. Re-anchor to the tile dir here. + cd "$work_dir" # Locate most recent existing run directory for the main prefix local run_dir @@ -421,16 +478,24 @@ function run_tile_job() { done fi - if [ "$is_complete" == "1" ] && [ -n "$complete_checks" ]; then + # force: always remove and rerun. + # retry: remove and rerun only if the existing run is incomplete/failed. + if [ "$force" == "1" ]; then + remove_run_dirs "$work_dir" "tile" "$run_prefixes" "Force-removing" + elif [ "$retry" == "1" ] && [ "$is_complete" == "0" ]; then + remove_run_dirs "$work_dir" "tile" "$run_prefixes" "Retry: removing incomplete run" + fi + + if [ "$force" != "1" ] && [ "$is_complete" == "1" ] && [ -n "$complete_checks" ]; then message "Complete: ( $check_desc)" "$debug_out" -1 return 0 fi - if [ "$check" == "1" ]; then + if [ "$check_only" == "1" ]; then if [ -n "$run_dir" ]; then - message "Incomplete: ($check_desc)" "$debug_out" -1 + message "Incomplete: ( $check_desc)" "$debug_out" -1 else - message "Missing: ($check_desc)" "$debug_out" -1 + message "Missing: ( $check_desc)" "$debug_out" -1 fi return 0 fi @@ -440,15 +505,85 @@ function run_tile_job() { [ -n "$debug_out" ] && debug_flag="--debug_out $debug_out" if [ ! -e "cfis" ]; then - ln -sf ~/shapepipe/example/cfis "cfis" + ln -sf $config_dir cfis fi command "update_runs_log_file.py" $dry_run # Run job script - command "job_sp_canfar_v2.0.bash -p $psf --tile_det $tile_det --tile_mask $tile_mask -j $tile_job --n_smp $N_SMP --nsh_jobs $N_SMP $debug_flag" $dry_run 2>&1 | tee -a "$log_file" + command "$SP_ROOT/scripts/sh/job_sp_canfar_v2.0.bash -c $config_dir -p $psf -r $retrieve --tile_det $tile_det --tile_mask $tile_mask -j $tile_job --n_smp $N_SMP --nsh_jobs $N_SMP $debug_flag" $dry_run 2>&1 | tee -a "$log_file" } + +# Wrap a run_tile_job or run_exp_job call with check-log and force-cleanup. +# Usage: run_job_logged BIT FUNC [FUNC_ARGS...] +# - With --force: removes log_job_BIT.txt before running FUNC +# - With --check: captures output of FUNC; writes log_job_BIT.txt only if no +# Missing/Incomplete lines appear +# - Otherwise: calls FUNC directly +function run_job_logged() { + local bit=$1 + shift + local log="$dir/logs/log_job_${ID}_${bit}.txt" + + [ "$force" == "1" ] && rm -f "$log" + + if [ "$check_only" == "1" ]; then + local out + out=$( "$@" 2>&1 ) + echo "$out" + echo "$out" | grep -qiE "(Missing|Incomplete)" || echo "$out" > "$log" + else + "$@" + local rc=$? + if [ $rc -eq 0 ]; then + if [ "$run_and_check" == "1" ]; then + check_only=1 + local saved_force=$force + force=0 + local out + out=$( "$@" 2>&1 ) + check_only=0 + force=$saved_force + echo "$out" + echo "$out" | grep -qiE "(Missing|Incomplete)" || echo "$out" > "$log" + else + echo "Completed job $bit $(date)" > "$log" + fi + fi + fi +} + + +if [ "$type" == "data" ]; then + + [ "$quiet" == "0" ] && echo "Running on data" + retrieve="vos" + config_dir="$SP_ROOT/example/cfis" + export SP_DIR=$dir + export SP_CONFIG=$config_dir + +elif [ "$type" == "image_sims" ]; then + + [ "$quiet" == "0" ] && echo "Running on image simulations" + retrieve="symlink" + config_dir="$SP_ROOT/example/cfis_image_sims" + # SP_DIR points to the run directory where input_tiles and input_exp live; + # configs use $SP_DIR/input_* so those dirs stay outside SP_RUN and are + # not found twice by ShapePipe's recursive glob scan. + export SP_DIR=$dir + export SP_CONFIG=$config_dir + tile_det='sx' + +else + + echo "Invalid input type $type" + +fi + +[ "$quiet" == "0" ] && echo "config_dir=$config_dir" + + # Init message if [ "$test_only" == "1" ]; then message "$(basename "$0") test mode, exiting." "$debug_out" 0 @@ -469,17 +604,28 @@ if [ "$ID" == "-1" ]; then fi if [ "$psf" != "psfex" ] && [ "$psf" != "mccd" ]; then - message "PSF (option -p) needs to be 'psfex' or 'mccd'" "$debug_out" 4 + message "PSF (option -p) needs to be 'psfex' or 'mccd', not '$psf'" "$debug_out" 4 fi if [ "$dry_run" != "0" ] && [ "$dry_run" != "1" ]; then message "dry_run must be 0 or 1, not $dry_run" "$debug_out" 8 fi +## Check input links +for link in "$dir/input_tiles" "$dir/input_exp" "$dir/cfis"; do + if [ -L "$link" ]; then + if [ ! -e "$link" ]; then + message "Broken symlink: $link" "$debug_out" 6 + fi + elif [ ! -e "$link" ]; then + message "Missing path: $link" "$debug_out" 6 + fi +done + # Start script -source $HOME/shapepipe/scripts/sh/functions.sh +source "$SP_ROOT/scripts/sh/functions.sh" message "Starting $(basename "$0")" "$debug_out" -1 message "`date`" "$debug_out" -1 @@ -495,18 +641,28 @@ Letter=${letter^} cd $dir -# Derive tile path components from ID (e.g. "000.227" -> IDra="000") +# Derive tile path components from ID (e.g. "000.227" -> IDra="000", IDdec="227") IDra=${ID%%.*} +IDdec=${ID##*.} +ID_DASHED="${IDra}-${IDdec}" work_dir="$dir/tiles/$IDra/$ID" log_file="$work_dir/job_sp_canfar_v2.0.log" # Create tile work directory [ ! -d "$work_dir" ] && command "mkdir -p $work_dir" $dry_run cd $work_dir +echo "$0 $@" > "$log_file" # Write ID to first input +# Image sims use dash format (e.g. 233-293); real data uses dot format (233.293) +# which ShapePipe's in2out_pattern converts to dashes for output naming only, +# not for input file lookup — so write the format that matches the actual files. if [ ! -e tile_numbers.txt ]; then - echo $ID > tile_numbers.txt + if [ "$type" == "image_sims" ]; then + echo ${ID//./-} > tile_numbers.txt + else + echo $ID > tile_numbers.txt + fi fi # Output directory @@ -514,7 +670,8 @@ if [ ! -d "output" ]; then command "mkdir output" $dry_run fi -echo -n "pwd: "; pwd + +[ "$quiet" == "0" ] && { echo -n "pwd: "; pwd; } # Avoid Qt error with setools @@ -537,44 +694,99 @@ IDdec=${ID##*.} (( do_job = job & 1 )) if [[ $do_job != 0 ]]; then # Job 1: download tile images and weights - run_tile_job 1 "Git" "get_images_runner:4" + if [ "$type" == "image_sims" ]; then + n_exp=2 + else + n_exp=4 + fi + run_job_logged 1 run_tile_job 1 "Git" "get_images_runner:${n_exp}" fi (( do_job = job & 2 )) if [[ $do_job != 0 ]]; then - # Job 2: uncompress tile weights - run_tile_job 2 "Uz" "uncompress_fits_runner:1" + log_2="$dir/logs/log_job_${ID}_2.txt" + [ "$force" == "1" ] && rm -f "$log_2" + if [ "$type" == "image_sims" ]; then + # Image sims weights are already uncompressed; fake the Uz output directory + # so downstream jobs can find the weight via last:uncompress_fits_runner. + weight_src="$dir/input_tiles/CFIS_simu_weight-${ID//./-}.fits" + if [ "$check_only" == "1" ]; then + uz_run_dir=$(ls -dt "$work_dir/output/run_sp_tile_Uz"* 2>/dev/null | head -1) + if [ -n "$uz_run_dir" ] && [ -e "$uz_run_dir/uncompress_fits_runner/output/$(basename $weight_src)" ]; then + msg="Complete: ( Uz/uncompress_fits_runner[fake] 1/1 )" + message "$msg" "$debug_out" -1 + echo "$msg" > "$log_2" + else + message "Missing: Uz $(basename $weight_src)" "$debug_out" -1 + fi + else + uz_out="$work_dir/output/run_sp_tile_Uz$(date +_%Y-%m-%d_%H-%M-%S)/uncompress_fits_runner/output" + command "mkdir -p $uz_out" $dry_run + if [ -e "$weight_src" ] && [ ! -e "$uz_out/$(basename $weight_src)" ]; then + command "ln -sf $weight_src $uz_out/$(basename $weight_src)" $dry_run + fi + if [ "$run_and_check" == "1" ]; then + uz_run_dir=$(ls -dt "$work_dir/output/run_sp_tile_Uz"* 2>/dev/null | head -1) + if [ -n "$uz_run_dir" ] && [ -e "$uz_run_dir/uncompress_fits_runner/output/$(basename $weight_src)" ]; then + msg="Complete: Uz $(basename $weight_src)" + message "$msg" "$debug_out" -1 + echo "$msg" > "$log_2" + else + message "Missing: Uz $(basename $weight_src)" "$debug_out" -1 + fi + else + echo "Completed job 2 $(date)" > "$log_2" + fi + fi + else + # Job 2: uncompress tile weights + run_job_logged 2 run_tile_job 2 "Uz" "uncompress_fits_runner:1" + fi fi (( do_job = job & 4 )) if [[ $do_job != 0 ]]; then # Job 4: find exposures - run_tile_job 4 "Fe" "find_exposures_runner:1" + run_job_logged 4 run_tile_job 4 "Fe" "find_exposures_runner:1" fi (( do_job = job & 8 )) if [[ $do_job != 0 ]]; then # Job 8: retrieve exposure images - run_exp_job 8 "Gie" "get_images_runner:6" + if [ "$type" == "image_sims" ]; then + n_exp=3 + else + n_exp=6 + fi + run_job_logged 8 run_exp_job 8 "Gie" "get_images_runner:${n_exp}" fi (( do_job = job & 16 )) if [[ $do_job != 0 ]]; then # Job 16: split exposures, get WCS headers - run_exp_job 16 "Sp" "split_exp_runner:121" + run_job_logged 16 run_exp_job 16 "Sp" "split_exp_runner:121" fi (( do_job = job & 32 )) if [[ $do_job != 0 ]]; then # Job 32: mask exposures - run_exp_job 32 "Ma" "mask_runner:40" + run_job_logged 32 run_exp_job 32 "Ma" "mask_runner:40" fi (( do_job = job & 64 )) if [[ $do_job != 0 ]]; then - # Job 64: process stars on exposures, PSF model - if [ "$psf" == "psfex" ]; then - run_exp_job 64 "SxSePsf${Letter}i" "sextractor_runner:80 psfex_runner:80 psfex_interp_runner:40::warn setools_runner:80:rand_split" + # Job 64: PSF model + # For image_sims: fake PSF runs as part of job 512 (requires sexcat from job 256) + # For data: run full exposure-level PSF modelling pipeline + if [ "$type" == "image_sims" ]; then + # Fake PSF is handled inside job 512; write placeholder log so the sequence is complete + log_64="$dir/logs/log_job_${ID}_64.txt" + [ "$force" == "1" ] && rm -f "$log_64" + msg="Complete: job 64 placeholder (fake PSF runs as part of job 512)" + message "$msg" "$debug_out" -1 + echo "$msg" > "$log_64" + elif [ "$psf" == "psfex" ]; then + run_job_logged 64 run_exp_job 64 "SxSePsf${Letter}i" "sextractor_runner:80 psfex_runner:80 psfex_interp_runner:40::warn setools_runner:80:rand_split" else message "MCCD not implemented yet for v2.0" "$debug_out" 10 fi @@ -583,31 +795,81 @@ fi (( do_job = job & 128 )) if [[ $do_job != 0 ]]; then # Job 128: merge exposure WCS headers into tile-level sqlite log - run_tile_job 128 "Mh_exp" "merge_headers_runner:1" + run_job_logged 128 run_tile_job 128 "Mh_exp" "merge_headers_runner:1" fi (( do_job = job & 256 )) if [[ $do_job != 0 ]]; then # Job 256: object selection on tiles if [ "$tile_det" == "uc" ]; then - run_tile_job 256 "Gic Uc" "get_images_runner:2 read_ext_sexcat_runner:1" + run_job_logged 256 run_tile_job 256 "Gic Uc" "get_images_runner:2 read_ext_sexcat_runner:1" else - run_tile_job 256 "Sx" "sextractor_runner:1" + n_exp=2 + run_job_logged 256 run_tile_job 256 "Sx" "sextractor_runner:$n_exp" fi fi (( do_job = job & 512 )) if [[ $do_job != 0 ]]; then - # Job 512: process tiles (PSF interp, vignet) - run_tile_job 512 "${Letter}iViVi ${Letter}iViVi ${Letter}iViVi" "psfex_interp_runner:1 vignetmaker_runner_run_1:1 vignetmaker_runner_run_2:4" + # Job 512: process tiles ([PSF interp,] vignets) + # For image_sims: fake PSF runs first (requires sexcat from job 256), then vignets + log_512="$dir/logs/log_job_${ID}_512.txt" + [ "$force" == "1" ] && rm -f "$log_512" + if [ "$type" == "data" ]; then + if [ "$check_only" == "1" ]; then + out=$(run_tile_job 512 "${Letter}iViVi ${Letter}iViVi ${Letter}iViVi" "psfex_interp_runner:1 vignetmaker_runner_run_1:1 vignetmaker_runner_run_2:4" 2>&1) + echo "$out" + echo "$out" | grep -qiE "(Missing|Incomplete)" || echo "$out" > "$log_512" + else + run_tile_job 512 "${Letter}iViVi ${Letter}iViVi ${Letter}iViVi" "psfex_interp_runner:1 vignetmaker_runner_run_1:1 vignetmaker_runner_run_2:4" + if [ "$run_and_check" == "1" ]; then + check_only=1; local saved_force_512=$force; force=0 + out=$(run_tile_job 512 "${Letter}iViVi ${Letter}iViVi ${Letter}iViVi" "psfex_interp_runner:1 vignetmaker_runner_run_1:1 vignetmaker_runner_run_2:4" 2>&1) + check_only=0; force=$saved_force_512 + echo "$out" + echo "$out" | grep -qiE "(Missing|Incomplete)" || echo "$out" > "$log_512" + else + echo "Completed job 512 $(date)" > "$log_512" + fi + fi + else + if [ "$check_only" == "1" ]; then + out1=$(run_tile_job 64 "fpsf" "fake_psf_runner:1" 2>&1) + out2=$(run_tile_job 512 "ViVi ViVi" "vignetmaker_runner_run_1:1 vignetmaker_runner_run_2:3" 2>&1) + echo "$out1"; echo "$out2" + { echo "$out1"; echo "$out2"; } | grep -qiE "(Missing|Incomplete)" || \ + { echo "$out1"; echo "$out2"; } > "$log_512" + else + run_tile_job 64 "fpsf" "fake_psf_runner:1" + run_tile_job 512 "ViVi ViVi" "vignetmaker_runner_run_1:1 vignetmaker_runner_run_2:3" + if [ "$run_and_check" == "1" ]; then + check_only=1; saved_force_512im=$force; force=0 + out1=$(run_tile_job 64 "fpsf" "fake_psf_runner:1" 2>&1) + out2=$(run_tile_job 512 "ViVi ViVi" "vignetmaker_runner_run_1:1 vignetmaker_runner_run_2:3" 2>&1) + check_only=0; force=$saved_force_512im + echo "$out1"; echo "$out2" + { echo "$out1"; echo "$out2"; } | grep -qiE "(Missing|Incomplete)" || \ + { echo "$out1"; echo "$out2"; } > "$log_512" + else + echo "Completed job 512 $(date)" > "$log_512" + fi + fi + fi fi (( do_job = job & 1024 )) if [[ $do_job != 0 ]]; then # Job 1024: shape measurement - run_tile_job 1024 "Ng" "ngmix_interp_runner:1" + run_job_logged 1024 run_tile_job 1024 "Ng" "ngmix_runner:1" fi +(( do_job = job & 2048 )) +if [[ $do_job != 0 ]]; then + # Job 2048: merge catalogues + run_job_logged 2048 run_tile_job 2048 "Mc_${psf}" "make_cat_runner:1" +fi + + if [ -n "$scratch" ]; then message "Syncing output from scratch back to permanent dir" "$debug_out" -1 command "rsync -a output/ $work_dir/output/" $dry_run diff --git a/src/shapepipe/modules/fake_psf_package/__init__.py b/src/shapepipe/modules/fake_psf_package/__init__.py new file mode 100644 index 000000000..c77dd5658 --- /dev/null +++ b/src/shapepipe/modules/fake_psf_package/__init__.py @@ -0,0 +1,28 @@ +"""FAKE PSF PACKAGE. + +This package creates fake PSF postage-stamp dictionaries for image simulations. + +:Author: Martin Kilbinger + +:Parent module: ``sextractor_runner`` + +:Input: Tile SExtractor catalogue (multi-epoch FITS format) + +:Output: Per-galaxy PSF sqlite dictionary (``galaxy_psf-XXX-XXX.sqlite``) + +Description +=========== + +For each galaxy in the tile catalogue the module looks up the pre-computed +PSF stamp for every contributing exposure-CCD combination from a pickled +PSF dictionary and writes the result to a SqliteDict file in the same format +produced by ``psfex_interp_runner``. + +Module-specific config file entries +==================================== + +PSF_DICT_PATH : str + Path to the pickled PSF dictionary (``Full_psf_dict.pickle``). +""" + +__all__ = ["fake_psf.py"] diff --git a/src/shapepipe/modules/fake_psf_package/fake_psf.py b/src/shapepipe/modules/fake_psf_package/fake_psf.py new file mode 100644 index 000000000..1313b0917 --- /dev/null +++ b/src/shapepipe/modules/fake_psf_package/fake_psf.py @@ -0,0 +1,112 @@ +"""FAKE PSF. + +Create fake PSF postage-stamp dictionaries for image simulations. + +:Author: Martin Kilbinger + +""" + +import pickle + +import numpy as np +from astropy.io import fits +from sqlitedict import SqliteDict + + +class FakePsf: + """Fake PSF. + + Parameters + ---------- + sexcat_path : str + Path to the tile SExtractor catalogue (multi-epoch FITS format). + psf_dict_path : str + Path to the pickled PSF dictionary. + output_path : str + Path for the output SqliteDict file. + w_log : logging.Logger + Pipeline logger. + """ + + def __init__(self, sexcat_path, psf_dict_path, output_path, w_log): + self._sexcat_path = sexcat_path + self._psf_dict_path = psf_dict_path + self._output_path = output_path + self._w_log = w_log + + def process(self): + """Run fake PSF creation.""" + self._w_log.info(f"Reading sexcat: {self._sexcat_path}") + try: + sex = fits.open(self._sexcat_path, ignore_missing_simple=True) + except Exception as e: + self._w_log.error(f"Error opening FITS file {self._sexcat_path}: {e}") + raise + + self._w_log.info(f"Catalogue has {len(sex)} HDUs") + + if len(sex) < 4: + self._w_log.error( + f"Catalogue file has insufficient HDUs: expected ≥4, got {len(sex)} " + f"(file: {self._sexcat_path})" + ) + raise ValueError(f"Invalid catalogue: only {len(sex)} HDUs (need ≥4)") + + self._w_log.info(f"Loading PSF dictionary: {self._psf_dict_path}") + try: + with open(self._psf_dict_path, "rb") as f: + psf_dict = pickle.load(f) + except Exception as e: + self._w_log.error(f"Error loading PSF dictionary {self._psf_dict_path}: {e}") + raise + + try: + n_gal = len(sex[3].data.field("NUMBER")) + except Exception as e: + self._w_log.error( + f"Error reading catalogue data from HDU 3 in {self._sexcat_path}: {e}" + ) + raise + n_exp = len(sex) - 3 + + self._w_log.info(f"Processing {n_gal} galaxies over {n_exp} epochs") + + # Build (n_gal, n_exp) array of "expname-ccdnum" strings + string_array = np.empty((n_gal, n_exp), dtype="object") + for exp_ind in range(3, len(sex)): + ccd_n = sex[exp_ind].data.field("CCD_N") + exp_name = sex[exp_ind].data.field("EXP_NAME") + string_array[:, exp_ind - 3] = [ + f"{int(exp)}-{int(ccd)}" + for exp, ccd in zip(exp_name, ccd_n) + ] + + # Mask invalid entries (exp or ccd == -99 produces strings with "--") + mask = np.array( + [["--" in item for item in row] for row in string_array] + ) + masked = np.ma.masked_array(string_array, mask) + + # Build per-galaxy PSF dictionaries + output_file = SqliteDict(self._output_path) + missing = 0 + for idx, gal_row in enumerate(masked): + galaxy_number = idx + 1 # 1-based, matches NUMBER field + gal_dict = {} + for exp_ccd in gal_row.compressed(): + if exp_ccd not in psf_dict: + missing += 1 + self._w_log.warning( + f"Galaxy {galaxy_number}: key '{exp_ccd}' not in PSF dict" + ) + continue + gal_dict[exp_ccd] = psf_dict[exp_ccd] + output_file[str(galaxy_number)] = gal_dict + + output_file.commit() + output_file.close() + sex.close() + + if missing: + self._w_log.warning(f"{missing} missing PSF entries across all galaxies") + self._w_log.info(f"Written: {self._output_path}") diff --git a/src/shapepipe/modules/fake_psf_runner.py b/src/shapepipe/modules/fake_psf_runner.py new file mode 100644 index 000000000..d3dcbbc0e --- /dev/null +++ b/src/shapepipe/modules/fake_psf_runner.py @@ -0,0 +1,36 @@ +"""FAKE PSF RUNNER. + +Module runner for ``fake_psf``. + +:Author: Martin Kilbinger + +""" + +from shapepipe.modules.fake_psf_package import fake_psf +from shapepipe.modules.module_decorator import module_runner + + +@module_runner( + version="1.0", + file_pattern=["sexcat"], + file_ext=".fits", + depends=["numpy", "astropy", "sqlitedict"], + numbering_scheme="-000-000", +) +def fake_psf_runner( + input_file_list, + run_dirs, + file_number_string, + config, + module_config_sec, + w_log, +): + """Define The Fake PSF Runner.""" + sexcat_path = input_file_list[0] + psf_dict_path = config.getexpanded(module_config_sec, "PSF_DICT_PATH") + output_path = f'{run_dirs["output"]}/galaxy_psf{file_number_string}.sqlite' + + inst = fake_psf.FakePsf(sexcat_path, psf_dict_path, output_path, w_log) + inst.process() + + return None, None diff --git a/src/shapepipe/modules/mask_package/mask.py b/src/shapepipe/modules/mask_package/mask.py index 17e052a6e..144a4ea34 100644 --- a/src/shapepipe/modules/mask_package/mask.py +++ b/src/shapepipe/modules/mask_package/mask.py @@ -429,10 +429,14 @@ def make_mask(self): ) # Handle stdout / stderr - general_stdout = f"\nCDSClient\n{self._CDS_stdout}" + # _CDS_stdout/_CDS_stderr are only set when find_stars ran, i.e. + # when HALO_MAKE or SPIKE_MAKE is True (False for image sims) + general_stdout = "" general_stderr = "" - if self._CDS_stderr != "": - general_stderr += f"\nCDSClient\n{self._CDS_stderr}" + if hasattr(self, "_CDS_stdout"): + general_stdout += f"\nCDSClient\n{self._CDS_stdout}" + if self._CDS_stderr != "": + general_stderr += f"\nCDSClient\n{self._CDS_stderr}" if hasattr(self, "_WW_stdout") or hasattr(self, "_WW_stdout"): general_stdout += f"\n\nWeightWatcher\n{self._WW_stdout}" if self._WW_stderr != "": diff --git a/src/shapepipe/modules/merge_headers_runner.py b/src/shapepipe/modules/merge_headers_runner.py index 9bee35159..6c3cfa38e 100644 --- a/src/shapepipe/modules/merge_headers_runner.py +++ b/src/shapepipe/modules/merge_headers_runner.py @@ -43,6 +43,9 @@ def merge_headers_runner( f"Tile-level merge: collecting headers from {exp_base_dir} " f"using {exp_numbers_file}" ) + warn_missing = config.getboolean( + module_config_sec, "WARN_MISSING_EXP", fallback=False + ) headers_file_list = get_exp_output_files( exp_base_dir, exp_numbers_file, @@ -50,6 +53,7 @@ def merge_headers_runner( "headers", ".npy", w_log=w_log, + warn_only=warn_missing, ) # Extract tile number from the exp_numbers filename, e.g. # "exp_numbers-284.272-1.000.txt" -> "-284.272-1.000" diff --git a/src/shapepipe/modules/ngmix_package/ngmix.py b/src/shapepipe/modules/ngmix_package/ngmix.py index 371264e23..af0eecd96 100644 --- a/src/shapepipe/modules/ngmix_package/ngmix.py +++ b/src/shapepipe/modules/ngmix_package/ngmix.py @@ -119,6 +119,64 @@ def get_prior(pixel_scale, rng, T_range=None, F_range=None): ) +def position_seed(ra, dec, ccd): + """Deterministic RNG seed from an object's sky position (ngmix#796). + + For image-simulation m-bias with the Pujol estimator, the same scene is + simulated under different applied shears ("image branches") and the shear + response is read from the branch difference of the SAME objects. Metacal's + ``fixnoise`` adds a counter-noise realisation drawn from an RNG; if that RNG + is seeded per tile, an object gets *different* added noise in different + branches (detection order differs), and the noise fails to cancel in the + difference, inflating ``sigma_m``. Seeding the per-object RNG from sky + position instead makes the same object draw the same added noise (and the + same fit guesses) in every branch, so both cancel and the m-bias error + shrinks. Off in production; a knob for the sim path only. + + Box math (kept exactly as Fabian's issue #796):: + + box_x = floor(ra * 3600 / 3) + (ccd + 1) + box_y = floor(dec * 3600 / 3) + (ccd + 2) + + The 3-arcsec boxes make the seed robust to detection-order changes and to + small centroid jitter within a box: the same physical object lands in the + same box across branches and so gets the same seed. The ``+ccd`` offsets + disambiguate exposure overlap. + + **Deviation from #796.** Fabian combines the boxes as ``box_x + box_y``, + which collides along anti-diagonals — e.g. boxes ``(10, 20)`` and + ``(11, 19)`` both give seed 30, so two distinct objects share a noise + stream. We instead combine them with a Cantor pairing (a bijection on + non-negative integers), after a zig-zag fold that maps signed box indices + (dec can be negative) onto non-negative ones. The result is reduced mod + ``2**32`` to land in the valid ``numpy.random.RandomState`` seed range + ``[0, 2**32)``. The box math is untouched; only the collision-prone sum is + replaced. + + Parameters + ---------- + ra, dec : float + Object sky position in degrees (first-epoch coordinates; all epochs of + one object share them). + ccd : int + CCD number of the first epoch. + + Returns + ------- + int + Seed in ``[0, 2**32)`` for ``numpy.random.RandomState``. + """ + box_x = int(np.floor((ra * 3600) / 3) + (ccd + 1)) + box_y = int(np.floor((dec * 3600) / 3) + (ccd + 2)) + # Zig-zag fold signed -> non-negative (0,-1,1,-2,... -> 0,1,2,3,...) so the + # Cantor pairing, which is a bijection on the non-negative integers, stays + # injective for southern (dec<0) boxes. + zx = 2 * box_x if box_x >= 0 else -2 * box_x - 1 + zy = 2 * box_y if box_y >= 0 else -2 * box_y - 1 + cantor = (zx + zy) * (zx + zy + 1) // 2 + zy + return cantor % (2 ** 32) + + class Tile_cat(): """Tile_cat. @@ -186,6 +244,10 @@ def __init__( self.wcs = [] self.ra = [] self.dec = [] + # CCD number of the first epoch, used only to build the per-object + # position seed (see :func:`position_seed`, ngmix#796). ``None`` when + # position seeding is off, so the field costs nothing on the hot path. + self.ccd = None self.bkg_sub = bkg_sub self.megacam_flip = megacam_flip @@ -213,11 +275,10 @@ def __init__( flag_vignet_path, f_wcs_path, bkg_rms_vignet_path=None, - ): self.f_wcs_file = SqliteDict(f_wcs_path) self.gal_vign_cat = SqliteDict(gal_vignet_path) - self.bkg_vign_cat = SqliteDict(bkg_vignet_path) + self.bkg_vign_cat = SqliteDict(bkg_vignet_path) if bkg_vignet_path is not None else None self.psf_vign_cat = SqliteDict(psf_vignet_path) self.weight_vign_cat = SqliteDict(weight_vignet_path) self.flag_vign_cat = SqliteDict(flag_vignet_path) @@ -230,7 +291,8 @@ def __init__( def close(self): self.f_wcs_file.close() self.gal_vign_cat.close() - self.bkg_vign_cat.close() + if self.bkg_vign_cat is not None: + self.bkg_vign_cat.close() self.flag_vign_cat.close() self.weight_vign_cat.close() self.psf_vign_cat.close() @@ -273,6 +335,13 @@ class Ngmix(object): (robust for galaxies); ``"wcs"`` uses the catalog sky position projected through the WCS (better for stars, whose HSM moments are noisy). See :func:`make_ngmix_observation`. + seed_from_position : bool, optional + If ``True``, replace the tile-level RNG with a per-object RNG seeded + from the object's sky position (:func:`position_seed`) inside the + object loop, so metacal's ``fixnoise`` counter-noise and the fit + guesses cancel across Pujol image-simulation branches (ngmix#796). The + default ``False`` leaves the production path byte-identical. See + :func:`position_seed` for the physics and the seed construction. Raises ------ @@ -293,24 +362,46 @@ def __init__( save_batch=-1, id_obj_min=-1, id_obj_max=-1, + bkg_sub=True, centroid_source="hsm", + seed_from_position=False, ): - if len(input_file_list) not in {6, 7}: + # Base count = catalogue + vignets, excluding the f_wcs headers (passed + # separately). One fewer when the background vignet is absent + # (``bkg_sub=False``, image sims). An extra trailing slot may carry the + # optional background-rms vignet (#779), so both counts are valid. + n_base = 6 if bkg_sub else 5 + if len(input_file_list) not in {n_base, n_base + 1}: raise IndexError( f"Input file list has length {len(input_file_list)}," - + " required is 6 or 7" + + f" required is {n_base} or {n_base + 1}" ) self._tile_cat_path = input_file_list[0] + if bkg_sub: + bkg_path, psf_path, weight_path, flag_path = ( + input_file_list[2], input_file_list[3], + input_file_list[4], input_file_list[5], + ) + else: + bkg_path, psf_path, weight_path, flag_path = ( + None, input_file_list[2], + input_file_list[3], input_file_list[4], + ) + bkg_rms_vignet_path = ( + input_file_list[n_base] + if len(input_file_list) == n_base + 1 + else None + ) self._vignet_cat = Vignet( input_file_list[1], - input_file_list[2], - input_file_list[3], - input_file_list[4], - input_file_list[5], + bkg_path, + psf_path, + weight_path, + flag_path, f_wcs_path, - input_file_list[6] if len(input_file_list) == 7 else None, + bkg_rms_vignet_path, ) self._output_dir = output_dir @@ -324,7 +415,9 @@ def __init__( self._save_batch = save_batch self._id_obj_min = id_obj_min self._id_obj_max = id_obj_max + self._bkg_sub = bkg_sub self._centroid_source = centroid_source + self._seed_from_position = seed_from_position self._w_log = w_log @@ -332,6 +425,11 @@ def __init__( seed = int(''.join(re.findall(r'\d+', self._file_number_string))) self._rng = np.random.RandomState(seed) self._w_log.info(f'Random generator initialisation seed = {seed}') + if self._seed_from_position: + self._w_log.info( + 'SEED_FROM_POSITION on: per-object RNG seeded from sky position' + ' for Pujol noise cancellation (image sims, ngmix#796)' + ) @classmethod def MegaCamFlip(self, vign, ccd_nb): @@ -477,9 +575,13 @@ def compile_results(self, results): ) # ngmix 2.x reports the solver's function-evaluation count # (nfev, ~tens-hundreds; -1 on some failures), not the v1 - # 1-5 retry count, so the column is named accordingly. + # 1-5 retry count, so the column is named accordingly. Fits + # that fail before nfev is ever set (e.g. the NaN-filled + # branch above) fall back to the same -1 sentinel, keeping + # the column int64 (FITS 'K') across every batch: nfev_fit + # is int64 with -1 meaning failed/absent. output_dict[name]["nfev_fit"].append( - fit.get("nfev", np.nan) + fit.get("nfev", -1) ) # The two PSF families are object-level (one value per # object, not per shear type) and self-named: every key @@ -629,6 +731,61 @@ def check_key(self, expccd_name_tmp, vign_cat, vignet_path): + f" file '{vignet_path}'" ) + def log_mean_ellipticity(self): + """Log mean ellipticity from NOSHEAR HDU to the run log. + + Reports , with standard errors for all objects and for + objects passing the default metacal cuts (flags==0, mcal_flags==0, + 10 < SNR < 500, T/Tpsf > 0.5). + """ + output_path = self.get_output_path(self._output_dir) + try: + with fits.open(output_path) as hdul: + d = hdul['NOSHEAR'].data + g1 = d['g1'].astype(float) + g2 = d['g2'].astype(float) + flags = d['flags'] + mcal_flags = d['mcal_flags'] + s2n = d['s2n'].astype(float) + T = d['T'].astype(float) + Tpsf = d['Tpsf'].astype(float) + except Exception as e: + self._w_log.warning(f"Could not compute mean ellipticity: {e}") + return + + n_total = len(g1) + if n_total == 0: + self._w_log.info("Mean ellipticity: no objects in output catalogue") + return + + def _log_stats(g1_sel, g2_sel, label): + n = len(g1_sel) + if n == 0: + self._w_log.info(f" {label}: 0 objects") + return + mean_g1 = g1_sel.mean() + mean_g2 = g2_sel.mean() + err_g1 = g1_sel.std() / np.sqrt(n) + err_g2 = g2_sel.std() / np.sqrt(n) + self._w_log.info( + f" {label} (N={n}):" + f" = {mean_g1:+.4e} +/- {err_g1:.4e}," + f" = {mean_g2:+.4e} +/- {err_g2:.4e}" + ) + + self._w_log.info(f"Mean ellipticity (NOSHEAR, N_total={n_total}):") + _log_stats(g1, g2, "no cuts") + + with np.errstate(invalid='ignore'): + mask = ( + (flags == 0) + & (mcal_flags == 0) + & (s2n >= 10.0) + & (s2n <= 500.0) + & (T / Tpsf >= 0.5) + ) + _log_stats(g1[mask], g2[mask], "SNR in [10, 500], T/Tpsf > 0.5") + def process(self): """Process. @@ -678,12 +835,31 @@ def process(self): n_empty_cat += 1 continue - stamp = prepare_postage_stamps(vignet_cat, obj_id, i_tile, tile_cat) + stamp = prepare_postage_stamps(vignet_cat, obj_id, i_tile, tile_cat, self._bkg_sub) if len(stamp.gals) == 0: n_no_epoch += 1 continue + # Position-seeded per-object RNG for Pujol noise cancellation in + # image sims (ngmix#796): the same object gets the same fixnoise + # counter-noise and fit guesses in every shear branch, so both + # cancel in the branch difference. The prior is rebuilt from the + # same per-object RNG because the guesser draws its initial guess + # via prior.sample() (ngmix guessers.py), which consumes the RNG the + # prior was CONSTRUCTED with — so a per-object rng alone would leave + # the guess drawing from the shared tile stream and break + # cancellation. Off in production, where the single tile-level + # self._rng and the tile-level prior carry the whole loop. + if self._seed_from_position: + obj_rng = np.random.RandomState( + position_seed(stamp.ra[0], stamp.dec[0], stamp.ccd) + ) + obj_prior = get_prior(self._pixel_scale, obj_rng) + else: + obj_rng = self._rng + obj_prior = prior + try: flux_guess = ( tile_cat.flux[i_tile] @@ -692,9 +868,9 @@ def process(self): ) res, psf_res, psf_orig_res = do_ngmix_metacal( stamp, - prior, + obj_prior, flux_guess, - self._rng, + obj_rng, centroid_source=self._centroid_source, ) except Exception as ee: @@ -769,9 +945,12 @@ def process(self): # Save results self.save_results(res_dict) -def prepare_postage_stamps(vignet, obj_id, i_tile, tile_cat): + # Log mean ellipticity statistics + self.log_mean_ellipticity() + +def prepare_postage_stamps(vignet, obj_id, i_tile, tile_cat, bkg_sub=True): # define per-object lists of individual exposures to go into ngmix - stamp = Postage_stamp() + stamp = Postage_stamp(bkg_sub=bkg_sub) #identify exposure and ccd number from psf catalog psf_expccd_names = list(vignet.psf_vign_cat[str(obj_id)].keys()) for expccd_name in psf_expccd_names: @@ -795,6 +974,13 @@ def prepare_postage_stamps(vignet, obj_id, i_tile, tile_cat): else: gal_vign_sub_bkg = gal_vign + # Skip epochs where sigma_mad=0 (CCD edge: mostly-zero stamp). + # prepare_ngmix_weights divides by sig_noise; zero sigma causes NaN/inf + # weights that corrupt GalSim's C-level FFT allocations. + # MKDEBUG TODO this needs to be checked! + if not sigma_mad(gal_vign_sub_bkg) > 0: + continue + tile_vign = ( np.copy(tile_cat.vign[i_tile]) if tile_cat.vign is not None @@ -856,6 +1042,11 @@ def prepare_postage_stamps(vignet, obj_id, i_tile, tile_cat): stamp.wcs.append(epoch_wcs) stamp.ra.append(tile_cat.ra[i_tile]) stamp.dec.append(tile_cat.dec[i_tile]) + # CCD of the first surviving epoch — Fabian's coord_list[0] convention + # for the position seed (ngmix#796). All epochs of one object share the + # ra/dec above, so first-epoch CCD pins one deterministic seed stream. + if stamp.ccd is None: + stamp.ccd = int(ccd_n) return stamp @@ -1040,6 +1231,14 @@ def prepare_ngmix_weights(gal, weight, flag, rng, bkg_rms=None): else sigma_mad(gal) ) + # Guard: sig_noise=0 means galaxy is at the CCD edge (mostly-zero stamp). + # Division by zero would make weight_map NaN/inf, crashing GalSim C code. + # np.all keeps the guard valid for the per-pixel bkg_rms path (any zero + # pixel would produce the same NaN weight the scalar case guards against). + # MKDEBUG TODO this needs to be checked! + if not np.all(sig_noise > 0): + return np.zeros_like(gal), np.zeros_like(weight_map), np.zeros_like(gal) + noise_img = rng.standard_normal(gal.shape) * sig_noise noise_img_gal = rng.standard_normal(gal.shape) * sig_noise diff --git a/src/shapepipe/modules/ngmix_runner.py b/src/shapepipe/modules/ngmix_runner.py index 31028c7f8..2a18d60c3 100644 --- a/src/shapepipe/modules/ngmix_runner.py +++ b/src/shapepipe/modules/ngmix_runner.py @@ -8,6 +8,8 @@ import os +from sqlitedict import SqliteDict + from shapepipe.modules.module_decorator import module_runner from shapepipe.modules.ngmix_package.ngmix import Ngmix @@ -49,9 +51,18 @@ def ngmix_runner( # Pixel scale pixel_scale = config.getfloat(module_config_sec, "PIXEL_SCALE") + # Background subtraction: disable for image sims, where there is no + # background image to subtract. The background vignet occupies one input + # slot, so the f_wcs headers -- and the optional background-rms vignet -- + # shift by one when it is absent. + bkg_sub = config.getboolean(module_config_sec, "BKG_SUB", fallback=True) + wcs_idx = 6 if bkg_sub else 5 + # Path to merged single-exposure single-HDU headers - f_wcs_path = input_file_list[6] + f_wcs_path = input_file_list[wcs_idx] + # Optional background-rms vignet -> ngmix inverse-variance weight map (#779); + # all-or-nothing. Appended as the final positional input after the wcs slice. if config.has_option(module_config_sec, "BKG_RMS_VIGNET_PATH"): bkg_rms_vignet_path = config.getexpanded( module_config_sec, @@ -61,9 +72,9 @@ def ngmix_runner( raise FileNotFoundError( f"Background RMS vignet file not found: {bkg_rms_vignet_path}" ) - input_file_list = input_file_list[:6] + [bkg_rms_vignet_path] + input_file_list = input_file_list[:wcs_idx] + [bkg_rms_vignet_path] else: - input_file_list = input_file_list[:6] + input_file_list = input_file_list[:wcs_idx] # Batch save option if config.has_option(module_config_sec, "SAVE_BATCH"): @@ -79,15 +90,54 @@ def ngmix_runner( id_obj_min = config.getint(module_config_sec, "ID_OBJ_MIN") id_obj_max = config.getint(module_config_sec, "ID_OBJ_MAX") - # Centroid source for the galaxy Jacobian origin: "wcs" (default — the + # Centroid source for the galaxy Jacobian origin: "wcs" (default -- the # catalog sky position projected through the WCS, trusting the astrometry) # or "hsm" (legacy HSM adaptive-moment centroid, being phased out: noisy - # for stars and flagged as incorrect by Fabian — see #767). + # for stars and flagged as incorrect by Fabian -- see #767). if config.has_option(module_config_sec, "CENTROID_SOURCE"): centroid_source = config.get(module_config_sec, "CENTROID_SOURCE") else: centroid_source = "wcs" + # Seed the per-object RNG from sky position instead of per tile, so + # metacal's fixnoise counter-noise (and the fit guesses) cancel across + # Pujol image-simulation shear branches (ngmix#796). Default False leaves + # the production path byte-identical. + if config.has_option(module_config_sec, "SEED_FROM_POSITION"): + seed_from_position = config.getboolean( + module_config_sec, "SEED_FROM_POSITION" + ) + else: + seed_from_position = False + + # Check PSF vignets first: if all are empty dicts {}, the exposures for this + # tile are absent from the PSF dictionary and no shape measurement is possible. + # This check must come before reading image vignets to avoid a C-level malloc + # crash that occurs when large numbers of numpy arrays are allocated then freed. + psf_idx = 3 if bkg_sub else 2 + psf_vignet_path = input_file_list[psf_idx] + with SqliteDict(psf_vignet_path) as db: + psf_keys = list(db.keys()) + n_empty_psf = sum(1 for k in psf_keys if len(db[k]) == 0) + if psf_keys and n_empty_psf == len(psf_keys): + w_log.warning( + f"All {len(psf_keys)} PSF vignet entries are empty in " + f"{psf_vignet_path} — no PSF coverage for this tile. Skipping ngmix." + ) + return None, None + + # Check that image vignets are not all empty before initialising ngmix + image_vignet_path = input_file_list[1] + with SqliteDict(image_vignet_path) as db: + keys = list(db.keys()) + n_empty = sum(1 for k in keys if db[k] == "empty") + if keys and n_empty == len(keys): + w_log.warning( + f"All {len(keys)} image vignets are 'empty' in {image_vignet_path} " + "— no valid CCD coverage for this tile. Skipping ngmix." + ) + return None, None + # Initialise class instance ngmix_inst = Ngmix( input_file_list, @@ -100,7 +150,9 @@ def ngmix_runner( save_batch=save_batch, id_obj_min=id_obj_min, id_obj_max=id_obj_max, + bkg_sub=bkg_sub, centroid_source=centroid_source, + seed_from_position=seed_from_position, ) # Process ngmix shape measurement and metacalibration diff --git a/src/shapepipe/modules/sextractor_package/sextractor_script.py b/src/shapepipe/modules/sextractor_package/sextractor_script.py index f124612b7..89e3be6df 100644 --- a/src/shapepipe/modules/sextractor_package/sextractor_script.py +++ b/src/shapepipe/modules/sextractor_package/sextractor_script.py @@ -98,9 +98,13 @@ def make_post_process(cat_path, f_wcs_path, pos_params, ccd_size): history.append(idx) exp_list = [] - pattern = r"([0-9]*)p\.(.*)" + pattern = r"([0-9]+)p?\.(.*)" for hist in history: m = re.search(pattern, hist) + if m is None: + raise ValueError( + f"Could not parse exposure ID from HISTORY entry: '{hist}'" + ) exp_list.append(m.group(1)) obj_id = np.copy(cat.get_data()["NUMBER"]) diff --git a/src/shapepipe/modules/vignetmaker_package/vignetmaker.py b/src/shapepipe/modules/vignetmaker_package/vignetmaker.py index 74e26915b..505c037dd 100644 --- a/src/shapepipe/modules/vignetmaker_package/vignetmaker.py +++ b/src/shapepipe/modules/vignetmaker_package/vignetmaker.py @@ -323,7 +323,8 @@ def _get_stamp_me(self, image_dirs, image_pattern): else: self._w_log.info( f"All CCDs skipped for exposure {exp_name} " - + "(ccd == -1 for all); epoch excluded from output" + + f"[{image_pattern}]" + + " (ccd == -1 for all); epoch excluded from output" ) cat.close() @@ -380,6 +381,7 @@ def process_me(self, image_dirs, image_pattern, f_wcs_path, rad): else: index = idx + self._w_log.info(f"Processing {image_pattern[idx]} vignets") output_dict = self._get_stamp_me( image_dirs[index], image_pattern[idx], diff --git a/src/shapepipe/pipeline/exp_utils.py b/src/shapepipe/pipeline/exp_utils.py index 8f0df784a..b934f0bb3 100644 --- a/src/shapepipe/pipeline/exp_utils.py +++ b/src/shapepipe/pipeline/exp_utils.py @@ -22,6 +22,7 @@ def get_exp_output_files( file_pattern, file_ext, w_log=None, + warn_only=False, ): """Collect output files from a per-exposure runner for all tile exposures. @@ -31,9 +32,9 @@ def get_exp_output_files( ///output/run_sp_*//output/ - where ``exp_prefix = exp_id[:2]`` and ``exp_base = exp_id[:-1]`` (the - trailing letter, typically ``p``, is stripped because directory names - do not carry it). + where ``exp_prefix = exp_id[:2]`` and ``exp_base`` is ``exp_id`` with the + trailing letter stripped if present (e.g. ``2113864p`` → ``2113864``), or + the full ``exp_id`` for numeric-only IDs (image simulations). Parameters ---------- @@ -52,6 +53,10 @@ def get_exp_output_files( File extension including the leading dot, e.g. ``.npy`` w_log : logging.Logger, optional Pipeline logger; ``None`` silences all logging + warn_only : bool, optional + If ``True``, log a warning for missing exposures and continue with + the files that were found instead of raising ``FileNotFoundError``. + Default is ``False``. Returns ------- @@ -88,9 +93,10 @@ def get_exp_output_files( for exp_id in exp_ids: # Directory structure mirrors run_job_canfar_v2.0.sh: # exp_prefix = first 2 chars of exp_id (e.g. "21") - # exp_base = exp_id without trailing letter (e.g. "2113864") + # exp_base = exp_id without trailing letter if present (e.g. "2113864"), + # or full exp_id for numeric-only ids (image sims) exp_prefix = exp_id[:2] - exp_base = exp_id[:-1] + exp_base = exp_id[:-1] if exp_id[-1].isalpha() else exp_id pattern = os.path.join( exp_base_dir, @@ -115,10 +121,12 @@ def get_exp_output_files( w_log.warning(f" {exp_id}: no match for {pattern}") if missing: - raise FileNotFoundError( - f"No {runner_name} output found for " - f"{len(missing)} exposure(s): {missing}" - ) + msg = f"No {runner_name} output found for {len(missing)} exposure(s): {missing}" + if warn_only: + if w_log: + w_log.warning(msg) + else: + raise FileNotFoundError(msg) if w_log: w_log.info(f"Found {len(file_list)} exposure output files") @@ -183,7 +191,7 @@ def get_exp_output_dirs( for exp_id in exp_ids: exp_prefix = exp_id[:2] - exp_base = exp_id[:-1] + exp_base = exp_id[:-1] if exp_id[-1].isalpha() else exp_id pattern = os.path.join( exp_base_dir, diff --git a/src/shapepipe/testing/simulate.py b/src/shapepipe/testing/simulate.py index 8ab5a06f0..ee62f43c6 100644 --- a/src/shapepipe/testing/simulate.py +++ b/src/shapepipe/testing/simulate.py @@ -103,3 +103,98 @@ def make_data( jacob_lists.append(wcs.jacobian()) return gals, psfs, psfs_sigmas, weights, flags, jacob_lists + + +def make_data_pujol( + rng, + shear_list, + noise=1e-5, + n_epochs=1, + share_shift=False, + gal_hlr=0.3, + gal_flux=1000.0, + psf_fwhm=0.55, + pixel_scale=0.1857, + img_size=201, +): + """Simulate galaxies at multiple shear values sharing the same noise per epoch. + + Implements the Pujol (2018) denoising approach: all shear variants receive + the identical noise and sub-pixel shift each epoch, so noise contributions + to the shear response R and to the m/c estimator cancel exactly. + + Parameters + ---------- + rng : numpy.random.RandomState + shear_list : list of tuple of float + List of (g1, g2) values to simulate simultaneously. + noise : float, optional + Per-pixel noise sigma. + n_epochs : int, optional + share_shift : bool, optional + If True all epochs share the same sub-pixel shift; the shift is + always shared across shear variants regardless of this flag. + gal_hlr, gal_flux, psf_fwhm, pixel_scale, img_size : float / int, optional + Same meaning as in make_data. + + Returns + ------- + list of tuple + One ``(gals, psfs, psfs_sigmas, weights, flags, jacob_lists)`` per + entry in ``shear_list``, all sharing identical per-epoch noise. + """ + psf_noise_sigma = 1.0e-6 + scale = pixel_scale + wcs = galsim.PixelScale(scale) + n_shears = len(shear_list) + + # output accumulators + all_gals = [[] for _ in range(n_shears)] + all_psfs = [[] for _ in range(n_shears)] + all_psfs_sigmas = [[] for _ in range(n_shears)] + all_weights = [[] for _ in range(n_shears)] + all_flags = [[] for _ in range(n_shears)] + all_jacobs = [[] for _ in range(n_shears)] + + if share_shift: + dy, dx = rng.uniform(low=-scale / 2, high=scale / 2, size=2) + + for epoch in range(n_epochs): + if not share_shift: + dy, dx = rng.uniform(low=-scale / 2, high=scale / 2, size=2) + + # One noise draw shared by all shear variants this epoch + noise_img = rng.normal(scale=noise, size=(img_size, img_size)) + psf_noise_img = rng.normal(scale=psf_noise_sigma, size=(img_size, img_size)) + + psf = galsim.Moffat(beta=2.5, fwhm=psf_fwhm) + psf_im_ = psf.drawImage(nx=img_size, ny=img_size, wcs=wcs) + psf_sigma = galsim.hsm.FindAdaptiveMom(psf_im_).moments_sigma + psf_im = psf_im_.array.astype(np.float64) + psf_noise_img + + jacob = wcs.jacobian() + weight = np.full((img_size, img_size), 1.0 / noise ** 2) + flag = np.zeros((img_size, img_size)) + + for j, shear in enumerate(shear_list): + obj = galsim.Convolve( + psf, + galsim.Exponential(half_light_radius=gal_hlr, flux=gal_flux).shear( + g1=shear[0], g2=shear[1] + ), + ).shift(dx, dy) + im = obj.drawImage(nx=img_size, ny=img_size, wcs=wcs).array.astype(np.float64) + im += noise_img # same noise for every shear variant + + all_gals[j].append(im) + all_psfs[j].append(psf_im) + all_psfs_sigmas[j].append(psf_sigma) + all_weights[j].append(weight.copy()) + all_flags[j].append(flag.copy()) + all_jacobs[j].append(jacob) + + return [ + (all_gals[j], all_psfs[j], all_psfs_sigmas[j], + all_weights[j], all_flags[j], all_jacobs[j]) + for j in range(n_shears) + ] diff --git a/tests/module/test_ngmix.py b/tests/module/test_ngmix.py index 37db4e3fc..23b700aad 100644 --- a/tests/module/test_ngmix.py +++ b/tests/module/test_ngmix.py @@ -1,5 +1,6 @@ """UNIT TESTS FOR MODULE PACKAGE: NGMIX.""" +from astropy.io import fits from hypothesis import given from hypothesis import strategies as st import numpy as np @@ -135,6 +136,84 @@ def test_metacal_is_reproducible_with_fixed_seed(): npt.assert_array_equal(_metacal_noshear_g(42), _metacal_noshear_g(42)) +# --- position_seed: the Pujol noise-cancellation seed (ngmix#796) --- +# The seed is a pure function of (ra, dec, ccd) that must be (a) invariant to +# sub-box centroid jitter, (b) distinct for neighbouring boxes, (c) in the +# valid RandomState range, and (d) free of the anti-diagonal collision of +# Fabian's original ``box_x + box_y`` sum. These properties are what make the +# per-object noise/guess draws cancel across image-simulation shear branches. + +# In-survey positions: RA in [0, 360), Dec in [-40, 40], MegaCam CCDs 0..39. +_ra = st.floats(min_value=0.0, max_value=359.999, allow_nan=False) +_dec = st.floats(min_value=-40.0, max_value=40.0, allow_nan=False) +_ccd = st.integers(min_value=0, max_value=39) + + +@given(_ra, _dec, _ccd) +def test_position_seed_in_valid_randomstate_range(ra, dec, ccd): + """Seed always lands in [0, 2**32), the RandomState-acceptable range.""" + from shapepipe.modules.ngmix_package.ngmix import position_seed + + seed = position_seed(ra, dec, ccd) + assert 0 <= seed < 2 ** 32 + + +@given(_ra, _dec, _ccd, st.floats(min_value=-1.4, max_value=1.4)) +def test_position_seed_invariant_to_sub_box_jitter(ra, dec, ccd, jitter_arcsec): + """Centroid jitter < half the 3-arcsec box leaves the seed unchanged. + + An object re-detected with a slightly different centroid across branches + stays in the same box and so must get the same seed — the robustness that + lets the added noise cancel in the branch difference. We nudge the position + to a box centre first so a bounded jitter cannot cross a boundary. + """ + from shapepipe.modules.ngmix_package.ngmix import position_seed + + # Snap ra/dec to their box centres (1.5 arcsec into the 3-arcsec box), so a + # |jitter| < 1.5 arcsec cannot leave the box. + ra_c = (np.floor(ra * 3600 / 3) * 3 + 1.5) / 3600 + dec_c = (np.floor(dec * 3600 / 3) * 3 + 1.5) / 3600 + base = position_seed(ra_c, dec_c, ccd) + jittered = position_seed( + ra_c + jitter_arcsec / 3600, dec_c + jitter_arcsec / 3600, ccd + ) + assert base == jittered + + +@given(_ra, _dec, _ccd) +def test_position_seed_distinguishes_adjacent_boxes(ra, dec, ccd): + """Neighbouring 3-arcsec boxes get different seeds (in RA and in Dec). + + Both positions are snapped to box centres and the neighbour is one full box + (3 arcsec) away, so the shift genuinely lands in the next box rather than + grazing a boundary where ``3 / 3600``'s float error could leave ``floor`` + in the same box. + """ + from shapepipe.modules.ngmix_package.ngmix import position_seed + + ra_c = (np.floor(ra * 3600 / 3) * 3 + 1.5) / 3600 + dec_c = (np.floor(dec * 3600 / 3) * 3 + 1.5) / 3600 + base = position_seed(ra_c, dec_c, ccd) + assert position_seed(ra_c + 3 / 3600, dec_c, ccd) != base + assert position_seed(ra_c, dec_c + 3 / 3600, ccd) != base + + +def test_position_seed_resolves_fabian_anti_diagonal_collision(): + """The Cantor pairing avoids the ``box_x + box_y`` anti-diagonal collision. + + Fabian's #796 formula sums the box indices, so boxes ``(10, 20)`` and + ``(11, 19)`` collide (both -> 30) and two distinct objects would share a + noise stream. The box centres below realise exactly that pair; our pairing + must separate them. + """ + from shapepipe.modules.ngmix_package.ngmix import position_seed + + # Box centres for boxes (box_x, box_y) at ccd=0: box_x = floor(ra*1200)+1. + ra_a, dec_a = (10 - 1 + 0.5) / 1200, (20 - 2 + 0.5) / 1200 + ra_b, dec_b = (11 - 1 + 0.5) / 1200, (19 - 2 + 0.5) / 1200 + assert position_seed(ra_a, dec_a, 0) != position_seed(ra_b, dec_b, 0) + + # The two PSF families carry distinct ellipticity AND size (shapepipe#749): # the original image PSF (-> *_psf_orig) and the metacal reconvolution kernel # (-> *_psf_reconv) are independent fits, so a regression to the old aliasing @@ -310,6 +389,46 @@ def test_compile_results_nan_fills_failed_fit_types(): assert out["noshear"]["mcal_flags"] == [0x8] +def test_save_results_batches_survive_a_later_failed_fit(tmp_path): + """A batch-1-all-success / batch-2-has-a-failure sequence must not crash. + + With ``SAVE_BATCH`` enabled, ``save_results`` writes batch 1 fresh via + ``save_as_fits`` (which locks column dtypes from that data) and appends + every later batch into a structured array typed to those locked dtypes + (shapepipe#795). ``nfev_fit`` was the only int-like column with a NaN + fallback: if batch 1 has zero failed fits, ``nfev_fit`` locks to int64, + and the first later batch containing a failed fit fed a NaN into that + int64 column, raising ``ValueError: cannot convert float NaN to integer`` + at append time. The fix uses the -1 sentinel (ngmix's own + failed/absent-nfev convention) instead of NaN, so the column stays + int64 across every batch. + """ + from shapepipe.modules.ngmix_package.ngmix import Ngmix + + inst = object.__new__(Ngmix) + inst._zero_point = 30.0 + inst._output_dir = str(tmp_path) + inst._file_number_string = "-000-000" + + # Batch 1: every fit type succeeds -> nfev_fit locks to int64 on write. + batch1 = [_fake_metacal_result(0.18, 0.02, 0.09, 0.001)] + inst.save_results(inst.compile_results(batch1)) + + # Batch 2: one fit type fails (no "nfev" key -> hits the fallback path). + batch2_res = _fake_metacal_result(0.18, 0.02, 0.09, 0.001) + batch2_res["1p"] = {"flags": 0x8} # failed fit: no nfev at all + batch2_res["mcal_types_fail"] = 1 + batch2_res["mcal_flags"] = 0x8 + + # Must not raise (pre-fix: ValueError: cannot convert float NaN to integer). + inst.save_results(inst.compile_results([batch2_res])) + + with fits.open(inst.get_output_path(str(tmp_path))) as hdul: + nfev = hdul["1P"].data["nfev_fit"] + assert nfev.dtype.kind == "i" + npt.assert_array_equal(nfev, [1, -1]) + + def test_get_mcal_flags_ors_per_type_fit_flags(): """mcal_flags is the bitwise OR of all per-type fit flags (v1 contract).