From 5476e6e1598e9e351645ee589156bc8b956c9d77 Mon Sep 17 00:00:00 2001 From: Dali Wang Date: Thu, 19 Mar 2026 23:17:33 -0400 Subject: [PATCH 1/5] Implement Africa-focused 5P bias correction and fine-tuning configurations - Updated `commands.txt` to include new commands for applying bias/scale corrections specifically for Africa and Amazon regions. - Added new configuration files: `finetune_phase2_africa_5p_only.json`, `training_config_africa_5p_box.json`, `training_config_africa_five_p_only.json`, and `training_config_amazon_5p_box.json` to support targeted training and fine-tuning of models focused on phosphorus variables in Africa and Amazon. - Introduced documentation detailing the Africa 5P bias issue and the proposed shrinkage plan in `AFRICA_5P_SHRINKAGE_PLAN.md`. - Enhanced existing documentation to clarify the workflow for training and applying the new configurations, ensuring users can effectively implement the changes. - Added scripts for merging bias-corrected predictions from both regions and updating restart files accordingly. --- commands.txt | 48 +++- config/finetune_phase2_africa_5p_only.json | 136 ++++++++++ config/training_config_africa_5p_box.json | 143 ++++++++++ .../training_config_africa_five_p_only.json | 169 ++++++++++++ config/training_config_amazon_5p_box.json | 151 +++++++++++ ...ing_config_phase2_tropical_soilp_only.json | 24 +- docs/AFRICA_5P_SHRINKAGE_PLAN.md | 65 +++++ docs/FINETUNE_PHASE2_AFRICA_5P_SUMMARY.md | 38 +++ ...TIONS_TESNORTH10PC_INFERENCE_TO_RESTART.md | 173 ++++++++++++ ...STRUCTIONS_TRENDY_1_AI_RESTART_CREATION.md | 154 +++++------ .../WORKFLOW_PHASE1_PHASE2_PHASE3_RESTARTS.md | 240 +++++++++++++---- scripts/ai_predictions_to_netcdf.py | 7 +- scripts/apply_5p_africa_shrinkage.py | 252 ++++++++++++++++++ scripts/apply_5p_bias_scale_correction.py | 69 ++++- .../generate_site_5p_restart_comparison.py | 46 ++-- .../merge_5p_bias_corrected_amazon_africa.py | 176 ++++++++++++ scripts/run_phase1_global.sh | 3 + scripts/run_phase3_tworegions.sh | 36 ++- 18 files changed, 1740 insertions(+), 190 deletions(-) create mode 100644 config/finetune_phase2_africa_5p_only.json create mode 100644 config/training_config_africa_5p_box.json create mode 100644 config/training_config_africa_five_p_only.json create mode 100644 config/training_config_amazon_5p_box.json create mode 100644 docs/AFRICA_5P_SHRINKAGE_PLAN.md create mode 100644 docs/FINETUNE_PHASE2_AFRICA_5P_SUMMARY.md create mode 100644 docs/INSTRUCTIONS_TESNORTH10PC_INFERENCE_TO_RESTART.md create mode 100644 scripts/apply_5p_africa_shrinkage.py create mode 100644 scripts/merge_5p_bias_corrected_amazon_africa.py diff --git a/commands.txt b/commands.txt index b33e87c..3a90c8c 100644 --- a/commands.txt +++ b/commands.txt @@ -509,7 +509,7 @@ python scripts/extract_elm_restart_point.py \ --restart-file cnp_results/run_20260311_204845_phase3_tworegions/updated_restart_phase3_tworegions_5P_bias_corrected_tropical.nc \ --lat 0.0 \ --lon 28.0 \ - --output-file cnp_results/run_20260311_204845_phase3_tworegions/Africa_phase3_tworegions_restart.nc + --output-filcnp_results/run_20260313_224805_phase2_tropical_soilp_amazon_africae cnp_results/run_20260311_204845_phase3_tworegions/Africa_phase3_tworegions_restart.nc python scripts/generate_site_5p_restart_comparison.py --lon 303.75 --lat -17.434553 --site-name amazon_303_17S python scripts/generate_site_5p_restart_comparison.py --lon 292.5 --lat -15.55 --site-name amazon_292_15S @@ -519,4 +519,50 @@ python scripts/generate_site_5p_restart_comparison.py --lon 28 --lat 0 --s python scripts/generate_site_5p_restart_comparison.py --lon 14 --lat 0 --site-name africa +cd /mnt/proj-shared/AI4BGC_7xw/AI4BGC + +# Amazon-only bias/scale +python scripts/apply_5p_bias_scale_correction.py \ + --run-dir cnp_results/run_20260311_204845_phase3_tworegions\ + --region-config-json config/training_config_amazon_5p_box.json \ + --output-subdir soil_2d_predictions_5P_bias_corrected_amazon + +# Africa-only bias/scale +python scripts/apply_5p_bias_scale_correction.py \ + --run-dir cnp_results/run_20260311_204845_phase3_tworegions\ + --region-config-json config/training_config_africa_5p_box.json \ + --output-subdir soil_2d_predictions_5P_bias_corrected_africa + +python scripts/merge_5p_bias_corrected_amazon_africa.py \ + --run-dir cnp_results/run_20260311_204845_phase3_tworegions + + +python scripts/ai_predictions_to_netcdf.py \ + --ai-predictions cnp_results/run_20260311_204845_phase3_tworegions/cnp_inference_entire_dataset/cnp_predictions/ \ + --soil-2d-bias-corrected-subdir soil_2d_predictions_5P_bias_corrected_amazon_africa \ + --variable-list CNP_IO_updated9_dev_dw.txt \ + --output cnp_results/run_20260311_204845_phase3_tworegions/comparison_results/ai_predictions_5P_bias_corrected_amazon_africa.nc + + + python scripts/ai_predictions_to_restart.py \ + --ai-predictions cnp_results/run_20260311_204845_phase3_tworegions/comparison_results/ai_predictions_5P_bias_corrected_amazon_africa.nc \ + --restart-file cnp_results/run_20260311_174843_phase2_tropical/updated_restart_phase2_tropical_5P_raw.nc \ + --output cnp_results/run_20260311_204845_phase3_tworegions/updated_restart_phase3_tworegions_5P_bias_corrected_amazon_africa.nc \ + --variable-list CNP_IO_updated9_dev_dw.txt \ + --variables-to-update labilep_vr,occlp_vr,solutionp_vr,secondp_vr,primp_vr \ + "--tropical-lat-range=-30,30" + + python scripts/generate_site_5p_restart_comparison.py \ + --lon 28.0 \ + --lat 0.0 \ + --site-name africa_28_0_amazon_africa_nc \ + --output-dir cnp_results/run_20260311_204845_phase3_tworegions/analysis/africa_28_0_amazon_africa_5p_restart_comparison + + + python scripts/generate_site_5p_restart_comparison.py \ + --lon 303.75 \ + --lat -17.43 \ + --site-name amazon303_tworegion_biascorrection \ + --output-dir cnp_results/run_20260311_204845_phase3_tworegions/analysis/amazon_5p_restart_comparison_twocorrections + python train_cnp_model.py --training-config-json config/training_config_phase2_tropical_soilp_only.json --variable-list CNP_IO_updated9_dev_dw.txt --tropical-lat-range -30,30 --output-dir-suffix phase2_tropical_soilp_amazon_africa \ No newline at end of file diff --git a/config/finetune_phase2_africa_5p_only.json b/config/finetune_phase2_africa_5p_only.json new file mode 100644 index 0000000..abe1f62 --- /dev/null +++ b/config/finetune_phase2_africa_5p_only.json @@ -0,0 +1,136 @@ +{ + "_description": "Fine-tune Phase2 P-focused model on Africa region only, 5P-focused. Use: python scripts/run_finetuning_json.py --config-json config/finetune_phase2_africa_5p_only.json --output-dir-suffix phase2_africa_5p. Run from repo root. Improves Africa 5P without changing Amazon or other variables.", + "finetuning_paths": { + "variable_list_path": "../CNP_IO_updated9_dev_dw.txt", + "pretrained_model_path": "../cnp_results/run_20260305_153217_phase2_pvariable_focus/cnp_predictions/model.pth", + "dataset_path": "/mnt/proj-shared/AI4BGC_7xw/TrainingData/Trendy_1_data_CNP", + "file_pattern": "training_data_batch_*.pkl", + "output_dir": "../cnp_results", + "finetuned_model_filename": "finetuned_phase2_africa_5p.pth" + }, + "training_hyperparameters": { + "num_epochs": 50, + "batch_size": 128, + "learning_rate": 1e-5, + "optimizer_type": "adamw", + "weight_decay": 0.01, + "use_scheduler": true, + "scheduler_type": "cosine", + "scheduler_step_size": 10, + "scheduler_gamma": 0.1, + "xsmrpool_loss_weight": 10, + "litter_c_loss_weight": 1, + "litter_n_loss_weight": 1, + "litter_p_loss_weight": 4, + "scalar_loss_weight": 0.06, + "vector_loss_weight": 0.4, + "matrix_loss_weight": 2.5 + }, + "reproducibility_config": { + "random_seed": 42, + "strict_determinism": false, + "train_split": 0.8, + "normalization": "individual", + "dropout_p": 0.12 + }, + "data_filtering_config": { + "tropical_only": false, + "tropical_lat_range": [-30.0, 30.0], + "longitudes_to_drop": [0, 358.75], + "natveg_only": true, + "natveg_filter_before_split": false, + "region_boxes": [ + [-15, 15, 0, 30] + ] + }, + "variable_weights": { + "pft1d_weights": { + "cpool": 3, + "npool": 4, + "ppool": 6, + "totvegc": 3, + "leafc": 3, + "leafn": 3, + "leafp": 4, + "leafc_storage": 2, + "leafn_storage": 2, + "leafp_storage": 3, + "deadstemc": 3, + "deadstemn": 3, + "deadstemp": 4, + "frootc": 2, + "frootn": 2, + "frootp": 3, + "frootc_storage": 1.5, + "frootn_storage": 1.5, + "frootp_storage": 2, + "tlai": 2 + }, + "soil2d_weights": { + "cwdc_vr": 0.1, + "cwdn_vr": 0.1, + "cwdp_vr": 0.1, + "primp_vr": 35, + "occlp_vr": 50, + "labilep_vr": 40, + "secondp_vr": 40, + "solutionp_vr": 40, + "litr2c_vr": 0.1, + "litr2p_vr": 0.1, + "litr2n_vr": 0.1, + "litr3c_vr": 0.1, + "litr3n_vr": 0.1, + "litr3p_vr": 0.1, + "soil1c_vr": 0.1, + "soil1n_vr": 0.1, + "soil1p_vr": 0.1, + "soil2c_vr": 0.1, + "soil2n_vr": 0.1, + "soil2p_vr": 0.1, + "soil3c_vr": 0.1, + "soil3n_vr": 0.1, + "soil3p_vr": 0.1, + "soil4c_vr": 0.1, + "soil4n_vr": 0.1, + "soil4p_vr": 0.1 + }, + "scalar_weights": { + "GPP": 0.08, + "NPP": 0.08, + "AR": 0.08, + "HR": 0.08 + } + }, + "tail_aware_config": { + "loss": "log1p_huber", + "epsilon": 1e-8, + "huber_delta": 0.7 + }, + "tail_aware_weights": { + "cpool": 3, + "deadstemc": 3, + "deadcrootc": 3, + "livestemc": 3, + "livecrootc": 3, + "npool": 4, + "ppool": 6, + "totvegc": 3, + "leafc": 3, + "leafn": 3, + "leafp": 4, + "labilep": 10, + "labilep_vr": 40, + "primp_vr": 35, + "occlp_vr": 50, + "secondp_vr": 40, + "solutionp_vr": 40 + }, + "pft_mask_config": { + "mask_absent_pfts": true, + "pft_presence_threshold": 0 + }, + "pft_zero_sparsity_config": { + "weight": 1, + "threshold": 1e-8 + } +} diff --git a/config/training_config_africa_5p_box.json b/config/training_config_africa_5p_box.json new file mode 100644 index 0000000..bde78d4 --- /dev/null +++ b/config/training_config_africa_5p_box.json @@ -0,0 +1,143 @@ +{ + "_description": "Africa-only 5P region config: used to drive bias/scale correction on Africa box only. Data filtering keeps only Africa natveg cells for fitting, but model training can still be global. Primary use is with apply_5p_bias_scale_correction.py via --region-config-json. Africa box: lat [-15, 15], lon [0, 30].", + "training_hyperparameters": { + "num_epochs": 150, + "batch_size": 128, + "learning_rate": 5e-5, + "optimizer_type": "adamw", + "weight_decay": 0.01, + "use_scheduler": true, + "scheduler_type": "cosine", + "scheduler_step_size": 10, + "scheduler_gamma": 0.1, + "xsmrpool_loss_weight": 10, + "litter_c_loss_weight": 1, + "litter_n_loss_weight": 1, + "litter_p_loss_weight": 4, + "scalar_loss_weight": 0.06, + "vector_loss_weight": 0.4, + "matrix_loss_weight": 2.5 + }, + "reproducibility_config": { + "random_seed": 42, + "strict_determinism": false, + "train_split": 0.8, + "normalization": "individual", + "dropout_p": 0.12 + }, + "data_filtering_config": { + "tropical_only": false, + "tropical_lat_range": [-30.0, 30.0], + "longitudes_to_drop": [0, 358.75], + "natveg_only": true, + "natveg_filter_before_split": false, + "region_boxes": [ + [-15.0, 15.0, 0.0, 30.0] + ] + }, + "variable_weights": { + "pft1d_weights": { + "cpool": 3, + "npool": 4, + "ppool": 6, + "totvegc": 3, + "leafc": 3, + "leafn": 3, + "leafp": 4, + "leafc_storage": 2, + "leafn_storage": 2, + "leafp_storage": 3, + "deadstemc": 3, + "deadstemn": 3, + "deadstemp": 4, + "frootc": 2, + "frootn": 2, + "frootp": 3, + "frootc_storage": 1.5, + "frootn_storage": 1.5, + "frootp_storage": 2, + "tlai": 2 + }, + "soil2d_weights": { + "cwdc_vr": 0.1, + "cwdn_vr": 0.1, + "cwdp_vr": 0.1, + "primp_vr": 35, + "occlp_vr": 50, + "labilep_vr": 40, + "secondp_vr": 40, + "solutionp_vr": 40, + "litr2c_vr": 0.1, + "litr2p_vr": 0.1, + "litr2n_vr": 0.1, + "litr3c_vr": 0.1, + "litr3n_vr": 0.1, + "litr3p_vr": 0.1, + "soil1c_vr": 0.1, + "soil1n_vr": 0.1, + "soil1p_vr": 0.1, + "soil2c_vr": 0.1, + "soil2n_vr": 0.1, + "soil2p_vr": 0.1, + "soil3c_vr": 0.1, + "soil3n_vr": 0.1, + "soil3p_vr": 0.1, + "soil4c_vr": 0.1, + "soil4n_vr": 0.1, + "soil4p_vr": 0.1 + }, + "scalar_weights": { + "GPP": 0.08, + "NPP": 0.08, + "AR": 0.08, + "HR": 0.08 + } + }, + "tail_aware_config": { + "loss": "log1p_huber", + "epsilon": 1e-8, + "huber_delta": 0.7 + }, + "tail_aware_weights": { + "cpool": 3, + "deadstemc": 3, + "deadcrootc": 3, + "livestemc": 3, + "livecrootc": 3, + "npool": 4, + "ppool": 6, + "totvegc": 3, + "leafc": 3, + "leafn": 3, + "leafp": 4, + "labilep": 10, + "labilep_vr": 40, + "primp_vr": 35, + "occlp_vr": 50, + "secondp_vr": 40, + "solutionp_vr": 40 + }, + "pft_mask_config": { + "mask_absent_pfts": true, + "pft_presence_threshold": 0 + }, + "pft_zero_sparsity_config": { + "weight": 1, + "threshold": 1e-8 + }, + "pft_zero_sparsity_weights": { + "cpool": 2, + "deadstemc": 1, + "deadcrootc": 1, + "livestemc": 1, + "livecrootc": 1 + }, + "pft1d_activation_overrides": { + "cpool": "abs", + "deadstemc": "abs", + "deadcrootc": "abs", + "livestemc": "abs", + "livecrootc": "abs" + } +} + diff --git a/config/training_config_africa_five_p_only.json b/config/training_config_africa_five_p_only.json new file mode 100644 index 0000000..5f7389e --- /dev/null +++ b/config/training_config_africa_five_p_only.json @@ -0,0 +1,169 @@ +{ + "_description": "Africa-only 5P model: train on Africa natveg grid cells (region box [-15, 15, 0, 30]) with loss dominated by the five soil P variables. Intended to provide better 5P predictions in Africa without changing Amazon or other regions. Train: python train_cnp_model.py --training-config-json config/training_config_africa_five_p_only.json --variable-list CNP_IO_updated9_dev_dw.txt --use-tva4km --output-dir-suffix africa_five_p", + "training_hyperparameters": { + "num_epochs": 120, + "batch_size": 128, + "learning_rate": 5e-5, + "optimizer_type": "adamw", + "weight_decay": 0.01, + "use_scheduler": true, + "scheduler_type": "cosine", + "scheduler_step_size": 10, + "scheduler_gamma": 0.1, + "xsmrpool_loss_weight": 10, + "litter_c_loss_weight": 1, + "litter_n_loss_weight": 1, + "litter_p_loss_weight": 4, + "scalar_loss_weight": 0.06, + "vector_loss_weight": 0.4, + "matrix_loss_weight": 2.5 + }, + "reproducibility_config": { + "random_seed": 42, + "strict_determinism": false, + "train_split": 0.8, + "normalization": "individual", + "dropout_p": 0.12 + }, + "data_filtering_config": { + "tropical_only": false, + "tropical_lat_range": [-30.0, 30.0], + "longitudes_to_drop": [0, 358.75], + "natveg_only": true, + "natveg_filter_before_split": false, + "region_boxes": [ + [-15.0, 15.0, 0.0, 30.0] + ] + }, + "variable_weights": { + "pft1d_weights": { + "cpool": 3, + "npool": 4, + "ppool": 6, + "totvegc": 3, + "leafc": 3, + "leafn": 3, + "leafp": 4, + "leafc_storage": 2, + "leafn_storage": 2, + "leafp_storage": 3, + "deadstemc": 3, + "deadstemn": 3, + "deadstemp": 4, + "frootc": 2, + "frootn": 2, + "frootp": 3, + "frootc_storage": 1.5, + "frootn_storage": 1.5, + "frootp_storage": 2, + "tlai": 2 + }, + "soil2d_weights": { + "cwdc_vr": 0.1, + "cwdn_vr": 0.1, + "cwdp_vr": 0.1, + "primp_vr": 35, + "occlp_vr": 50, + "labilep_vr": 40, + "secondp_vr": 40, + "solutionp_vr": 40, + "litr2c_vr": 0.1, + "litr2p_vr": 0.1, + "litr2n_vr": 0.1, + "litr3c_vr": 0.1, + "litr3n_vr": 0.1, + "litr3p_vr": 0.1, + "soil1c_vr": 0.1, + "soil1n_vr": 0.1, + "soil1p_vr": 0.1, + "soil2c_vr": 0.1, + "soil2n_vr": 0.1, + "soil2p_vr": 0.1, + "soil3c_vr": 0.1, + "soil3n_vr": 0.1, + "soil3p_vr": 0.1, + "soil4c_vr": 0.1, + "soil4n_vr": 0.1, + "soil4p_vr": 0.1 + }, + "scalar_weights": { + "GPP": 0.08, + "NPP": 0.08, + "AR": 0.08, + "HR": 0.08 + } + }, + "tail_aware_config": { + "loss": "log1p_huber", + "epsilon": 1e-8, + "huber_delta": 0.7 + }, + "tail_aware_weights": { + "cpool": 3, + "deadstemc": 3, + "deadcrootc": 3, + "livestemc": 3, + "livecrootc": 3, + "npool": 4, + "ppool": 6, + "totvegc": 3, + "leafc": 3, + "leafn": 3, + "leafp": 4, + "leafc_storage": 2, + "leafn_storage": 2, + "leafp_storage": 3, + "deadstemp": 4, + "deadstemn": 3, + "frootc": 2, + "frootn": 2, + "frootp": 3, + "labilep": 10, + "labilep_vr": 40, + "primp_vr": 35, + "occlp_vr": 50, + "secondp_vr": 40, + "solutionp_vr": 40, + "litr2c_vr": 0.1, + "litr2p_vr": 0.1, + "litr2n_vr": 0.1, + "litr3c_vr": 0.1, + "litr3p_vr": 0.1, + "litr3n_vr": 0.1, + "soil1c_vr": 0.1, + "soil1n_vr": 0.1, + "soil1p_vr": 0.1, + "soil2c_vr": 0.1, + "soil2n_vr": 0.1, + "soil2p_vr": 0.1, + "soil3c_vr": 0.1, + "soil3n_vr": 0.1, + "soil3p_vr": 0.1, + "soil4c_vr": 0.1, + "soil4n_vr": 0.1, + "soil4p_vr": 0.1 + }, + "pft_mask_config": { + "mask_absent_pfts": true, + "pft_presence_threshold": 0 + }, + "pft_zero_sparsity_config": { + "weight": 1, + "threshold": 1e-8 + }, + "pft_zero_sparsity_weights": { + "cpool": 2, + "deadstemc": 1, + "deadcrootc": 1, + "livestemc": 1, + "livecrootc": 1 + }, + "pft1d_activation_overrides": { + "cpool": "abs", + "deadstemc": "abs", + "deadcrootc": "abs", + "livestemc": "abs", + "livecrootc": "abs" + } +} + diff --git a/config/training_config_amazon_5p_box.json b/config/training_config_amazon_5p_box.json new file mode 100644 index 0000000..64347b1 --- /dev/null +++ b/config/training_config_amazon_5p_box.json @@ -0,0 +1,151 @@ +{ + "_description": "Amazon-only 5P region config: used to drive bias/scale correction on Amazon box only. Data filtering keeps only Amazon natveg cells for fitting, but model training can still be global. Primary use is with apply_5p_bias_scale_correction.py via --region-config-json. Amazon box: lat [-30, 10], lon [270, 330].", + "training_hyperparameters": { + "num_epochs": 150, + "batch_size": 128, + "learning_rate": 5e-5, + "optimizer_type": "adamw", + "weight_decay": 0.01, + "use_scheduler": true, + "scheduler_type": "cosine", + "scheduler_step_size": 10, + "scheduler_gamma": 0.1, + "xsmrpool_loss_weight": 10, + "litter_c_loss_weight": 1, + "litter_n_loss_weight": 1, + "litter_p_loss_weight": 4, + "scalar_loss_weight": 0.06, + "vector_loss_weight": 0.4, + "matrix_loss_weight": 2.5 + }, + "reproducibility_config": { + "random_seed": 42, + "strict_determinism": false, + "train_split": 0.8, + "normalization": "individual", + "dropout_p": 0.12 + }, + "data_filtering_config": { + "tropical_only": false, + "tropical_lat_range": [-30.0, 30.0], + "longitudes_to_drop": [0, 358.75], + "natveg_only": true, + "natveg_filter_before_split": false, + "region_boxes": [ + [-30.0, 10.0, 270.0, 330.0] + ] + }, + "variable_weights": { + "pft1d_weights": { + "cpool": 3, + "npool": 4, + "ppool": 6, + "totvegc": 3, + "leafc": 3, + "leafn": 3, + "leafp": 4, + "leafc_storage": 2, + "leafn_storage": 2, + "leafp_storage": 3, + "deadstemc": 3, + "deadstemn": 3, + "deadstemp": 4, + "frootc": 2, + "frootn": 2, + "frootp": 3, + "frootc_storage": 1.5, + "frootn_storage": 1.5, + "frootp_storage": 2, + "tlai": 2 + }, + "soil2d_weights": { + "cwdc_vr": 0.1, + "cwdn_vr": 0.1, + "cwdp_vr": 0.1, + "primp_vr": 35, + "occlp_vr": 50, + "labilep_vr": 40, + "secondp_vr": 40, + "solutionp_vr": 40, + "litr2c_vr": 0.1, + "litr2p_vr": 0.1, + "litr2n_vr": 0.1, + "litr3c_vr": 0.1, + "litr3n_vr": 0.1, + "litr3p_vr": 0.1, + "soil1c_vr": 0.1, + "soil1n_vr": 0.1, + "soil1p_vr": 0.1, + "soil2c_vr": 0.1, + "soil2n_vr": 0.1, + "soil2p_vr": 0.1, + "soil3c_vr": 0.1, + "soil3n_vr": 0.1, + "soil3p_vr": 0.1, + "soil4c_vr": 0.1, + "soil4n_vr": 0.1, + "soil4p_vr": 0.1 + }, + "scalar_weights": { + "GPP": 0.08, + "NPP": 0.08, + "AR": 0.08, + "HR": 0.08 + } + }, + "tail_aware_config": { + "loss": "log1p_huber", + "epsilon": 1e-8, + "huber_delta": 0.7 + }, + "tail_aware_weights": { + "cpool": 3, + "deadstemc": 3, + "deadcrootc": 3, + "livestemc": 3, + "livecrootc": 3, + "npool": 4, + "ppool": 6, + "totvegc": 3, + "leafc": 3, + "leafn": 3, + "leafp": 4, + "leafc_storage": 2, + "leafn_storage": 2, + "leafp_storage": 3, + "deadstemp": 4, + "deadstemn": 3, + "frootc": 2, + "frootn": 2, + "frootp": 3, + "labilep": 10, + "labilep_vr": 40, + "primp_vr": 35, + "occlp_vr": 50, + "secondp_vr": 40, + "solutionp_vr": 40 + }, + "pft_mask_config": { + "mask_absent_pfts": true, + "pft_presence_threshold": 0 + }, + "pft_zero_sparsity_config": { + "weight": 1, + "threshold": 1e-8 + }, + "pft_zero_sparsity_weights": { + "cpool": 2, + "deadstemc": 1, + "deadcrootc": 1, + "livestemc": 1, + "livecrootc": 1 + }, + "pft1d_activation_overrides": { + "cpool": "abs", + "deadstemc": "abs", + "deadcrootc": "abs", + "livestemc": "abs", + "livecrootc": "abs" + } +} + diff --git a/config/training_config_phase2_tropical_soilp_only.json b/config/training_config_phase2_tropical_soilp_only.json index f0e0a21..9ebe4cb 100644 --- a/config/training_config_phase2_tropical_soilp_only.json +++ b/config/training_config_phase2_tropical_soilp_only.json @@ -1,5 +1,5 @@ { - "_description": "Phase 2 tropical-only model: Train on tropical natveg grid cells with EXTREME focus on soil P variables. Use with Phase 1 global restart: Phase 1 produces global restart; Phase 2 predictions replace variables in tropical cells only (via ai_predictions_to_restart --tropical-lat-range). Train: python train_cnp_model.py --training-config-json config/training_config_phase2_tropical_soilp_only.json --variable-list CNP_IO_updated9_dev_dw.txt --use-tva4km --tropical-only --tropical-lat-range -30,30 --output-dir-suffix phase2_tropical_soilp_only", + "_description": "Phase 2 tropical-only model (Amazon + Africa): Train on tropical natveg grid cells with strong focus on soil P variables; solutionp_vr weighted highest. Use with Phase 1 global restart. Train: python train_cnp_model.py --training-config-json config/training_config_phase2_tropical_soilp_only.json --variable-list CNP_IO_updated9_dev_dw.txt --tropical-lat-range -30,30 --output-dir-suffix phase2_tropical_soilp_only", "training_hyperparameters": { "num_epochs": 150, "batch_size": 128, @@ -15,7 +15,7 @@ "litter_n_loss_weight": 1, "litter_p_loss_weight": 4, "scalar_loss_weight": 0.06, - "vector_loss_weight": 0.4, + "vector_loss_weight": 0.5, "matrix_loss_weight": 2.5 }, "reproducibility_config": { @@ -59,11 +59,11 @@ "cwdc_vr": 1, "cwdn_vr": 1, "cwdp_vr": 2, - "primp_vr": 26, - "occlp_vr": 45, - "labilep_vr": 18, - "secondp_vr": 14, - "solutionp_vr": 18, + "primp_vr": 28, + "occlp_vr": 42, + "labilep_vr": 28, + "secondp_vr": 24, + "solutionp_vr": 50, "litr2c_vr": 4, "litr2p_vr": 16, "litr2n_vr": 4, @@ -116,11 +116,11 @@ "frootn": 2, "frootp": 3, "labilep": 10, - "labilep_vr": 18, - "primp_vr": 26, - "occlp_vr": 45, - "secondp_vr": 14, - "solutionp_vr": 18, + "labilep_vr": 28, + "primp_vr": 28, + "occlp_vr": 42, + "secondp_vr": 24, + "solutionp_vr": 50, "litr2c_vr": 4, "litr2p_vr": 16, "litr2n_vr": 4, diff --git a/docs/AFRICA_5P_SHRINKAGE_PLAN.md b/docs/AFRICA_5P_SHRINKAGE_PLAN.md new file mode 100644 index 0000000..c924e86 --- /dev/null +++ b/docs/AFRICA_5P_SHRINKAGE_PLAN.md @@ -0,0 +1,65 @@ +## Africa 5P bias issue and shrinkage plan + +### Diagnosis + +- **Problem**: In central Africa, all 5P predictions (`labilep_vr`, `occlp_vr`, `solutionp_vr`, `secondp_vr`, `primp_vr`) are systematically too high in both the global and Phase2 tropical models. +- **Magnitude**: + - `occlp_vr` ground truth max at key Africa sites is ~150, but predictions are ~400. + - `solutionp_vr` ground truth max is ~0.00125, but predictions after bias correction can be ~0.03 (10–30× larger). +- **Fine-tuning result**: Africa-only 5P fine-tuning of Phase2 reduces loss but does **not** fix the magnitude issue at the Africa validation site (27.5E, 2.36N); predictions remain orders of magnitude above GT. +- **Bias/scale limitation**: Region-wise linear corrections (a * pred + b) cannot: + - Pull down extreme Africa values into the GT range, **and** + - Preserve Phase2’s good behavior in Amazon and other regions. + +### Constraints + +- Preserve **Amazon** 5P performance (already excellent after bias/scale). +- Do **not** change non-P variables. +- Keep existing Phase2 / Phase3 restart pipeline as much as possible. + +### Proposed solution: Africa-only 5P shrinkage layer + +1. **Keep existing pipeline as-is outside Africa**: + - Global natveg_improved model → Phase1 restart. + - Phase2 P-focused model → full-grid/tropical predictions. + - Amazon + Africa bias/scale correction step for 5P (current `apply_5p_bias_scale_correction.py`). + - `ai_predictions_to_restart.py` to write the final restart. + +2. **Add an Africa-only shrinkage correction for 5P variables**: + - Restrict to the **Africa region box**: lat ∈ [-15, 15], lon ∈ [0, 30]. + - For each 5P variable, use Africa GT and Phase2 predictions to compute a **shrinkage factor per layer**: + - For cells with positive GT and prediction, compute `ratio = gt / pred` in the Africa box. + - For each layer, take a robust statistic such as `median(ratio)` over Africa. + - Clamp the factor for that layer to a safe range, e.g. `factor ∈ [0.0, 1.0]` with lower bound > 0. + - Apply `corrected = pred * factor` **only in Africa**, leaving Amazon and other regions untouched. + - Optionally clip corrected Africa values per-layer to an upper bound such as `p99(GT_Africa_layer)` to avoid any remaining outliers. + +3. **Implement as a separate script** `scripts/apply_5p_africa_shrinkage.py`: + - Inputs: + - `--run-dir`: a `cnp_results` run that already has + `cnp_inference_entire_dataset/cnp_predictions/soil_2d_ground_truth` and + `soil_2d_predictions_5P_bias_corrected_phase2/`. + - For each 5P variable: + - Read `ground_truth_Y_{var}.csv` and `predictions_Y_{var}_bias_corrected.csv`. + - Filter rows to the Africa box. + - Compute per-layer shrinkage factor as described above. + - Apply `pred * factor` for rows in Africa region, leave others unchanged. + - Write corrected CSVs to a new subdir, e.g. `soil_2d_predictions_5P_bias_corrected_africa_shrinkage/`. + - Save the shrinkage factors and basic diagnostics (per-layer factors, min/max GT and pred in Africa) to a JSON file under `run_dir/analysis/`. + +4. **Use shrinkage output in restart creation**: + - Modify the Phase3 workflow so that: + 1. Run `apply_5p_bias_scale_correction.py` (produces `soil_2d_predictions_5P_bias_corrected_phase2`). + 2. Run `apply_5p_africa_shrinkage.py` to produce Africa-shrunk 5P CSVs. + 3. Run `ai_predictions_to_netcdf.py` **pointing at the Africa-shrunk subdir** instead of the original bias-corrected subdir. + 4. Run `ai_predictions_to_restart.py` as before with `--tropical-lat-range=-30,30`. + +5. **Validation plan**: + - For selected Africa sites (including the problematic site near 28E, 0N and others), compare vertical profiles and MAE/RMSE for the 5P variables across: + - Phase2 raw, + - Phase2 + bias/scale, + - Phase2 + bias/scale + Africa shrinkage. + - Confirm: + - Africa 5P magnitudes are closer to GT (within a factor ~2 instead of 10–100×). + - Amazon plots are unchanged relative to the current Phase3 restart. + diff --git a/docs/FINETUNE_PHASE2_AFRICA_5P_SUMMARY.md b/docs/FINETUNE_PHASE2_AFRICA_5P_SUMMARY.md new file mode 100644 index 0000000..6061c42 --- /dev/null +++ b/docs/FINETUNE_PHASE2_AFRICA_5P_SUMMARY.md @@ -0,0 +1,38 @@ +# Summary and reason for fine-tuning Phase2 on Africa only (5P-focused) + +## Reason for fine-tuning + +- **Africa 5P predictions are far too high** from both the global and the tropical (Phase2) models. Examples: + - **occlp_vr:** ground truth max ~150 at Africa sites; predictions ~400. + - **solutionp_vr:** ground truth max ~0.00125; predictions after bias correction ~0.03 (order-of-magnitude too high). +- **Bias/scale correction alone cannot fix this:** the underlying model bias in Africa is too large; a linear correction either fails to bring values into a plausible range or overfits and degrades elsewhere. +- **We do not want to change other parts or variables:** Amazon and the rest of the tropics are working well with Phase2 + bias/scale. The goal is to improve **only** Africa 5P without altering Amazon or non-P variables. + +Therefore we **fine-tune the Phase2 model on Africa-only data**, with loss focused on the five P variables, so that the model learns Africa-specific 5P structure. The resulting checkpoint can then be used for Africa-region inference and its 5P predictions spliced into the restart in the Africa box only. + +## Summary + +| Item | Description | +|------|-------------| +| **Base model** | Phase2 P-focused: `cnp_results/run_20260305_153217_phase2_pvariable_focus/cnp_predictions/model.pth` | +| **Fine-tune data** | Africa region only: lat ∈ [-15, 15], lon ∈ [0, 30] (single region box). | +| **Loss focus** | Five P variables: labilep_vr, occlp_vr, solutionp_vr, secondp_vr, primp_vr (high weights; other variables down-weighted). | +| **Config** | `config/finetune_phase2_africa_5p_only.json` | +| **Output** | Timestamped run under `cnp_results/` with suffix `phase2_africa_5p`; checkpoint e.g. `finetuned_phase2_africa_5p.pth`. | +| **Use after training** | Run inference with the Africa fine-tuned model in the Africa box only; merge those 5P values into the restart (Phase2 + bias/scale elsewhere). | + +## Two-region model not used + +The two-region 5P model (`run_20260308_211537_two_region_five_p`) was evaluated and performs **worse** than Phase2 at the Amazon site. We do not use it for production. See `docs/AMAZON_5P_COMPARISON_TWO_REGION_VS_PHASE2_REPORT.md`. + +## How to run the fine-tuning + +From the repo root: + +```bash +cd /mnt/proj-shared/AI4BGC_7xw/AI4BGC + +python scripts/run_finetuning_json.py \ + --config-json config/finetune_phase2_africa_5p_only.json \ + --output-dir-suffix phase2_africa_5p +``` diff --git a/docs/INSTRUCTIONS_TESNORTH10PC_INFERENCE_TO_RESTART.md b/docs/INSTRUCTIONS_TESNORTH10PC_INFERENCE_TO_RESTART.md new file mode 100644 index 0000000..3646491 --- /dev/null +++ b/docs/INSTRUCTIONS_TESNORTH10PC_INFERENCE_TO_RESTART.md @@ -0,0 +1,173 @@ +# Instructions: TESNORTHERA5 restart update from a trained CNP + +This document walks through: +1) running inference with your trained model, +2) exporting predictions to a NetCDF (`.nc`), +3) updating the target ELM restart NetCDF using those predictions. + +All commands below are based on the repo scripts: +- `scripts/run_inference_all.py` +- `scripts/ai_predictions_to_netcdf.py` +- `scripts/ai_predictions_to_restart.py` + +--- + +## Inputs (your paths) + +### Trained model run (source of checkpoint) +- `TRAIN_RUN_DIR=/mnt/proj-shared/AI4BGC_7xw/AI4BGC/cnp_results/run_20260319_134216_tesnorth10pc` +- Model checkpoint: + - `$TRAIN_RUN_DIR/cnp_predictions/model.pth` + +### CNP IO variable list (drives variables + dataset defaults) +- `VAR_LIST=/mnt/proj-shared/AI4BGC_7xw/AI4BGC/CNP_IO_tesnorth10pct.txt` + +### Target restart NetCDF to update +- `RESTART_FILE=/mnt/proj-shared/AI4BGC_7xw/AI4BGC/ELM_data/TESNORTHERA5_20adspin_restartfile/uELM_NORTHERA5_ERA5REF_I1850uELMCNPRDCTCBC.elm.r.0021-01-01-00000.nc` + +--- + +## Step 0: Set up run output paths + +```bash +cd /mnt/proj-shared/AI4BGC_7xw/AI4BGC + +export TRAIN_RUN_DIR="/mnt/proj-shared/AI4BGC_7xw/AI4BGC/cnp_results/run_20260319_134216_tesnorth10pc" +export VAR_LIST="/mnt/proj-shared/AI4BGC_7xw/AI4BGC/CNP_IO_tesnorth10pct.txt" +export RESTART_FILE="/mnt/proj-shared/AI4BGC_7xw/AI4BGC/ELM_data/TESNORTHERA5_20adspin_restartfile/uELM_NORTHERA5_ERA5REF_I1850uELMCNPRDCTCBC.elm.r.0021-01-01-00000.nc" + +export OUT_ROOT="$TRAIN_RUN_DIR/inference_restart_update" +mkdir -p "$OUT_ROOT" + +# Where inference predictions (CSV) will land +export INFER_DIR="$OUT_ROOT/cnp_inference_entire_dataset" + +# Where the converter will write the NetCDF needed by the restart updater +export AI_PRED_NC="$OUT_ROOT/ai_predictions_for_restart.nc" + +# Where the restart updater will write the updated restart +export OUT_RESTART="$OUT_ROOT/updated_restart_tesnorth10pc.nc" +``` + +--- + +## Step 1: Run inference with your trained checkpoint + +### Recommended mode (use dataset defaults referenced by `CNP_IO_tesnorth10pct.txt` / training config) + +```bash +python scripts/run_inference_all.py \ + --model "$TRAIN_RUN_DIR/cnp_predictions/model.pth" \ + --variable-list "$VAR_LIST" \ + --output-dir "$INFER_DIR" +``` + +This runs inference on the entire dataset available to the script (it will use `data_paths` / `file_pattern` from the training run config if available; otherwise it falls back to defaults from the provided IO file). + +### If your "new dataset" is stored in different PKL directories + +If the inference data (`training_data_batch_*.pkl`) live somewhere else, pass: +- `--data-paths` (comma-separated directories) +- `--file-pattern` (e.g. `training_data_batch_*.pkl`) + +Example: +```bash +python scripts/run_inference_all.py \ + --model "$TRAIN_RUN_DIR/cnp_predictions/model.pth" \ + --variable-list "$VAR_LIST" \ + --data-paths "/path/to/new/pkl_root_1,/path/to/new/pkl_root_2" \ + --file-pattern "training_data_batch_*.pkl" \ + --output-dir "$INFER_DIR" +``` + +--- + +## Step 2: Convert predictions to NetCDF + +The restart updater expects the NetCDF layout produced by `ai_predictions_to_netcdf.py`. + +```bash +python scripts/ai_predictions_to_netcdf.py \ + --ai-predictions "$INFER_DIR/cnp_predictions" \ + --variable-list "$VAR_LIST" \ + --output "$AI_PRED_NC" +``` + +### Longitude wrapping (only if mapping fails) + +`ai_predictions_to_restart.py` maps AI gridcells to model gridcells by nearest-neighbor in `(lon,lat)`. +So the longitude convention must match. + +If your AI predictions NetCDF uses `0..360` longitudes but your restart uses `-180..180`, re-run Step 2 with: + +```bash +python scripts/ai_predictions_to_netcdf.py \ + --ai-predictions "$INFER_DIR/cnp_predictions" \ + --variable-list "$VAR_LIST" \ + --wrap-longitude \ + --output "$AI_PRED_NC" +``` + +For your specific restart file (`NORTHERA5`), it already uses negative longitudes, and your existing TESNORTH10PC predictions appear to use negative longitudes too—so `--wrap-longitude` is usually not needed. + +--- + +## Step 3: Update the ELM restart NetCDF + +### (A) Preview first (no file is modified) + +Update only the 5 soil-P variables that correspond to the common “5P” update: +`labilep_vr,occlp_vr,solutionp_vr,secondp_vr,primp_vr` + +```bash +python scripts/ai_predictions_to_restart.py \ + --ai-predictions "$AI_PRED_NC" \ + --restart-file "$RESTART_FILE" \ + --output "$OUT_RESTART" \ + --variable-list "$VAR_LIST" \ + --variables-to-update labilep_vr,occlp_vr,solutionp_vr,secondp_vr,primp_vr \ + --preview-only +``` + +Check the console output for: +- the variables detected as “to update” +- the coordinate mapping details + +### (B) Write the updated restart (recommended: keep a backup) + +```bash +python scripts/ai_predictions_to_restart.py \ + --ai-predictions "$AI_PRED_NC" \ + --restart-file "$RESTART_FILE" \ + --output "$OUT_RESTART" \ + --variable-list "$VAR_LIST" \ + --variables-to-update labilep_vr,occlp_vr,solutionp_vr,secondp_vr,primp_vr \ + --backup +``` + +This produces: +- `$OUT_RESTART` +- plus a backup next to `$RESTART_FILE` if `--backup` is set. + +--- + +## Step 4: Sanity checks (recommended) + +1) Confirm the updated restart exists: +```bash +ls -lah "$OUT_RESTART" +``` + +2) (Optional) Visual comparison / debugging: +- `scripts/restart_variable_plot.py` (often requires editing its `FILE_NEW` path) +- other comparison scripts under `scripts/` (e.g. scatter / profile comparisons) + +--- + +## Notes / gotchas + +- `ai_predictions_to_restart.py` updates: + - PFT1D variables in the model’s first “column mapping” slot (skips PFT0) + - soil2D variables in the first column and first 10 layers +- If your “new dataset” changes the variable set or ordering, you must keep using the correct `--variable-list` for TESNORTHERA5 (`CNP_IO_tesnorth10pct.txt`). + diff --git a/docs/INSTRUCTIONS_TRENDY_1_AI_RESTART_CREATION.md b/docs/INSTRUCTIONS_TRENDY_1_AI_RESTART_CREATION.md index 2ee8b3d..42d1f1e 100644 --- a/docs/INSTRUCTIONS_TRENDY_1_AI_RESTART_CREATION.md +++ b/docs/INSTRUCTIONS_TRENDY_1_AI_RESTART_CREATION.md @@ -9,8 +9,10 @@ This document gives step-by-step instructions to repeat the process in **docs/WO This workflow assumes you have **two trained CNP models**: 1. A **global natveg_improved-like model** with a global restart file. -2. A **tropical Phase2 P-focused model** (`phase2_pvariable_focus`-like) with a - trained checkpoint and inference outputs. +2. A **tropical-only Phase2 model** trained with the updated + `config/training_config_phase2_tropical_soilp_only.json` (emphasis on the + five soil P variables), with a trained checkpoint and inference outputs + (e.g. a run like `run_20260313_224805_phase2_tropical_soilp_amazon_africa`). If you are a **new user** and do not yet have these models, first follow Section 0 in `docs/WORKFLOW_PHASE1_PHASE2_PHASE3_RESTARTS.md` to train them @@ -34,21 +36,22 @@ with `train_cnp_model.py`. In short: export BASE_RESTART_FILE="$NATVEG_RUN_DIR/updated_restart_...your_file.nc" ``` -- Tropical P-focused model (phase2_tropical-like): +- Tropical-only Phase2 model (updated config, emphasis on 5 P variables): ```bash cd /mnt/proj-shared/AI4BGC_7xw/AI4BGC python train_cnp_model.py \ - --config config/training_config_phase2_tropical_soilp_only.json \ - --run-dir cnp_results/run_YYYYMMDD_HHMMSS_phase2_pvariable_focus_custom \ + --training-config-json config/training_config_phase2_tropical_soilp_only.json \ + --output-dir cnp_results \ + --output-dir-suffix phase2_tropical_soilp_amazon_africa \ --variable-list CNP_IO_updated9_dev_dw.txt ``` Then set: ```bash - export PHASE2_RUN_DIR="cnp_results/run_YYYYMMDD_HHMMSS_phase2_pvariable_focus_custom" + export PHASE2_RUN_DIR="cnp_results/run_YYYYMMDD_HHMMSS_phase2_tropical_soilp_amazon_africa" ``` If you already have trained runs, you can reuse them. The paths used below @@ -57,7 +60,9 @@ assume the **existing** experiments (adjust if your run IDs differ): - **Natveg run:** `cnp_results/run_20260228_214757_natveg_improved` - **Base restart file:** `cnp_results/run_20260228_214757_natveg_improved/updated_restart_CNP_IO_updated9_dev_dw_20251201_TRENDY2024_default_ICB1850CNRDCTCBC_ad_spinup.elm.r.0021-01-01-00000.nc` -- **Phase2 run:** `cnp_results/run_20260305_153217_phase2_pvariable_focus` +- **Phase2 run:** Tropical-only run with the updated Phase2 config (5 P emphasis), e.g. `cnp_results/run_20260313_224805_phase2_tropical_soilp_amazon_africa`. This is the basis run for Phase3. + +For the full **three-phase workflow** (phase1_global, phase2_tropical, phase3_tworegions): Phase2 is tropical-only with `config/training_config_phase2_tropical_soilp_only.json`; Phase3 is **two-region (Amazon + Africa) bias correction and merge** of the 5 P variables. See **docs/WORKFLOW_PHASE1_PHASE2_PHASE3_RESTARTS.md**. --- @@ -81,16 +86,18 @@ Use `TRENDY_RUN_DIR` in the steps below so the new restart and any logs live und --- -## Path A: Tropical-only inference → 4P or 5P restart (no bias/scale, or with bias/scale) +## Path A: Tropical-only inference → Phase2 restart (raw 5P in tropics) -Use this if you only have (or only want) **tropical** Phase2 predictions. +Use this if you only have (or only want) **tropical** Phase2 predictions from +the updated Phase2 run (tropical-only, 5 P emphasis). Set `PHASE2_RUN_DIR` to +your run (e.g. `cnp_results/run_20260313_224805_phase2_tropical_soilp_amazon_africa`). ### A1. Phase2 tropical-only inference (if not already done) ```bash -cd cnp_results/run_20260305_153217_phase2_pvariable_focus +cd "$PHASE2_RUN_DIR" python ../../scripts/run_inference_all.py \ - --model cnp_predictions/model.pth \ + --model cnp_model.pt \ --output-dir cnp_inference_tropical_only cd ../.. ``` @@ -99,126 +106,108 @@ cd ../.. ```bash python scripts/ai_predictions_to_netcdf.py \ - --ai-predictions cnp_results/run_20260305_153217_phase2_pvariable_focus/cnp_inference_tropical_only/cnp_predictions \ - --output cnp_results/run_20260305_153217_phase2_pvariable_focus/comparison_results/ai_predictions_tropical_only.nc + --ai-predictions "$PHASE2_RUN_DIR/cnp_inference_tropical_only/cnp_predictions" \ + --output "$PHASE2_RUN_DIR/comparison_results/ai_predictions_tropical_only.nc" \ + --variable-list CNP_IO_updated9_dev_dw.txt ``` -### A3. Update natveg restart (4 P variables, raw Phase2) +### A3. Update natveg restart (5 P variables, raw Phase2) ```bash python scripts/ai_predictions_to_restart.py \ - --ai-predictions cnp_results/run_20260305_153217_phase2_pvariable_focus/comparison_results/ai_predictions_tropical_only.nc \ - --restart-file cnp_results/run_20260228_214757_natveg_improved/updated_restart_CNP_IO_updated9_dev_dw_20251201_TRENDY2024_default_ICB1850CNRDCTCBC_ad_spinup.elm.r.0021-01-01-00000.nc \ - --output "$TRENDY_RUN_DIR/updated_restart_phase2_tropical_4p.nc" \ + --ai-predictions "$PHASE2_RUN_DIR/comparison_results/ai_predictions_tropical_only.nc" \ + --restart-file "$BASE_RESTART_FILE" \ + --output "$TRENDY_RUN_DIR/updated_restart_phase2_tropical_5P_raw.nc" \ --variable-list CNP_IO_updated9_dev_dw.txt \ - --variables-to-update occlp_vr,labilep_vr,solutionp_vr,secondp_vr \ + --variables-to-update labilep_vr,occlp_vr,solutionp_vr,secondp_vr,primp_vr \ "--tropical-lat-range=-30,30" ``` -**New restart file:** `$TRENDY_RUN_DIR/updated_restart_phase2_tropical_4p.nc` +**New restart file:** `$TRENDY_RUN_DIR/updated_restart_phase2_tropical_5P_raw.nc` -To use **5P with bias/scale** on tropical-only: first run the bias/scale step (Section B2 below) on the **tropical** CSVs (script expects `cnp_inference_entire_dataset` layout; if your tropical run writes elsewhere, either symlink or run bias/scale from a copy that matches that layout), then build a NetCDF from the bias-corrected CSVs and run `ai_predictions_to_restart.py` with `--variables-to-update labilep_vr,occlp_vr,solutionp_vr,secondp_vr,primp_vr` and `--output "$TRENDY_RUN_DIR/updated_restart_phase2_5P_bias_corrected_tropical.nc"`. +For **5P with two-region bias/scale** (Phase3): use Path B below — full-grid inference, then **Amazon and Africa bias correction separately**, then **merge** with `scripts/merge_5p_bias_corrected_amazon_africa.py`, then build NetCDF and restart. See **docs/WORKFLOW_PHASE1_PHASE2_PHASE3_RESTARTS.md** §3 Path 3B. --- -## Path B: Full-grid inference → 5P bias/scale → tropical restart (recommended for 5P bias/scale) +## Path B: Full-grid inference → two-region 5P bias correction and merge → global restart -Use this to get **bias-corrected 5P in Amazon + Africa** and a single global NetCDF. +Use this to get **bias-corrected 5P in Amazon + Africa**: apply bias/scale for +Amazon only and Africa only, **merge** with `scripts/merge_5p_bias_corrected_amazon_africa.py`, +then build a global restart from the Phase1 base. This is Phase3 of the +three-phase workflow. ### B1. Phase2 full-grid inference -From repo root. You can set the output directory with `--output-dir` (default is `cnp_inference_entire_dataset` under the current directory). Example: write outputs into the Phase2 run directory, or into your trendy run directory. - -** Output into your trendy run directory** (then bias/scale script must be pointed at this run’s paths or you copy/symlink the expected structure): +From repo root. Run full-grid inference into your **Phase2 run** (the tropical-only run with updated config, e.g. `run_20260313_224805_phase2_tropical_soilp_amazon_africa`). Set `PHASE2_RUN_DIR` to that path. ```bash -# TRENDY_RUN_DIR set in Step 0 python scripts/run_inference_all.py \ - --model cnp_results/run_20260305_153217_phase2_pvariable_focus/cnp_predictions/model.pth \ - --output-dir "$TRENDY_RUN_DIR/cnp_inference_entire_dataset" \ + --model "$PHASE2_RUN_DIR/cnp_model.pt" \ + --output-dir "$PHASE2_RUN_DIR/cnp_inference_entire_dataset" \ + --variable-list CNP_IO_updated9_dev_dw.txt \ --inference-full-grid ``` -This creates `.../cnp_predictions/` and `.../soil_2d_ground_truth/` under the given `--output-dir`. The bias/scale script (B2) expects `--run-dir` to be the run that contains `cnp_inference_entire_dataset`; if you use Option 2, set `--run-dir "$TRENDY_RUN_DIR"` in B2 and ensure the run has that subdir. +This creates `cnp_predictions/` and `soil_2d_ground_truth/` under `$PHASE2_RUN_DIR/cnp_inference_entire_dataset`. The next steps use `--run-dir "$PHASE2_RUN_DIR"`. -### B2. Apply 5P bias/scale correction (Amazon + Africa) +### B2. Apply 5P bias/scale for Amazon only, then Africa only -Use as `--run-dir` the directory that **contains** `cnp_inference_entire_dataset/` (where B1 wrote the CSVs). If you used **Option 1** in B1, that is the Phase2 run dir; if you used **Option 2**, that is `$TRENDY_RUN_DIR`. +Apply bias correction **separately** for the two regions, then merge (B3). Use `--run-dir "$PHASE2_RUN_DIR"`. -**If you used Option 1 (output under Phase2 run):** +**Amazon only:** ```bash python scripts/apply_5p_bias_scale_correction.py \ - --run-dir cnp_results/run_20260305_153217_phase2_pvariable_focus \ - --region-config-json config/training_config_two_region_five_p.json \ - --output-subdir soil_2d_predictions_5P_bias_corrected_phase2 + --run-dir "$PHASE2_RUN_DIR" \ + --region-config-json config/training_config_amazon_5p_box.json \ + --output-subdir soil_2d_predictions_5P_bias_corrected_amazon ``` -**If you used Option 2 (output under TRENDY_RUN_DIR):** +**Africa only:** ```bash python scripts/apply_5p_bias_scale_correction.py \ - --run-dir "$TRENDY_RUN_DIR" \ - --region-config-json config/training_config_two_region_five_p.json \ - --output-subdir soil_2d_predictions_5P_bias_corrected_phase2 + --run-dir "$PHASE2_RUN_DIR" \ + --region-config-json config/training_config_africa_5p_box.json \ + --output-subdir soil_2d_predictions_5P_bias_corrected_africa ``` -Corrected CSVs end up under `/cnp_inference_entire_dataset/cnp_predictions/soil_2d_predictions_5P_bias_corrected_phase2/`. Bias/scale params are written to `/analysis/bias_scale_params_5P_two_regions.json`. - -Optional: copy bias/scale params into the trendy run for provenance (if Phase2 run was used for B2, copy from there): +### B3. Merge Amazon and Africa corrected 5P ```bash -cp cnp_results/run_20260305_153217_phase2_pvariable_focus/analysis/bias_scale_params_5P_two_regions.json \ - "$TRENDY_RUN_DIR/bias_scale_params_5P_two_regions_phase2.json" +python scripts/merge_5p_bias_corrected_amazon_africa.py \ + --run-dir "$PHASE2_RUN_DIR" ``` -### B3. NetCDF from bias-corrected 5P predictions +This creates `$PHASE2_RUN_DIR/cnp_inference_entire_dataset/cnp_predictions/soil_2d_predictions_5P_bias_corrected_amazon_africa/`. -Use the same run directory you used in B2 (Phase2 run dir or `$TRENDY_RUN_DIR`). +### B4. NetCDF from merged bias-corrected 5P (amazon_africa) ```bash -# If you used Option 1 in B1/B2 (Phase2 run dir): -python scripts/ai_predictions_to_netcdf.py \ - --ai-predictions cnp_results/run_20260305_153217_phase2_pvariable_focus/cnp_inference_entire_dataset/cnp_predictions/ \ - --soil-2d-bias-corrected-subdir soil_2d_predictions_5P_bias_corrected_phase2 \ - --variable-list CNP_IO_updated9_dev_dw.txt \ - --output cnp_results/run_20260305_153217_phase2_pvariable_focus/comparison_results/ai_predictions_5P_bias_corrected_phase2.nc - -# If you used Option 2, write NetCDF into your run dir: -mkdir -p "$TRENDY_RUN_DIR/comparison_results" +mkdir -p "$PHASE2_RUN_DIR/comparison_results" python scripts/ai_predictions_to_netcdf.py \ - --ai-predictions "$TRENDY_RUN_DIR/cnp_inference_entire_dataset/cnp_predictions/" \ - --soil-2d-bias-corrected-subdir soil_2d_predictions_5P_bias_corrected_phase2 \ + --ai-predictions "$PHASE2_RUN_DIR/cnp_inference_entire_dataset/cnp_predictions/" \ + --soil-2d-bias-corrected-subdir soil_2d_predictions_5P_bias_corrected_amazon_africa \ --variable-list CNP_IO_updated9_dev_dw.txt \ - --output "$TRENDY_RUN_DIR/comparison_results/ai_predictions_5P_bias_corrected_phase2.nc" + --output "$PHASE2_RUN_DIR/comparison_results/ai_predictions_5P_bias_corrected_amazon_africa.nc" ``` -### B4. Generate new restart (5P in tropics, from natveg base) +### B5. Generate global restart (Phase1 base + 5P in tropics) -Use the path to the NetCDF you produced in B3 (Phase2 run dir or `$TRENDY_RUN_DIR/comparison_results/...`). +Use the Phase1 **base restart** (`$BASE_RESTART_FILE`) so the result is a full global restart with 5P updated only in the tropics (Amazon and Africa bias-corrected; other tropics raw Phase2). ```bash -# If Option 1 was used (NetCDF in Phase2 run dir): -python scripts/ai_predictions_to_restart.py \ - --ai-predictions cnp_results/run_20260305_153217_phase2_pvariable_focus/comparison_results/ai_predictions_5P_bias_corrected_phase2.nc \ - --restart-file cnp_results/run_20260228_214757_natveg_improved/updated_restart_CNP_IO_updated9_dev_dw_20251201_TRENDY2024_default_ICB1850CNRDCTCBC_ad_spinup.elm.r.0021-01-01-00000.nc \ - --output "$TRENDY_RUN_DIR/updated_restart_phase2_5P_bias_corrected_tropical.nc" \ - --variable-list CNP_IO_updated9_dev_dw.txt \ - --variables-to-update labilep_vr,occlp_vr,solutionp_vr,secondp_vr,primp_vr \ - "--tropical-lat-range=-30,30" - -# If Option 2 was used (NetCDF in trendy run dir): python scripts/ai_predictions_to_restart.py \ - --ai-predictions "$TRENDY_RUN_DIR/comparison_results/ai_predictions_5P_bias_corrected_phase2.nc" \ - --restart-file cnp_results/run_20260228_214757_natveg_improved/updated_restart_CNP_IO_updated9_dev_dw_20251201_TRENDY2024_default_ICB1850CNRDCTCBC_ad_spinup.elm.r.0021-01-01-00000.nc \ - --output "$TRENDY_RUN_DIR/updated_restart_phase2_5P_bias_corrected_tropical.nc" \ + --ai-predictions "$PHASE2_RUN_DIR/comparison_results/ai_predictions_5P_bias_corrected_amazon_africa.nc" \ + --restart-file "$BASE_RESTART_FILE" \ + --output "$TRENDY_RUN_DIR/updated_restart_global_5P_bias_corrected_amazon_africa.nc" \ --variable-list CNP_IO_updated9_dev_dw.txt \ --variables-to-update labilep_vr,occlp_vr,solutionp_vr,secondp_vr,primp_vr \ "--tropical-lat-range=-30,30" ``` -**New restart file:** `$TRENDY_RUN_DIR/updated_restart_phase2_5P_bias_corrected_tropical.nc` +**New restart file:** `$TRENDY_RUN_DIR/updated_restart_global_5P_bias_corrected_amazon_africa.nc` --- @@ -232,21 +221,24 @@ python scripts/ai_predictions_to_restart.py \ Example README: ```bash - echo "Restart created: $(date). Path B (full-grid + 5P bias/scale). Base: natveg_improved. Phase2: run_20260305_153217_phase2_pvariable_focus." > "$TRENDY_RUN_DIR/README_restart.txt" + echo "Restart created: $(date). Path B (two-region 5P bias correction + merge). Base: natveg_improved. Phase2: run_*_phase2_tropical_soilp_amazon_africa." > "$TRENDY_RUN_DIR/README_restart.txt" ``` 2. **Run the land model** using the new restart and validate (e.g. Amazon/Africa 5P profiles, NPP, fluxes) as in **docs/WORKFLOW_5P_TWO_REGIONS_BIAS_SCALE_AND_RESTARTS.md**. --- -## Summary checklist (Path B – 5P bias/scale) +## Summary checklist (Path B – two-region 5P bias correction and merge) - [ ] `git checkout -b trendy_1_ai_restart_creation` - [ ] Create `TRENDY_RUN_DIR=cnp_results/run__trendy_1_ai_restart_creation` -- [ ] Run Phase2 **full-grid** inference (`--inference-full-grid`) -- [ ] Run **apply_5p_bias_scale_correction.py** (Phase2 run-dir, two-region config) -- [ ] Run **ai_predictions_to_netcdf.py** on bias-corrected 5P CSVs -- [ ] Run **ai_predictions_to_restart.py** (natveg base restart → `$TRENDY_RUN_DIR/updated_restart_phase2_5P_bias_corrected_tropical.nc`, 5P, `--tropical-lat-range=-30,30`) +- [ ] Set `PHASE2_RUN_DIR` to your tropical-only Phase2 run (e.g. `run_*_phase2_tropical_soilp_amazon_africa`) +- [ ] Run Phase2 **full-grid** inference into `$PHASE2_RUN_DIR/cnp_inference_entire_dataset` +- [ ] Run **apply_5p_bias_scale_correction.py** for Amazon only (`training_config_amazon_5p_box.json`) +- [ ] Run **apply_5p_bias_scale_correction.py** for Africa only (`training_config_africa_5p_box.json`) +- [ ] Run **merge_5p_bias_corrected_amazon_africa.py** (`--run-dir "$PHASE2_RUN_DIR"`) +- [ ] Run **ai_predictions_to_netcdf.py** on merged `soil_2d_predictions_5P_bias_corrected_amazon_africa` +- [ ] Run **ai_predictions_to_restart.py** (Phase1 base → `updated_restart_global_5P_bias_corrected_amazon_africa.nc`, 5P, `--tropical-lat-range=-30,30`) - [ ] Document and validate -Reference: **docs/WORKFLOW_5P_TWO_REGIONS_BIAS_SCALE_AND_RESTARTS.md**. +Reference: **docs/WORKFLOW_PHASE1_PHASE2_PHASE3_RESTARTS.md**. diff --git a/docs/WORKFLOW_PHASE1_PHASE2_PHASE3_RESTARTS.md b/docs/WORKFLOW_PHASE1_PHASE2_PHASE3_RESTARTS.md index 81b46e6..4d52591 100644 --- a/docs/WORKFLOW_PHASE1_PHASE2_PHASE3_RESTARTS.md +++ b/docs/WORKFLOW_PHASE1_PHASE2_PHASE3_RESTARTS.md @@ -34,12 +34,11 @@ The script will: ## Workflow: phase1_global, phase2_tropical, phase3_tworegions restarts -This document summarizes how to create three restart files from existing -`natveg_improved` and `phase2_pvariable_focus` runs: +This document summarizes how to create three restart files: -- **phase1_global**: copy of the `natveg_improved` restart (baseline). -- **phase2_tropical**: tropical restart with raw Phase2 5P in \([-30°, 30°]\). -- **phase3_tworegions**: tropical restart with Phase2 5P plus Amazon + Africa bias/scale correction. +- **phase1_global**: global baseline restart (copy of the natveg_improved / Phase1 restart). +- **phase2_tropical**: tropical-only training with the **updated** `config/training_config_phase2_tropical_soilp_only.json` (emphasis on the five soil P variables). Produces a global restart with **raw** Phase2 5P in the tropical band \([-30°, 30°]\); extratropics unchanged from Phase1. The Phase2 run is the basis for Phase3 (e.g. `run_*_phase2_tropical_soilp_amazon_africa`). +- **phase3_tworegions**: **two-region (Amazon + Africa) bias correction and merge** of the 5 P variables. Apply bias/scale correction for Amazon and Africa separately, merge with `scripts/merge_5p_bias_corrected_amazon_africa.py`, then build a global restart from the Phase1 base with bias-corrected 5P in the tropics; extratropics unchanged from Phase1. It is based on: @@ -71,11 +70,10 @@ with your paths: - **Natveg run**: `cnp_results/run_20260228_214757_natveg_improved` - **Base restart**: - `cnp_results/run_20260228_214757_natveg_improved/updated_restart_CNP_IO_updated9_dev_dw_20251201_TRENDY2024_default_ICB1850CNRDCTCBC_ad_spinup.elm.r.0021-01-01-00000.nc` -- **Phase2 run**: `cnp_results/run_20260305_153217_phase2_pvariable_focus` +- **Phase2 run**: tropical-only run with the **updated** Phase2 config (emphasis on 5 P variables), e.g. `cnp_results/run_20260313_224805_phase2_tropical_soilp_amazon_africa`. This is the basis run for Phase3. - **Variable list file**: `CNP_IO_updated9_dev_dw.txt` in repo root. -If you are a new user and want to **train your own pair of models** instead of -reusing these runs, follow these steps first. +If you are a new user, train your own models as in §0.1 and §0.2. Phase2 must use **tropical-only** training with `config/training_config_phase2_tropical_soilp_only.json` (the updated config with more emphasis on the five P variables). #### 0.1 Train a natveg_improved‑like global model (for phase1_global) @@ -102,36 +100,61 @@ After training: export BASE_RESTART="$NATVEG_RUN_DIR/updated_restart_...your_file.nc" ``` -#### 0.2 Train a phase2_tropical‑like P‑focused tropical model +#### 0.2 Train the Phase2 tropical model (updated config, emphasis on 5 P variables) -From the repo root: +Phase2 uses **tropical-only** training with the **updated** +`config/training_config_phase2_tropical_soilp_only.json`, which emphasizes the +five soil P variables (labilep_vr, occlp_vr, solutionp_vr, secondp_vr, primp_vr). +The resulting run (e.g. `run_*_phase2_tropical_soilp_amazon_africa`) is the basis +for Phase3 two-region bias correction and merge. + +**Option A – use the automation script:** + +```bash +cd /mnt/proj-shared/AI4BGC_7xw/AI4BGC + +export CONFIG_TROPICAL=config/training_config_phase2_tropical_soilp_only.json +bash scripts/run_phase2_tropical.sh +``` + +The script writes to `run_*_phase2_tropical`. To match the basis run naming +(`run_*_phase2_tropical_soilp_amazon_africa`), use Option B. + +**Option B – train manually with suffix `phase2_tropical_soilp_amazon_africa`:** ```bash cd /mnt/proj-shared/AI4BGC_7xw/AI4BGC python train_cnp_model.py \ - --config config/training_config_phase2_tropical_soilp_only.json \ - --run-dir cnp_results/run_YYYYMMDD_HHMMSS_phase2_pvariable_focus_custom \ + --training-config-json config/training_config_phase2_tropical_soilp_only.json \ + --output-dir cnp_results \ + --output-dir-suffix phase2_tropical_soilp_amazon_africa \ --variable-list CNP_IO_updated9_dev_dw.txt ``` Then set: ```bash -export PHASE2_RUN_DIR="cnp_results/run_YYYYMMDD_HHMMSS_phase2_pvariable_focus_custom" +export PHASE2_RUN_DIR="cnp_results/run_YYYYMMDD_HHMMSS_phase2_tropical_soilp_amazon_africa" ``` -The remaining inference, bias/scale, and restart steps in this document are the -same; only `NATVEG_RUN_DIR`, `BASE_RESTART`, and `PHASE2_RUN_DIR` change to -point at your newly trained models. +The remaining inference, bias/scale, and restart steps are the same; only +`NATVEG_RUN_DIR`, `BASE_RESTART`, and `PHASE2_RUN_DIR` change to your run paths. --- -Create three run directories (example naming): +Create three run directories (example naming) and set run-specific variables. +**Existing users** can set `NATVEG_RUN_DIR`, `BASE_RESTART`, and `PHASE2_RUN_DIR` to +the example paths above (or to their own run paths) and skip training. ```bash cd /mnt/proj-shared/AI4BGC_7xw/AI4BGC +# For existing users: point to your natveg and Phase2 runs and base restart +# export NATVEG_RUN_DIR=cnp_results/run_20260228_214757_natveg_improved +# export BASE_RESTART="$NATVEG_RUN_DIR/updated_restart_CNP_IO_updated9_dev_dw_20251201_TRENDY2024_default_ICB1850CNRDCTCBC_ad_spinup.elm.r.0021-01-01-00000.nc" +# export PHASE2_RUN_DIR=cnp_results/run_20260313_224805_phase2_tropical_soilp_amazon_africa # tropical-only, updated Phase2 config (5 P emphasis) + TS=$(date +%Y%m%d_%H%M%S) RUN1_DIR="cnp_results/run_${TS}_phase1_global" @@ -145,11 +168,11 @@ mkdir -p "$RUN1_DIR" "$RUN2_DIR" "$RUN3_DIR" ### 1. phase1_global: global baseline restart (natveg_improved copy) -**Goal:** create a clearly labeled baseline restart identical to `natveg_improved`. +**Goal:** create a clearly labeled baseline restart identical to the natveg/Phase1 +global restart. Use `BASE_RESTART` (from Phase1 or `NATVEG_RUN_DIR`). ```bash -cp cnp_results/run_20260228_214757_natveg_improved/updated_restart_CNP_IO_updated9_dev_dw_20251201_TRENDY2024_default_ICB1850CNRDCTCBC_ad_spinup.elm.r.0021-01-01-00000.nc \ - "$RUN1_DIR/updated_restart_phase1_global_natveg_improved.nc" +cp "$BASE_RESTART" "$RUN1_DIR/updated_restart_phase1_global_natveg_improved.nc" ``` Optional README: @@ -161,18 +184,26 @@ echo "phase1_global: copy of natveg_improved restart on $(date)" \ --- -### 2. phase2_tropical: tropical restart from raw Phase2 predictions +### 2. phase2_tropical: tropical-only run (updated config) → restart with raw Phase2 5P + +**Goal:** Use the **Phase2 run** trained with the updated +`config/training_config_phase2_tropical_soilp_only.json` (tropical-only, emphasis +on the five P variables). Overwrite 5P in the tropical band \([-30°, 30°]\) with +**raw Phase2 predictions**; extratropics stay as in the Phase1 base. -**Goal:** overwrite all 5P variables in the tropical band \([-30°, 30°]\) -using **raw Phase2 predictions**. +Set `PHASE2_RUN_DIR` to your Phase2 run (e.g. +`cnp_results/run_20260313_224805_phase2_tropical_soilp_amazon_africa`). The +commands below use `$PHASE2_RUN_DIR`. #### 2.1 Tropical-only Phase2 inference (if not already done) +Use the checkpoint in your Phase2 run (`cnp_model.pt` or `cnp_predictions/model.pth`): + ```bash -cd cnp_results/run_20260305_153217_phase2_pvariable_focus +cd "$PHASE2_RUN_DIR" python ../../scripts/run_inference_all.py \ - --model cnp_predictions/model.pth \ + --model cnp_model.pt \ --output-dir cnp_inference_tropical_only cd - # back to repo root @@ -182,16 +213,17 @@ cd - # back to repo root ```bash python scripts/ai_predictions_to_netcdf.py \ - --ai-predictions cnp_results/run_20260305_153217_phase2_pvariable_focus/cnp_inference_tropical_only/cnp_predictions \ - --output cnp_results/run_20260305_153217_phase2_pvariable_focus/comparison_results/ai_predictions_tropical_only.nc + --ai-predictions "$PHASE2_RUN_DIR/cnp_inference_tropical_only/cnp_predictions" \ + --output "$PHASE2_RUN_DIR/comparison_results/ai_predictions_tropical_only.nc" \ + --variable-list CNP_IO_updated9_dev_dw.txt ``` #### 2.3 Create phase2_tropical restart (5P in tropics, raw Phase2) ```bash python scripts/ai_predictions_to_restart.py \ - --ai-predictions cnp_results/run_20260305_153217_phase2_pvariable_focus/comparison_results/ai_predictions_tropical_only.nc \ - --restart-file cnp_results/run_20260228_214757_natveg_improved/updated_restart_CNP_IO_updated9_dev_dw_20251201_TRENDY2024_default_ICB1850CNRDCTCBC_ad_spinup.elm.r.0021-01-01-00000.nc \ + --ai-predictions "$PHASE2_RUN_DIR/comparison_results/ai_predictions_tropical_only.nc" \ + --restart-file "$BASE_RESTART" \ --output "$RUN2_DIR/updated_restart_phase2_tropical_5P_raw.nc" \ --variable-list CNP_IO_updated9_dev_dw.txt \ --variables-to-update labilep_vr,occlp_vr,solutionp_vr,secondp_vr,primp_vr \ @@ -207,25 +239,118 @@ echo "phase2_tropical: natveg base + raw Phase2 5P in tropics on $(date)" \ --- -### 3. phase3_tworegions: tropical restart with two-region 5P bias/scale +### 3. phase3_tworegions: two-region bias correction and merge of 5 P variables + +**Goal:** Apply **Amazon and Africa 5P bias/scale correction** separately, **merge** +the corrected 5 P variables with `scripts/merge_5p_bias_corrected_amazon_africa.py`, +then build a **global** restart from the Phase1 base with bias-corrected 5P in the +tropics. This is the recommended flow (Path 3B, using your Phase2 run). + +**Path 3B** (recommended): Use the Phase2 run that already has (or will have) +full-grid inference; apply Amazon-only and Africa-only bias correction, merge, +then NetCDF and restart. Basis run: same as Phase2 (e.g. +`run_20260313_224805_phase2_tropical_soilp_amazon_africa`). + +**Path 3A** (alternative): Single two-region correction with inference in a +dedicated RUN3_DIR; use `scripts/run_phase3_tworegions.sh` or run steps manually. + +--- + +#### Path 3B (recommended): Amazon and Africa corrections separately, then merge + +Use your **Phase2 run** (same as in §2, e.g. +`run_20260313_224805_phase2_tropical_soilp_amazon_africa`). It must contain +full-grid inference (`cnp_inference_entire_dataset`). Apply bias/scale for +Amazon only and Africa only, merge with +`scripts/merge_5p_bias_corrected_amazon_africa.py`, then build the NetCDF and +global restart. Region configs: `config/training_config_amazon_5p_box.json` and +`config/training_config_africa_5p_box.json`. + +**3B.1 Ensure full-grid inference exists under Phase2 run** + +If not already done, run full-grid inference into your Phase2 run: + +```bash +python scripts/run_inference_all.py \ + --model "$PHASE2_RUN_DIR/cnp_model.pt" \ + --output-dir "$PHASE2_RUN_DIR/cnp_inference_entire_dataset" \ + --variable-list CNP_IO_updated9_dev_dw.txt \ + --inference-full-grid +``` + +**3B.2 Apply 5P bias/scale for Amazon only** + +```bash +python scripts/apply_5p_bias_scale_correction.py \ + --run-dir "$PHASE2_RUN_DIR" \ + --region-config-json config/training_config_amazon_5p_box.json \ + --output-subdir soil_2d_predictions_5P_bias_corrected_amazon +``` + +**3B.3 Apply 5P bias/scale for Africa only** -**Goal:** use Phase2 predictions plus **Amazon + Africa 5P bias/scale correction** -and map those corrected 5P values into the tropical band of the natveg restart. +```bash +python scripts/apply_5p_bias_scale_correction.py \ + --run-dir "$PHASE2_RUN_DIR" \ + --region-config-json config/training_config_africa_5p_box.json \ + --output-subdir soil_2d_predictions_5P_bias_corrected_africa +``` -This follows **Path B** from -`WORKFLOW_5P_TWO_REGIONS_BIAS_SCALE_AND_RESTARTS.md`, but writes outputs into -`RUN3_DIR`. +**3B.4 Merge Amazon and Africa corrected 5P** -#### 3.1 Full-grid Phase2 inference into RUN3_DIR +```bash +python scripts/merge_5p_bias_corrected_amazon_africa.py \ + --run-dir "$PHASE2_RUN_DIR" +``` + +This creates +`$PHASE2_RUN_DIR/cnp_inference_entire_dataset/cnp_predictions/soil_2d_predictions_5P_bias_corrected_amazon_africa/`. + +**3B.5 Build global NetCDF and create global restart (Phase1 base + 5P in tropics)** + +```bash +mkdir -p "$PHASE2_RUN_DIR/comparison_results" + +python scripts/ai_predictions_to_netcdf.py \ + --ai-predictions "$PHASE2_RUN_DIR/cnp_inference_entire_dataset/cnp_predictions/" \ + --soil-2d-bias-corrected-subdir soil_2d_predictions_5P_bias_corrected_amazon_africa \ + --variable-list CNP_IO_updated9_dev_dw.txt \ + --output "$PHASE2_RUN_DIR/comparison_results/ai_predictions_5P_bias_corrected_amazon_africa.nc" + +python scripts/ai_predictions_to_restart.py \ + --ai-predictions "$PHASE2_RUN_DIR/comparison_results/ai_predictions_5P_bias_corrected_amazon_africa.nc" \ + --restart-file "$BASE_RESTART" \ + --output "$RUN3_DIR/updated_restart_global_5P_bias_corrected_amazon_africa.nc" \ + --variable-list CNP_IO_updated9_dev_dw.txt \ + --variables-to-update labilep_vr,occlp_vr,solutionp_vr,secondp_vr,primp_vr \ + "--tropical-lat-range=-30,30" +``` + +For Path 3B you can write the restart into `RUN3_DIR` (as above) or into +`$PHASE2_RUN_DIR`; adjust `--output` accordingly. The **base restart must be +the Phase1 global restart** (`BASE_RESTART`) so the result is a full global +restart with 5P updated only in the tropics (Amazon and Africa use +bias-corrected values; other tropics use raw Phase2 from the merged NetCDF). + +--- + +#### Path 3A (alternative): Single two-region correction (inference in RUN3_DIR) + +Use when you want a dedicated Phase3 run directory or when using +`scripts/run_phase3_tworegions.sh`. Full-grid inference is run into `RUN3_DIR` +and one bias/scale correction is applied for both regions. + +**3A.1 Full-grid Phase2 inference into RUN3_DIR** ```bash python scripts/run_inference_all.py \ - --model cnp_results/run_20260305_153217_phase2_pvariable_focus/cnp_predictions/model.pth \ + --model "$PHASE2_RUN_DIR/cnp_model.pt" \ --output-dir "$RUN3_DIR/cnp_inference_entire_dataset" \ + --variable-list CNP_IO_updated9_dev_dw.txt \ --inference-full-grid ``` -#### 3.2 Apply 5P bias/scale correction (Amazon + Africa only) +**3A.2 Apply 5P bias/scale correction (Amazon + Africa in one call)** ```bash python scripts/apply_5p_bias_scale_correction.py \ @@ -234,14 +359,11 @@ python scripts/apply_5p_bias_scale_correction.py \ --output-subdir soil_2d_predictions_5P_bias_corrected_phase2 ``` -Outputs: - -- Corrected 5P CSVs under - `"$RUN3_DIR/cnp_inference_entire_dataset/cnp_predictions/soil_2d_predictions_5P_bias_corrected_phase2/"` -- Bias/scale parameters under - `"$RUN3_DIR/analysis/bias_scale_params_5P_two_regions.json"` +Outputs: corrected 5P CSVs under +`$RUN3_DIR/cnp_inference_entire_dataset/cnp_predictions/soil_2d_predictions_5P_bias_corrected_phase2/` +and `$RUN3_DIR/analysis/bias_scale_params_5P_two_regions.json`. -#### 3.3 Build global NetCDF from bias-corrected 5P predictions +**3A.3 Build global NetCDF and create phase3 restart** ```bash mkdir -p "$RUN3_DIR/comparison_results" @@ -251,21 +373,19 @@ python scripts/ai_predictions_to_netcdf.py \ --soil-2d-bias-corrected-subdir soil_2d_predictions_5P_bias_corrected_phase2 \ --variable-list CNP_IO_updated9_dev_dw.txt \ --output "$RUN3_DIR/comparison_results/ai_predictions_5P_bias_corrected_phase2.nc" -``` - -#### 3.4 Create phase3_tworegions restart (bias-corrected 5P in tropics) -```bash python scripts/ai_predictions_to_restart.py \ --ai-predictions "$RUN3_DIR/comparison_results/ai_predictions_5P_bias_corrected_phase2.nc" \ - --restart-file cnp_results/run_20260228_214757_natveg_improved/updated_restart_CNP_IO_updated9_dev_dw_20251201_TRENDY2024_default_ICB1850CNRDCTCBC_ad_spinup.elm.r.0021-01-01-00000.nc \ + --restart-file "$BASE_RESTART" \ --output "$RUN3_DIR/updated_restart_phase3_tworegions_5P_bias_corrected_tropical.nc" \ --variable-list CNP_IO_updated9_dev_dw.txt \ --variables-to-update labilep_vr,occlp_vr,solutionp_vr,secondp_vr,primp_vr \ "--tropical-lat-range=-30,30" ``` -Optional README: +--- + +Optional README for Phase3: ```bash echo "phase3_tworegions: natveg base + Phase2 5P with Amazon+Africa bias/scale in tropics on $(date)" \ @@ -291,16 +411,22 @@ inspecting results between phases: - Exports `NATVEG_RUN_DIR` and `BASE_RESTART` in its log output. 2. `scripts/run_phase2_tropical.sh` - - Trains the tropical Phase2 P‑focused model (or reuses an existing one). + - Trains the **tropical-only** Phase2 model with the updated + `config/training_config_phase2_tropical_soilp_only.json` (emphasis on 5 P + variables), or reuses an existing run (e.g. + `run_*_phase2_tropical_soilp_amazon_africa`). - Runs tropical‑only inference and creates the **phase2_tropical** restart with raw Phase2 5P in the tropics, using `BASE_RESTART` from phase 1. - Exports `PHASE2_RUN_DIR` and prints the path to the tropical restart. -3. `scripts/run_phase3_tworegions.sh` - - Runs full‑grid inference with the Phase2 model. - - Applies two‑region 5P bias/scale correction (Amazon + Africa). - - Creates the **phase3_tworegions** tropical restart with bias‑corrected - Phase2 5P, using the same `BASE_RESTART`. +3. **Phase3 – two-region bias correction and merge:** + - **Recommended (Path 3B):** Run the Phase3 steps in §3 manually: full-grid + inference in `PHASE2_RUN_DIR` (if needed), apply Amazon-only and + Africa-only bias correction, run `scripts/merge_5p_bias_corrected_amazon_africa.py`, + then build NetCDF and global restart from `BASE_RESTART`. + - **Alternative (Path 3A):** `scripts/run_phase3_tworegions.sh` runs + full‑grid inference into a new RUN3_DIR, applies a single two‑region + correction, and creates the phase3_tworegions restart. Between scripts you can: diff --git a/scripts/ai_predictions_to_netcdf.py b/scripts/ai_predictions_to_netcdf.py index da7104a..59521a2 100644 --- a/scripts/ai_predictions_to_netcdf.py +++ b/scripts/ai_predictions_to_netcdf.py @@ -248,9 +248,10 @@ def load_ai_predictions( if 'soil_2d' not in preds: preds['soil_2d'] = {} preds['soil2d_coords'] = {} - for p in sorted(bias_dir.glob('predictions_Y_*_bias_corrected.csv')): - # predictions_Y__bias_corrected.csv -> - var_name = p.stem.replace('predictions_Y_', '').replace('_bias_corrected', '') + for p in sorted(bias_dir.glob('predictions_Y_*_bias_corrected*.csv')): + # predictions_Y__bias_corrected.csv or ..._bias_corrected_amazon_africa.csv -> + stem = p.stem.replace('predictions_Y_', '') + var_name = stem.split('_bias_corrected')[0] df = pd.read_csv(p) lon, lat = _extract_coords(df) preds['soil2d_coords'][var_name] = (lon, lat) diff --git a/scripts/apply_5p_africa_shrinkage.py b/scripts/apply_5p_africa_shrinkage.py new file mode 100644 index 0000000..f7070a6 --- /dev/null +++ b/scripts/apply_5p_africa_shrinkage.py @@ -0,0 +1,252 @@ +#!/usr/bin/env python3 +""" +Africa-only 5P shrinkage correction for bias-corrected Phase2 predictions. + +This script reads ground truth and bias-corrected 5P prediction CSVs from a +CNP run directory, computes per-layer shrinkage factors for the Africa region, +and writes corrected CSVs where 5P values in Africa are scaled down toward the +observed range without touching Amazon or other regions. + +Expected directory structure under --run-dir: + + /cnp_inference_entire_dataset/cnp_predictions/ + soil_2d_ground_truth/ground_truth_Y_.csv + soil_2d_predictions_5P_bias_corrected_phase2/predictions_Y__bias_corrected.csv + +Output: + /cnp_inference_entire_dataset/cnp_predictions/ + soil_2d_predictions_5P_bias_corrected_africa_shrinkage/ + predictions_Y__bias_corrected_africa_shrinkage.csv + +Usage (example): + + python scripts/apply_5p_africa_shrinkage.py \ + --run-dir cnp_results/run_20260311_204845_phase3_tworegions +""" + +import argparse +import json +from pathlib import Path +from typing import Dict, List + +import numpy as np +import pandas as pd + +# Africa region box (lat_min, lat_max, lon_min, lon_max) +AFRICA_BOX = (-15.0, 15.0, 0.0, 30.0) + +# Per-variable maximum shrinkage factors (Africa only). Values in (0, 1]. +# These enforce stronger shrinking where we know overprediction is severe. +MAX_FACTOR_BY_VAR = { + "solutionp_vr": 0.2, # very aggressive shrink + "secondp_vr": 0.5, + "primp_vr": 0.5, + # occlp_vr and labilep_vr use default cap of 1.0 +} + +# Quantile used to cap corrected values relative to Africa GT (per layer). +AFRICA_GT_CAP_QUANTILE = 90.0 + +FIVE_P: List[str] = [ + "labilep_vr", + "occlp_vr", + "solutionp_vr", + "secondp_vr", + "primp_vr", +] + + +def _africa_mask(df: pd.DataFrame) -> np.ndarray: + """Return boolean mask selecting rows within the Africa lat/lon box.""" + if "Longitude" not in df.columns or "Latitude" not in df.columns: + raise ValueError("DataFrame must contain 'Longitude' and 'Latitude' columns.") + lon = df["Longitude"].to_numpy() + lat = df["Latitude"].to_numpy() + lat_min, lat_max, lon_min, lon_max = AFRICA_BOX + return (lat >= lat_min) & (lat <= lat_max) & (lon >= lon_min) & (lon <= lon_max) + + +def _per_layer_columns(df: pd.DataFrame, var: str) -> List[str]: + """Return ordered list of layer columns for a given 5P var (Y__col1_layer*).""" + prefix = f"Y_{var}_col1_layer" + cols = [c for c in df.columns if c.startswith(prefix)] + cols = sorted(cols, key=lambda c: int(c.split("layer")[-1])) + return cols + + +def compute_africa_shrinkage_factors( + gt_df: pd.DataFrame, pred_df: pd.DataFrame, var: str +) -> Dict[str, float]: + """Compute per-layer shrinkage factors for Africa for one variable. + + For each layer, we look at Africa rows where pred > 0 and gt > 0 and compute + ratio = gt / pred. The shrinkage factor is median(ratio), clamped to [0, 1] + and also to a per-variable MAX_FACTOR if provided (to avoid overly weak + shrinkage in regimes with extreme overprediction). + If there are too few valid points, we fall back to factor=1.0 (no change). + """ + if gt_df.shape != pred_df.shape: + raise ValueError(f"Shape mismatch for {var}: GT {gt_df.shape}, PRED {pred_df.shape}") + + mask_af = _africa_mask(gt_df) + layer_cols = _per_layer_columns(gt_df, var) + factors: Dict[str, float] = {} + + for col in layer_cols: + gt_vals = gt_df.loc[mask_af, col].to_numpy(dtype=float) + pred_vals = pred_df.loc[mask_af, col].to_numpy(dtype=float) + + valid = (gt_vals > 0.0) & (pred_vals > 0.0) & np.isfinite(gt_vals) & np.isfinite(pred_vals) + if valid.sum() < 10: + factors[col] = 1.0 + continue + + ratios = gt_vals[valid] / pred_vals[valid] + ratios = ratios[np.isfinite(ratios)] + if ratios.size == 0: + factors[col] = 1.0 + continue + + median_ratio = float(np.median(ratios)) + # Base clamp to [0, 1] + factor = max(0.0, min(1.0, median_ratio)) + # Apply per-variable maximum factor where defined + max_factor = MAX_FACTOR_BY_VAR.get(var) + if max_factor is not None: + factor = min(factor, max_factor) + factors[col] = factor + + return factors + + +def apply_africa_shrinkage( + gt_df: pd.DataFrame, + pred_df: pd.DataFrame, + var: str, + factors: Dict[str, float], +) -> pd.DataFrame: + """Apply per-layer shrinkage to Africa rows for one variable. + + corrected = pred * factor in Africa; outside Africa, predictions unchanged. + We also clip corrected values to [0, q_GT_Africa_layer] as a safety cap, + where q is AFRICA_GT_CAP_QUANTILE (e.g. 90th percentile). + """ + out_df = pred_df.copy() + mask_af = _africa_mask(pred_df) + layer_cols = _per_layer_columns(pred_df, var) + + for col in layer_cols: + factor = factors.get(col, 1.0) + if factor >= 0.9999: + continue + + gt_vals = gt_df.loc[mask_af, col].to_numpy(dtype=float) + if gt_vals.size == 0 or not np.isfinite(gt_vals).any(): + cap = None + else: + cap = float( + np.nanpercentile( + gt_vals[np.isfinite(gt_vals)], AFRICA_GT_CAP_QUANTILE + ) + ) + + vals = out_df.loc[mask_af, col].to_numpy(dtype=float) + corrected = vals * factor + if cap is not None and cap > 0.0: + corrected = np.minimum(corrected, cap) + corrected = np.maximum(corrected, 0.0) + out_df.loc[mask_af, col] = corrected + + return out_df + + +def main() -> None: + parser = argparse.ArgumentParser( + description=( + "Apply Africa-only shrinkage to 5P bias-corrected Phase2 predictions " + "(reduces overly large Africa 5P values while leaving other regions unchanged)." + ) + ) + parser.add_argument( + "--run-dir", + required=True, + help=( + "CNP run directory containing cnp_inference_entire_dataset/cnp_predictions/ " + "with soil_2d_ground_truth and soil_2d_predictions_5P_bias_corrected_phase2." + ), + ) + parser.add_argument( + "--output-subdir", + default="soil_2d_predictions_5P_bias_corrected_africa_shrinkage", + help=( + "Subdirectory name under cnp_inference_entire_dataset/cnp_predictions where " + "Africa-shrunk 5P CSVs will be written." + ), + ) + args = parser.parse_args() + + run_dir = Path(args.run_dir).resolve() + predictions_root = run_dir / "cnp_inference_entire_dataset" / "cnp_predictions" + gt_dir = predictions_root / "soil_2d_ground_truth" + bias_dir = predictions_root / "soil_2d_predictions_5P_bias_corrected_phase2" + out_dir = predictions_root / args.output_subdir + out_dir.mkdir(parents=True, exist_ok=True) + + if not gt_dir.is_dir(): + raise SystemExit(f"Ground-truth dir not found: {gt_dir}") + if not bias_dir.is_dir(): + raise SystemExit(f"Bias-corrected predictions dir not found: {bias_dir}") + + print(f"Run dir: {run_dir}") + print(f"GT dir: {gt_dir}") + print(f"Bias-corrected 5P dir: {bias_dir}") + print(f"Output (Africa shrinkage) dir: {out_dir}") + + all_params = {} + + for var in FIVE_P: + gt_path = gt_dir / f"ground_truth_Y_{var}.csv" + pred_path = bias_dir / f"predictions_Y_{var}_bias_corrected.csv" + if not gt_path.is_file() or not pred_path.is_file(): + print(f"Skipping {var}: missing GT or prediction CSV ({gt_path}, {pred_path})") + continue + + print(f"\nVariable: {var}") + print(f" GT: {gt_path}") + print(f" Pred: {pred_path}") + + gt_df = pd.read_csv(gt_path) + pred_df = pd.read_csv(pred_path) + if gt_df.shape != pred_df.shape: + raise SystemExit(f"Shape mismatch for {var}: GT {gt_df.shape}, PRED {pred_df.shape}") + + mask_af = _africa_mask(gt_df) + n_af = int(mask_af.sum()) + print(f" Africa rows: {n_af}") + if n_af == 0: + print(" No Africa rows found; skipping shrinkage for this variable.") + out_df = pred_df.copy() + else: + factors = compute_africa_shrinkage_factors(gt_df, pred_df, var) + out_df = apply_africa_shrinkage(gt_df, pred_df, var, factors) + all_params[var] = { + "factors": factors, + "africa_rows": n_af, + } + + out_path = out_dir / f"predictions_Y_{var}_bias_corrected_africa_shrinkage.csv" + out_df.to_csv(out_path, index=False) + print(f" Wrote Africa-shrunk predictions to {out_path}") + + if all_params: + analysis_dir = run_dir / "analysis" + analysis_dir.mkdir(parents=True, exist_ok=True) + params_path = analysis_dir / "africa_5p_shrinkage_params.json" + with params_path.open("w", encoding="utf-8") as f: + json.dump(all_params, f, indent=2) + print(f"\nSaved Africa shrinkage parameters to {params_path}") + + +if __name__ == "__main__": + main() + diff --git a/scripts/apply_5p_bias_scale_correction.py b/scripts/apply_5p_bias_scale_correction.py index a746e83..816c654 100644 --- a/scripts/apply_5p_bias_scale_correction.py +++ b/scripts/apply_5p_bias_scale_correction.py @@ -26,17 +26,31 @@ DECIMAL_PLACES_SCALED_VARS: int = 8 # keep small solutionp_vr values (e.g. 3e-6) visible # Variables where we minimize relative (percent) error so that layer-wise error is <10%. -# Uses weighted least squares with weight 1/(gt+eps)^2. -RELATIVE_ERROR_WEIGHTED_VARS: Set[str] = {"solutionp_vr"} +# Uses weighted least squares with weight 1/(gt+eps)^2. solutionp_vr and occlp_vr are critical. +RELATIVE_ERROR_WEIGHTED_VARS: Set[str] = {"solutionp_vr", "occlp_vr"} RELATIVE_ERROR_EPS: float = 1e-12 # Reference sites (lon, lat) to prioritize in relative-error fit so layer error <10% there. -REFERENCE_SITES: List[Tuple[float, float]] = [ +# We keep separate lists for Amazon and Africa so each region's fit is anchored by its own sites. +AMAZON_REFERENCE_SITES: List[Tuple[float, float]] = [ (303.75, -17.434553), # default Amazon site - (300.0, 4.240838), # Site A - (292.5, -15.549738), # Site B + (300.0, 4.240838), # Amazon Site A + (292.5, -15.549738), # Amazon Site B ] -REFERENCE_SITE_EXTRA_WEIGHT: float = 500.0 # extra weight for cells at these sites + +AFRICA_REFERENCE_SITES: List[Tuple[float, float]] = [ + (27.5, 0.471204), # nearest to (28, 0) in phase2 cnp_inference_entire_dataset + (22.5, 5.183247), # Mid-Africa + (28.0, 0.0), # Africa 28E, 0N (may not be on grid) + (20.0, -5.0), # Additional central tropical Africa site +] + +# Default combined list (used when no region-specific override is needed) +REFERENCE_SITES: List[Tuple[float, float]] = AMAZON_REFERENCE_SITES + AFRICA_REFERENCE_SITES + +REFERENCE_SITE_EXTRA_WEIGHT_AMAZON: float = 500.0 # extra weight for Amazon reference cells +REFERENCE_SITE_EXTRA_WEIGHT_AFRICA: float = 250.0 # slightly softer extra weight for Africa +REFERENCE_SITE_EXTRA_WEIGHT: float = REFERENCE_SITE_EXTRA_WEIGHT_AMAZON REFERENCE_SITE_ATOL: float = 1e-4 @@ -350,9 +364,27 @@ def compute_and_apply_corrections( scale = scale_factor_vars.get(var) if scale_factor_vars else None if scale is not None and scale != 1.0: mult_only = False # use linear fit on scaled data + + # Use relative-error objective only for variables that request it. rel_err = var in relative_error_weighted_vars - ref_lon = lon[region_mask] if rel_err and reference_sites else None - ref_lat = lat[region_mask] if rel_err and reference_sites else None + + # Choose region-specific reference sites and extra weights + if rel_err: + if region.name.startswith("amazon"): + region_ref_sites = AMAZON_REFERENCE_SITES + region_extra_weight = REFERENCE_SITE_EXTRA_WEIGHT_AMAZON + elif region.name.startswith("africa"): + region_ref_sites = AFRICA_REFERENCE_SITES + region_extra_weight = REFERENCE_SITE_EXTRA_WEIGHT_AFRICA + else: + region_ref_sites = reference_sites + region_extra_weight = REFERENCE_SITE_EXTRA_WEIGHT + else: + region_ref_sites = None + region_extra_weight = REFERENCE_SITE_EXTRA_WEIGHT + + ref_lon = lon[region_mask] if rel_err and region_ref_sites else None + ref_lat = lat[region_mask] if rel_err and region_ref_sites else None a, b = fit_bias_scale( pred_vals, gt_vals, @@ -362,12 +394,19 @@ def compute_and_apply_corrections( relative_eps=RELATIVE_ERROR_EPS, ref_site_lon=ref_lon, ref_site_lat=ref_lat, - ref_sites=reference_sites if rel_err else None, - ref_site_extra_weight=REFERENCE_SITE_EXTRA_WEIGHT, + ref_sites=region_ref_sites, + ref_site_extra_weight=region_extra_weight, ref_site_atol=REFERENCE_SITE_ATOL, ) - # For relative-error vars: if ref site still has >10% error, nudge (a,b) so ref site is exact (minimal change) - if rel_err and reference_sites and ref_lon is not None and ref_lat is not None: + # For relative-error vars: if ref site still has >10% error, nudge (a,b) so ref site is exact (minimal change). + # We only enforce this hard constraint for Amazon regions to avoid overfitting Africa to noisy GT. + if ( + rel_err + and region_ref_sites + and region.name.startswith("amazon") + and ref_lon is not None + and ref_lat is not None + ): for rlon, rlat in reference_sites: at_site = np.isclose(ref_lon, rlon, atol=REFERENCE_SITE_ATOL) & np.isclose( ref_lat, rlat, atol=REFERENCE_SITE_ATOL @@ -437,6 +476,12 @@ def compute_and_apply_corrections( raw = (a * pred_vals * scale + b) / scale else: raw = a * pred_vals + b + + # For Africa regions, apply a modest blend between raw prediction and fully corrected value + # to avoid over-correcting in noisier GT regimes. + if region.name.startswith("africa"): + alpha = 0.7 # 70% corrected, 30% original + raw = alpha * raw + (1.0 - alpha) * pred_vals # Concentrations must be non-negative corrected_df.loc[region_mask, col_name] = np.maximum(raw, 0.0) diff --git a/scripts/generate_site_5p_restart_comparison.py b/scripts/generate_site_5p_restart_comparison.py index 6a5324a..464d2a3 100644 --- a/scripts/generate_site_5p_restart_comparison.py +++ b/scripts/generate_site_5p_restart_comparison.py @@ -1,6 +1,6 @@ import argparse import os -from typing import Dict +from typing import Dict, Optional import matplotlib.pyplot as plt import numpy as np @@ -12,15 +12,10 @@ # Fixed run directories and global restart paths (same as current workflow) -PHASE1_RUN_DIR = ( - "/mnt/proj-shared/AI4BGC_7xw/AI4BGC/cnp_results/run_20260311_165538_phase1_global" -) -PHASE2_RUN_DIR = ( - "/mnt/proj-shared/AI4BGC_7xw/AI4BGC/cnp_results/run_20260311_174843_phase2_tropical" -) -PHASE3_RUN_DIR = ( - "/mnt/proj-shared/AI4BGC_7xw/AI4BGC/cnp_results/run_20260311_204845_phase3_tworegions" -) +# Default run directories (can be overridden via --phase3-run-dir for phase3) +PHASE1_RUN_DIR = "/mnt/proj-shared/AI4BGC_7xw/AI4BGC/cnp_results/run_20260315_113900_phase1_global" +PHASE2_RUN_DIR = "/mnt/proj-shared/AI4BGC_7xw/AI4BGC/cnp_results/run_20260315_175250_phase2_tropical" +PHASE3_RUN_DIR = "/mnt/proj-shared/AI4BGC_7xw/AI4BGC/cnp_results/run_20260315_210715_phase3_tworegions" PHASE1_GLOBAL_RESTART = os.path.join(PHASE1_RUN_DIR, "updated_restart_base.nc") PHASE2_GLOBAL_RESTART = os.path.join( @@ -107,9 +102,18 @@ def _load_profiles_from_restart(path: str) -> Dict[str, np.ndarray]: return profs -def generate_site_restarts(lon: float, lat: float, site_name: str) -> Dict[str, str]: +def generate_site_restarts( + lon: float, + lat: float, + site_name: str, + phase3_run_dir: Optional[str] = None, +) -> Dict[str, str]: """Extract phase1/phase2/phase3 single-point restarts for this site.""" site_id = _safe_site_name(site_name, lon, lat) + p3_dir = phase3_run_dir if phase3_run_dir else PHASE3_RUN_DIR + phase3_restart = os.path.join( + p3_dir, "updated_restart_phase3_tworegions_5P_bias_corrected_tropical.nc" + ) out_paths = { "phase1_global": os.path.join(PHASE1_RUN_DIR, f"{site_id}_phase1_global_restart.nc"), @@ -117,7 +121,7 @@ def generate_site_restarts(lon: float, lat: float, site_name: str) -> Dict[str, PHASE2_RUN_DIR, f"{site_id}_phase2_tropical_restart.nc" ), "phase3_tworegions": os.path.join( - PHASE3_RUN_DIR, f"{site_id}_phase3_tworegions_restart.nc" + p3_dir, f"{site_id}_phase3_tworegions_restart.nc" ), } @@ -135,7 +139,7 @@ def generate_site_restarts(lon: float, lat: float, site_name: str) -> Dict[str, target_lon=lon, ) extract_single_point_elm( - source_nc=PHASE3_GLOBAL_RESTART, + source_nc=phase3_restart, output_nc=out_paths["phase3_tworegions"], target_lat=lat, target_lon=lon, @@ -275,6 +279,15 @@ def main() -> None: default="site", help="Optional human-readable site name used in filenames and plot titles.", ) + parser.add_argument( + "--phase3-run-dir", + type=str, + default="", + help=( + "Phase3 run directory (e.g. cnp_results/run_YYYYMMDD_HHMMSS_phase3_tworegions). " + "If set, restarts and plots use this run instead of the script default." + ), + ) parser.add_argument( "--output-dir", type=str, @@ -290,12 +303,13 @@ def main() -> None: lat = args.lat site_name = args.site_name site_id = _safe_site_name(site_name, lon, lat) + phase3_run_dir = os.path.abspath(args.phase3_run_dir) if args.phase3_run_dir else None if args.output_dir: output_dir = os.path.abspath(args.output_dir) else: output_dir = os.path.join( - PHASE3_RUN_DIR, + phase3_run_dir if phase3_run_dir else PHASE3_RUN_DIR, "analysis", f"{site_id}_5p_restart_comparison", ) @@ -303,7 +317,9 @@ def main() -> None: print(f"Site: {site_id} (lon={lon}, lat={lat})") print(f"Output directory: {output_dir}") - restart_paths = generate_site_restarts(lon=lon, lat=lat, site_name=site_name) + restart_paths = generate_site_restarts( + lon=lon, lat=lat, site_name=site_name, phase3_run_dir=phase3_run_dir + ) compare_5p_for_site( lon=lon, lat=lat, diff --git a/scripts/merge_5p_bias_corrected_amazon_africa.py b/scripts/merge_5p_bias_corrected_amazon_africa.py new file mode 100644 index 0000000..a2eabc0 --- /dev/null +++ b/scripts/merge_5p_bias_corrected_amazon_africa.py @@ -0,0 +1,176 @@ +#!/usr/bin/env python3 +""" +Merge Amazon-only and Africa-only 5P bias-corrected predictions into one set. + +Given a Phase2 run directory with: + + /cnp_inference_entire_dataset/cnp_predictions/ + soil_2d_predictions/ (raw Phase2 predictions) + soil_2d_predictions_5P_bias_corrected_amazon/ (Amazon-only correction) + soil_2d_predictions_5P_bias_corrected_africa/ (Africa-only correction) + +this script creates: + + soil_2d_predictions_5P_bias_corrected_amazon_africa/ + +such that for each 5P variable and grid cell: + - If (lat, lon) in Amazon box [-30, 10] x [270, 330], + use the Amazon-corrected value. + - Else if (lat, lon) in Africa box [-15, 15] x [0, 30], + use the Africa-corrected value. + - Else: + keep the raw Phase2 prediction. + +Usage (example): + + python scripts/merge_5p_bias_corrected_amazon_africa.py \\ + --run-dir cnp_results/run_20260305_153217_phase2_pvariable_focus +""" + +import argparse +from pathlib import Path +from typing import List + +import numpy as np +import pandas as pd + +FIVE_P: List[str] = [ + "labilep_vr", + "occlp_vr", + "solutionp_vr", + "secondp_vr", + "primp_vr", +] + +# Boxes: (lat_min, lat_max, lon_min, lon_max) +AMAZON_BOX = (-30.0, 10.0, 270.0, 330.0) +AFRICA_BOX = (-15.0, 15.0, 0.0, 30.0) + + +def _region_mask(df: pd.DataFrame, box) -> np.ndarray: + """Return mask for rows inside a given (lat_min, lat_max, lon_min, lon_max) box.""" + lat_min, lat_max, lon_min, lon_max = box + if "Longitude" not in df.columns or "Latitude" not in df.columns: + raise ValueError("CSV must contain 'Longitude' and 'Latitude' columns.") + lon = df["Longitude"].to_numpy() + lat = df["Latitude"].to_numpy() + return (lat >= lat_min) & (lat <= lat_max) & (lon >= lon_min) & (lon <= lon_max) + + +def _per_layer_columns(df: pd.DataFrame, var: str) -> List[str]: + """Return ordered list of layer columns for a given 5P var (Y__col1_layer*).""" + prefix = f"Y_{var}_col1_layer" + cols = [c for c in df.columns if c.startswith(prefix)] + return sorted(cols, key=lambda c: int(c.split("layer")[-1])) + + +def main() -> None: + parser = argparse.ArgumentParser( + description=( + "Merge Amazon-only and Africa-only 5P bias-corrected CSVs into a single " + "prediction set: Amazon box uses Amazon-corrected values, Africa box uses " + "Africa-corrected, elsewhere keep raw Phase2 predictions." + ) + ) + parser.add_argument( + "--run-dir", + required=True, + help="Phase2 run directory containing cnp_inference_entire_dataset/cnp_predictions.", + ) + parser.add_argument( + "--raw-subdir", + default="soil_2d_predictions", + help="Subdir under cnp_predictions with raw Phase2 soil2d predictions.", + ) + parser.add_argument( + "--amazon-subdir", + default="soil_2d_predictions_5P_bias_corrected_amazon", + help="Subdir with Amazon-only bias-corrected 5P CSVs.", + ) + parser.add_argument( + "--africa-subdir", + default="soil_2d_predictions_5P_bias_corrected_africa", + help="Subdir with Africa-only bias-corrected 5P CSVs.", + ) + parser.add_argument( + "--output-subdir", + default="soil_2d_predictions_5P_bias_corrected_amazon_africa", + help="Output subdir for merged 5P CSVs.", + ) + args = parser.parse_args() + + run_dir = Path(args.run_dir).resolve() + pred_root = run_dir / "cnp_inference_entire_dataset" / "cnp_predictions" + raw_dir = pred_root / args.raw_subdir + amazon_dir = pred_root / args.amazon_subdir + africa_dir = pred_root / args.africa_subdir + out_dir = pred_root / args.output_subdir + + if not raw_dir.is_dir(): + raise SystemExit(f"Raw predictions dir not found: {raw_dir}") + if not amazon_dir.is_dir(): + raise SystemExit(f"Amazon-corrected dir not found: {amazon_dir}") + if not africa_dir.is_dir(): + raise SystemExit(f"Africa-corrected dir not found: {africa_dir}") + + out_dir.mkdir(parents=True, exist_ok=True) + + print(f"Run dir: {run_dir}") + print(f"Raw dir: {raw_dir}") + print(f"Amazon-corrected dir: {amazon_dir}") + print(f"Africa-corrected dir: {africa_dir}") + print(f"Output dir: {out_dir}") + + for var in FIVE_P: + raw_path = raw_dir / f"predictions_Y_{var}.csv" + amz_path = amazon_dir / f"predictions_Y_{var}_bias_corrected.csv" + afr_path = africa_dir / f"predictions_Y_{var}_bias_corrected.csv" + if not raw_path.is_file() or not amz_path.is_file() or not afr_path.is_file(): + print(f"Skipping {var}: missing CSV ({raw_path}, {amz_path}, {afr_path})") + continue + + print(f"\nVariable: {var}") + print(f" Raw: {raw_path}") + print(f" Amazon: {amz_path}") + print(f" Africa: {afr_path}") + + raw_df = pd.read_csv(raw_path) + amz_df = pd.read_csv(amz_path) + afr_df = pd.read_csv(afr_path) + + if raw_df.shape != amz_df.shape or raw_df.shape != afr_df.shape: + raise SystemExit( + f"Shape mismatch for {var}: raw {raw_df.shape}, amazon {amz_df.shape}, africa {afr_df.shape}" + ) + + # Sanity check coordinates + for name, df in [("amazon", amz_df), ("africa", afr_df)]: + if not np.allclose(raw_df["Longitude"], df["Longitude"]) or not np.allclose( + raw_df["Latitude"], df["Latitude"] + ): + raise SystemExit(f"Longitude/Latitude mismatch between raw and {name} for {var}") + + mask_amz = _region_mask(raw_df, AMAZON_BOX) + mask_afr = _region_mask(raw_df, AFRICA_BOX) + + print(f" Amazon rows: {int(mask_amz.sum())}") + print(f" Africa rows: {int(mask_afr.sum())}") + + out_df = raw_df.copy() + layer_cols = _per_layer_columns(raw_df, var) + + # Overwrite Amazon region from amazon-corrected + for col in layer_cols: + out_df.loc[mask_amz, col] = amz_df.loc[mask_amz, col].to_numpy() + # Overwrite Africa region from africa-corrected + for col in layer_cols: + out_df.loc[mask_afr, col] = afr_df.loc[mask_afr, col].to_numpy() + + out_path = out_dir / f"predictions_Y_{var}_bias_corrected_amazon_africa.csv" + out_df.to_csv(out_path, index=False) + print(f" Wrote merged predictions to {out_path}") + + +if __name__ == "__main__": + main() + diff --git a/scripts/run_phase1_global.sh b/scripts/run_phase1_global.sh index bac08f1..98cf62a 100644 --- a/scripts/run_phase1_global.sh +++ b/scripts/run_phase1_global.sh @@ -16,6 +16,9 @@ set -euo pipefail # - CONFIG_GLOBAL training config JSON for global model # - VARIABLE_LIST CNP_IO variable list file # - RESTART_TEMPLATE path to original ELM restart .nc (required) + +# e.g. /mnt/proj-shared/AI4BGC_7xw/AI4BGC/ELM_data/ +# 20251201_TRENDY2024_default_ICB1850CNRDCTCBC_ad_spinup.elm.r.0021-01-01-00000.nc # # Outputs: # - NATVEG_RUN_DIR (exported): global model run directory diff --git a/scripts/run_phase3_tworegions.sh b/scripts/run_phase3_tworegions.sh index b953799..a0423a7 100644 --- a/scripts/run_phase3_tworegions.sh +++ b/scripts/run_phase3_tworegions.sh @@ -6,7 +6,8 @@ set -euo pipefail # # Phase 3 for new users: # - Full-grid inference with Phase2 model -# - Apply 5P bias/scale correction (Amazon + Africa) +# - Apply 5P bias/scale correction: Amazon-only, then Africa-only, then merge +# (same workflow as run_20260313_224805_phase2_tropical_soilp_amazon_africa) # - Create a tropical restart with bias-corrected Phase2 5P # # Requirements: @@ -17,9 +18,12 @@ set -euo pipefail # - TS timestamp suffix (optional; only affects output dir name) # # Usage: -# export PHASE2_RUN_DIR=... +# export PHASE2_RUN_DIR=... # optional; default prefers run_*_phase2_tropical_soilp* then latest phase2_tropical # export BASE_RESTART=... # bash scripts/run_phase3_tworegions.sh +# +# For best Africa/Amazon 5P, use: PHASE2_RUN_DIR=cnp_results/run_20260313_224805_phase2_tropical_soilp_amazon_africa +# See docs/PHASE3_TWOREGIONS_VS_GOOD_RUN_DIFFERENCE.md for why phase3 results can differ from that run. ############################################################################### ROOT_DIR="$(cd "$(dirname "${BASH_SOURCE[0]}")/.." && pwd)" @@ -31,9 +35,13 @@ if [[ -z "${PHASE2_RUN_DIR:-}" ]] || [[ -z "${BASE_RESTART:-}" ]]; then local d="$1" [[ -f "$d/cnp_model.pt" ]] || [[ -f "$d/cnp_predictions/model.pth" ]] } + # Prefer phase2_tropical_soilp* (e.g. phase2_tropical_soilp_amazon_africa) for best Africa/Amazon 5P. _latest_phase2() { local d - for d in $(ls -td cnp_results/run_*_phase2_tropical 2>/dev/null); do + for d in $(ls -td cnp_results/run_*_phase2_tropical_soilp* 2>/dev/null); do + [[ -d "$d" ]] && _has_model "$d" && { echo "$d"; return 0; } + done + for d in $(ls -td cnp_results/run_*_phase2_tropical* 2>/dev/null); do [[ -d "$d" ]] && _has_model "$d" && { echo "$d"; return 0; } done return 1 @@ -76,19 +84,29 @@ if [[ ! -d "$INF_FULL_DIR/cnp_predictions" ]]; then --inference-full-grid fi -echo "Applying 5P bias/scale correction (Amazon + Africa)..." +# Match good-run workflow: Amazon-only correction, Africa-only correction, then merge. +echo "Applying 5P bias/scale correction (Amazon only)..." python scripts/apply_5p_bias_scale_correction.py \ --run-dir "$RUN3_DIR" \ - --region-config-json config/training_config_two_region_five_p.json \ - --output-subdir soil_2d_predictions_5P_bias_corrected_phase2 + --region-config-json config/training_config_amazon_5p_box.json \ + --output-subdir soil_2d_predictions_5P_bias_corrected_amazon + +echo "Applying 5P bias/scale correction (Africa only)..." +python scripts/apply_5p_bias_scale_correction.py \ + --run-dir "$RUN3_DIR" \ + --region-config-json config/training_config_africa_5p_box.json \ + --output-subdir soil_2d_predictions_5P_bias_corrected_africa + +echo "Merging Amazon and Africa bias-corrected 5P..." +python scripts/merge_5p_bias_corrected_amazon_africa.py --run-dir "$RUN3_DIR" mkdir -p "$RUN3_DIR/comparison_results" -BC_NETCDF="$RUN3_DIR/comparison_results/ai_predictions_5P_bias_corrected_phase2.nc" +BC_NETCDF="$RUN3_DIR/comparison_results/ai_predictions_5P_bias_corrected_amazon_africa.nc" echo "Converting bias-corrected predictions to NetCDF..." python scripts/ai_predictions_to_netcdf.py \ --ai-predictions "$INF_FULL_DIR/cnp_predictions/" \ - --soil-2d-bias-corrected-subdir soil_2d_predictions_5P_bias_corrected_phase2 \ + --soil-2d-bias-corrected-subdir soil_2d_predictions_5P_bias_corrected_amazon_africa \ --variable-list "$VARIABLE_LIST" \ --output "$BC_NETCDF" @@ -101,7 +119,7 @@ python scripts/ai_predictions_to_restart.py \ --variables-to-update labilep_vr,occlp_vr,solutionp_vr,secondp_vr,primp_vr \ "--tropical-lat-range=-30,30" -echo "phase3_tworegions: natveg base + Phase2 5P with Amazon+Africa bias/scale in tropics on $(date)" > "$RUN3_DIR/README_phase3_tworegions.txt" +echo "phase3_tworegions: natveg base + Phase2 5P with Amazon-only + Africa-only bias correction merged (same workflow as good run) in tropics on $(date)" > "$RUN3_DIR/README_phase3_tworegions.txt" echo From 68373abd56decf45bc2a18c1963df183f34bc29e Mon Sep 17 00:00:00 2001 From: Dali Wang Date: Thu, 28 May 2026 14:02:20 -0400 Subject: [PATCH 2/5] Add comprehensive report on Phase3 five-P bias correction - Introduced `REPORT_PHASE3_5P_BIAS_CORRECTION_REVIEW.md` detailing the bias correction process for soil phosphorus in the Amazon and Africa regions. - Summarized findings from empirical comparisons between older and newer Phase3 outputs, highlighting improvements and trade-offs in model performance. - Documented related workflows and scripts to enhance user understanding of the bias correction methodology and its implementation. - Provided recommendations for optimizing the bias correction process to achieve more uniform quality across regions. --- ...REPORT_PHASE3_5P_BIAS_CORRECTION_REVIEW.md | 156 +++++++++ .../compare_5p_gt_two_regions_inference.py | 300 ++++++++++++++++++ 2 files changed, 456 insertions(+) create mode 100644 docs/REPORT_PHASE3_5P_BIAS_CORRECTION_REVIEW.md create mode 100644 scripts/compare_5p_gt_two_regions_inference.py diff --git a/docs/REPORT_PHASE3_5P_BIAS_CORRECTION_REVIEW.md b/docs/REPORT_PHASE3_5P_BIAS_CORRECTION_REVIEW.md new file mode 100644 index 0000000..1824e4d --- /dev/null +++ b/docs/REPORT_PHASE3_5P_BIAS_CORRECTION_REVIEW.md @@ -0,0 +1,156 @@ +# Phase3 five-P bias correction: review and recommendations + +This report summarizes how the Phase3 soil phosphorus bias / scale correction works in this repository, what we observed when comparing older and newer Phase3 outputs, and concrete directions to improve **`solutionp_vr`** and **`occlp_vr`** across the **entire** Amazon and Africa target boxes—not only at a small set of hand-picked reference coordinates. + +**Related workflow documentation** + +- [WORKFLOW_5P_TWO_REGIONS_BIAS_SCALE_AND_RESTARTS.md](./WORKFLOW_5P_TWO_REGIONS_BIAS_SCALE_AND_RESTARTS.md) — end-to-end 5P bias/scale and restart steps. +- [WORKFLOW_PHASE1_PHASE2_PHASE3_RESTARTS.md](./WORKFLOW_PHASE1_PHASE2_PHASE3_RESTARTS.md) — Phase1 / Phase2 / Phase3 restart orchestration. +- [PHASE3_TWOREGIONS_VS_GOOD_RUN_DIFFERENCE.md](./PHASE3_TWOREGIONS_VS_GOOD_RUN_DIFFERENCE.md) — Phase3 aligned with Amazon-only → Africa-only → merge workflow. +- [SITE_5P_RESTART_COMPARISON.md](./SITE_5P_RESTART_COMPARISON.md) — site-level 5P profile plots (`generate_site_5p_restart_comparison.py`). +- [CNP_pipeline_runbook.md](./CNP_pipeline_runbook.md) — inference, validation plots, and quality reports. + +**Scripts** + +- `scripts/apply_5p_bias_scale_correction.py` — fits and applies per-layer affine maps \(\mathrm{GT} \approx a\cdot\mathrm{pred} + b\) inside each region box. +- `scripts/merge_5p_bias_corrected_amazon_africa.py` — merges **separate** Amazon-only and Africa-only corrected CSVs so each geographic box uses its own calibration. +- `scripts/compare_5p_gt_two_regions_inference.py` — `gt-summary` and `pred-eval` for Amazon/Africa boxes (median RMSE across layers per cell, multi-run CSVs). + +--- + +## 1. Executive summary + +- The pipeline uses **one affine map per variable, per soil layer, per region** (Amazon box vs Africa box), fitted on **all gridcells** falling in that box. That is already an **overall regional** fit in the sense of pooling every cell in the box. + +- For **`solutionp_vr`** and **`occlp_vr`**, the fitter switches to a **relative-error–weighted** least squares objective and can **up-weight** a fixed list of **reference lon/lat pairs**. Those sites were chosen from visualization convenience; they are **not** a substitute for optimizing error over the whole region, and they can **skew** the fit if weights are large or if sites do not lie exactly on the inference grid. + +- A **post-fit “snap”** toward reference sites (forcing the corrected value toward GT at selected cells) is applied today for **Amazon only** for the relative-error variables, not symmetrically for Africa. That asymmetry conflicts with the product goal of **near-uniform quality in both regions**. + +- Empirical comparison (see Section 3) showed that **merging independent Amazon and Africa corrections** (newer Phase3 style) fixes severe **Africa `occlp_vr`** failure modes seen when Africa inherited a correction dominated elsewhere, at the cost of small regressions on some other variables in Amazon. That tradeoff is structurally plausible for a **single** \((a,b)\) per layer per box. + +- **Ideal direction** for “almost perfect” **`solutionp_vr`** / **`occlp_vr`** on the **full** boxes: prioritize **region-wide objectives** (e.g. minimize median or Huber loss over all cells, explicit tail control), **reduce or remove dependence** on a few reference coordinates unless they are replaced by **grid-anchored** targets, and add **stability constraints** (weight caps, log-space, or sub-regions) rather than relying on five visual picks. + +--- + +## 2. What the code does today + +### 2.1 Regional affine calibration (`apply_5p_bias_scale_correction.py`) + +For each of the five variables (`labilep_vr`, `occlp_vr`, `solutionp_vr`, `secondp_vr`, `primp_vr`), and each soil layer column `Y__col1_layer{1..10}`: + +1. Select all inference rows whose \((\mathrm{lat}, \mathrm{lon})\) lie in the active **region box** (from the training JSON’s `data_filtering_config.region_boxes`, optionally splitting Amazon into south/north by latitude). + +2. Fit **one** scalar \(a\) and intercept \(b\) so that \(\mathrm{GT} \approx a\cdot\mathrm{pred} + b\) on those rows. + +3. Special cases already in code: + + - **`solutionp_vr`**: optional internal scaling (`SCALE_FACTOR_VARS`) for numerical stability when fitting small concentrations. + + - **`solutionp_vr`** and **`occlp_vr`**: **`RELATIVE_ERROR_WEIGHTED_VARS`** — weights \(w_i \propto 1 / (\mathrm{GT}_i + \varepsilon)^2\) in the weighted normal equations, plus **extra weight** at configured reference sites. + +4. After fitting, for relative-error variables in **Amazon-shaped** region names only, an additional loop can **nudge** \((a,b)\) so a reference cell is closer to exact match, subject to a relative error cap. + +The asymmetry (Amazon-only nudge) is called out in-code as avoiding “overfitting Africa to noisy GT”; the side effect is that **Africa does not receive the same hard anchoring** for the two most sensitive variables. + +### 2.2 Merge of two regional products (`merge_5p_bias_corrected_amazon_africa.py`) + +The merged product assigns: + +- Amazon box cells → predictions from the **Amazon-only** bias correction run. + +- Africa box cells → predictions from the **Africa-only** bias correction run. + +- All other cells → **raw** Phase2 neural-network predictions. + +This preserves **two independent regional fits** and avoids a single compromised regression over Amazon \(\cup\) Africa. + +--- + +## 3. Empirical context (Phase3 old vs new) + +Structured evaluation used `scripts/compare_5p_gt_two_regions_inference.py` `pred-eval` on the Amazon and Africa boxes (same definitions as the merge script: Amazon lat \([-30,10]\), lon \([270,330]\); Africa lat \([-15,15]\), lon \([0,30]\), lon in \(0\)–\(360^\circ\)). + +Summary of outcomes (median RMSE across ten layers per gridcell): + +- **Older Phase3** (`run_20260315_202304_phase3_tworegions`, `soil_2d_predictions_5P_bias_corrected_phase2`): **Africa `occlp_vr`** showed catastrophic bulk error (median RMSE on the order of \(10^2\)), consistent with restart shock in African simulations. + +- **Newer Phase3** (`run_20260315_210715_phase3_tworegions`, merged `soil_2d_predictions_5P_bias_corrected_amazon_africa`): **Africa `occlp_vr`** median RMSE improved by roughly two orders of magnitude, with modest **Amazon** regressions on some variables (e.g. `primp_vr`) and small **solutionp_vr** median changes—aligned with sacrificing a little Amazon quality to rescue Africa. + +Detailed tables and CSVs live under the newer run’s `analysis/` directory (for example `5p_phase3_old_vs_new_pred_eval_summary.csv` and `5p_phase3_old_vs_new_report.txt`). + +--- + +## 4. Reference sites vs “entire Amazon and Africa” fitting + +### 4.1 What the reference list is doing today + +In `apply_5p_bias_scale_correction.py`, **`AMAZON_REFERENCE_SITES`** and **`AFRICA_REFERENCE_SITES`** are used to: + +1. **Inflate weights** in the weighted least squares fit for **`solutionp_vr`** and **`occlp_vr`** at rows whose lon/lat match a site within a tight tolerance (`REFERENCE_SITE_ATOL`). + +2. Optionally drive the **Amazon-only** post-fit projection toward a match at one of those coordinates. + +The coordinates were chosen to align with **site plots** and manual inspection. That is useful for **debugging** and for **forcing** a known anchor, but it is **not** the same as defining the objective “minimize error over every cell in the box.” + +### 4.2 Why this can conflict with region-wide goals + +- **Coverage**: Reference sites may not coincide with inference grid nodes; tight `isclose` checks can **miss** the intended cell, so extra weight never applies. + +- **Optimality**: Heavy weight on a few points **reweights** the global minimum of the pooled loss. If the goal is **uniform** accuracy, the optimum is closer to **unweighted** or **percentile-balanced** losses over **all** regional cells (optionally with **capped** weights so tiny GT does not dominate). + +- **Asymmetry**: Hard projection for Amazon but not Africa biases the pipeline toward “Amazon looks good on the map” while Africa still depends only on the weighted LS—misaligned with **equal** standards for **`solutionp_vr`** and **`occlp_vr`**. + +### 4.3 What “overall fitting across the entire regions” should mean in practice + +Operationally, for each region \(R \in \{\mathrm{Amazon}, \mathrm{Africa}\}\) and each critical variable \(v \in \{\mathrm{solutionp\_vr}, \mathrm{occlp\_vr}\}\), a region-wide objective could be stated as: + +- Minimize a **robust** scalar over all cells \(i \in R\) and layers \(\ell\), e.g. + \(\sum_{i,\ell} \rho\bigl((a_\ell p_{i,\ell} + b_\ell) - g_{i,\ell}\bigr)\) + where \(\rho\) is Huber or a **clipped** relative loss, **without** concentrating mass on five lon/lat pairs unless those pairs are formally part of the specification. + +- Report **distribution** of errors (median, p90, worst 1%) per region after fit—not only at reference locations. + +Reference locations can remain as **diagnostics** or as **soft** priors with **small** weight, but they should not dominate the objective if the stated goal is **whole-box** performance. + +--- + +## 5. Recommendations (prioritized) + +1. **Treat `solutionp_vr` and `occlp_vr` with symmetric, region-wide objectives** + Remove or drastically lower reference-site **extra weight** for these variables, **or** replace “exact lon/lat” with **nearest inference gridcell** in the region so anchors always exist. If hard projection is kept, apply **the same rule in Africa as in Amazon** for these two variables only. + +2. **Stabilize relative-error fitting for small concentrations** + Cap weights \(1/(\mathrm{GT}+\varepsilon)^2\), blend with absolute-error terms below a GT threshold, or fit in **log-space** with back-transform so a thin tail of near-zero GT does not dictate \((a,b)\) for the whole layer. + +3. **Keep separate regional fits + merge** + Continue using **independent** Amazon and Africa calibration followed by `merge_5p_bias_corrected_amazon_africa.py`. Avoid a single pooled \((a,b)\) over both boxes unless you introduce explicit spatial features (sub-boxes, covariates). + +4. **Richer within-region structure (if affine is not enough)** + Subdivide Amazon (north/south is already supported via `split_first_region_lat` in `load_region_boxes`) or add a **shallow residual** model of corrections vs static soil / PFT predictors, still fit only on withheld or training-consistent data to avoid double-dipping if that becomes a concern. + +5. **Physical post-processing** + Enforce non-negativity and plausible bounds after correction; optionally enforce mild **monotonicity across depth** where ELM semantics support it, to reduce restart instability when the affine map wobbles layer-to-layer. + +6. **Evaluation loop** + Use `compare_5p_gt_two_regions_inference.py pred-eval` routinely on **median, mean, p90, and tail** (e.g. 99th percentile RMSE) for **`solutionp_vr`** and **`occlp_vr`** in each box, not only at reference sites. + +--- + +## 6. Reproducibility pointers + +Workflow cross-references: [WORKFLOW_5P_TWO_REGIONS_BIAS_SCALE_AND_RESTARTS.md](./WORKFLOW_5P_TWO_REGIONS_BIAS_SCALE_AND_RESTARTS.md) (Section 4: apply + merge), [WORKFLOW_PHASE1_PHASE2_PHASE3_RESTARTS.md](./WORKFLOW_PHASE1_PHASE2_PHASE3_RESTARTS.md), [PHASE3_TWOREGIONS_VS_GOOD_RUN_DIFFERENCE.md](./PHASE3_TWOREGIONS_VS_GOOD_RUN_DIFFERENCE.md). + +- Bias application: `scripts/apply_5p_bias_scale_correction.py` +- Merge: `scripts/merge_5p_bias_corrected_amazon_africa.py` +- Regional pred vs GT metrics: `scripts/compare_5p_gt_two_regions_inference.py pred-eval` + +Example Phase3 comparison artifacts (paths on the shared filesystem used in this project): + +- `cnp_results/run_20260315_210715_phase3_tworegions/analysis/5p_phase3_old_vs_new_runs.json` +- `cnp_results/run_20260315_210715_phase3_tworegions/analysis/5p_phase3_old_vs_new_report.txt` + +--- + +## Document history + +- **2026-05-11**: Initial report from code review and Phase3 empirical comparison; emphasis on full-box objectives for `solutionp_vr` / `occlp_vr` and on reducing over-reliance on visualization-driven reference sites. diff --git a/scripts/compare_5p_gt_two_regions_inference.py b/scripts/compare_5p_gt_two_regions_inference.py new file mode 100644 index 0000000..85dcf77 --- /dev/null +++ b/scripts/compare_5p_gt_two_regions_inference.py @@ -0,0 +1,300 @@ +#!/usr/bin/env python3 +""" +Compare ground truth of the five soil P variables between the Amazon and Africa +boxes used in the Phase3 / 5P bias workflow. + +Region bounds are read from the same JSON configs as apply_5p_bias_scale_correction +(default: config/training_config_amazon_5p_box.json and +config/training_config_africa_5p_box.json), i.e. the same definitions as +scripts/merge_5p_bias_corrected_amazon_africa.py: + + Amazon: lat [-30, 10], lon [270, 330] (degrees, lon 0–360) + Africa: lat [-15, 15], lon [0, 30] + +Ground truth CSVs are read from: + + //cnp_predictions/soil_2d_ground_truth/ground_truth_Y_.csv + +Each file must have Longitude, Latitude and layer columns Y__col1_layer{1..10}. + +Usage (from repo root): + + python scripts/compare_5p_gt_two_regions_inference.py \\ + --run-dir cnp_results/run_20260315_202304_phase3_tworegions + + python scripts/compare_5p_gt_two_regions_inference.py \\ + --run-dir cnp_results/run_20260315_113900_phase1_global \\ + --output-csv cnp_results/run_20260315_113900_phase1_global/analysis/5p_gt_amazon_vs_africa.csv +""" + +from __future__ import annotations + +import argparse +import json +import sys +from dataclasses import dataclass +from pathlib import Path +from typing import Dict, List, Optional, Tuple + +import numpy as np +import pandas as pd + +REPO_ROOT = Path(__file__).resolve().parents[1] + +FIVE_P_VARS: List[str] = [ + "labilep_vr", + "occlp_vr", + "solutionp_vr", + "secondp_vr", + "primp_vr", +] + + +@dataclass +class RegionBox: + lat_min: float + lat_max: float + lon_min: float + lon_max: float + name: str + + def mask(self, lat: np.ndarray, lon: np.ndarray) -> np.ndarray: + lon = np.asarray(lon, dtype=np.float64) + lat = np.asarray(lat, dtype=np.float64) + lon360 = lon.copy() + lon360[~np.isfinite(lon360)] = np.nan + lon360[lon360 < 0] += 360.0 + return ( + np.isfinite(lat) + & np.isfinite(lon360) + & (lat >= self.lat_min) + & (lat <= self.lat_max) + & (lon360 >= self.lon_min) + & (lon360 <= self.lon_max) + ) + + +def _first_region_box_from_config(config_path: Path) -> Tuple[float, float, float, float]: + with open(config_path, "r", encoding="utf-8") as f: + cfg = json.load(f) + boxes = cfg.get("data_filtering_config", {}).get("region_boxes") + if not boxes or len(boxes[0]) != 4: + raise ValueError(f"No valid data_filtering_config.region_boxes in {config_path}") + lat_min, lat_max, lon_min, lon_max = (float(x) for x in boxes[0]) + return lat_min, lat_max, lon_min, lon_max + + +def _layer_columns(df: pd.DataFrame, var: str) -> List[str]: + prefix = f"Y_{var}_col1_layer" + cols = [c for c in df.columns if c.startswith(prefix)] + cols.sort(key=lambda c: int(c.replace(prefix, ""))) + return cols + + +def _stats(x: np.ndarray) -> Dict[str, float]: + x = np.asarray(x, dtype=np.float64).ravel() + x = x[np.isfinite(x)] + if x.size == 0: + return {"n": 0, "mean": np.nan, "std": np.nan, "p05": np.nan, "p50": np.nan, "p95": np.nan, "min": np.nan, "max": np.nan} + return { + "n": int(x.size), + "mean": float(np.mean(x)), + "std": float(np.std(x)), + "p05": float(np.percentile(x, 5)), + "p50": float(np.percentile(x, 50)), + "p95": float(np.percentile(x, 95)), + "min": float(np.min(x)), + "max": float(np.max(x)), + } + + +def _optional_natveg_mask(gt: pd.DataFrame, static_path: Path) -> Optional[np.ndarray]: + if not static_path.is_file(): + return None + st = pd.read_csv(static_path) + if "Latitude" not in st.columns or "Longitude" not in st.columns: + return None + pct_nv = None + for c in ("PCT_NATVEG", "pct_natveg"): + if c in st.columns: + pct_nv = pd.to_numeric(st[c], errors="coerce").to_numpy() + break + pct0 = None + for c in ("PCT_NAT_PFT_0",): + if c in st.columns: + pct0 = pd.to_numeric(st[c], errors="coerce").to_numpy() + break + if pct_nv is None: + return None + keep = pct_nv > 0 + if pct0 is not None: + keep &= pct0 < 100.0 + # Align rows: merge on rounded lon/lat (grid snap) + gt_k = gt.assign( + _lon=np.round(gt["Longitude"].to_numpy(dtype=np.float64), 5), + _lat=np.round(gt["Latitude"].to_numpy(dtype=np.float64), 5), + ) + st_k = st.assign( + _lon=np.round(st["Longitude"].to_numpy(dtype=np.float64), 5), + _lat=np.round(st["Latitude"].to_numpy(dtype=np.float64), 5), + ) + flags = pd.DataFrame({"_lon": gt_k["_lon"], "_lat": gt_k["_lat"]}) + meta = st_k[["_lon", "_lat"]].copy() + meta["_natveg_ok"] = keep.astype(bool) + meta = meta.drop_duplicates(subset=["_lon", "_lat"], keep="first") + merged = flags.merge(meta, on=["_lon", "_lat"], how="left") + if merged["_natveg_ok"].isna().all(): + return None + filled = merged["_natveg_ok"].fillna(True).to_numpy(dtype=bool) + return filled + + +def main() -> None: + p = argparse.ArgumentParser(description="Compare 5P ground truth: Amazon vs Africa (inference CSVs).") + p.add_argument( + "--run-dir", + type=str, + required=True, + help="Run directory containing cnp_inference_entire_dataset (or set --inference-dir).", + ) + p.add_argument( + "--inference-subdir", + type=str, + default="cnp_inference_entire_dataset", + help="Subdirectory under run-dir with cnp_predictions/soil_2d_ground_truth (default: cnp_inference_entire_dataset).", + ) + p.add_argument( + "--inference-dir", + type=str, + default=None, + help="If set, use this path directly as the inference folder (overrides run-dir + inference-subdir).", + ) + p.add_argument( + "--amazon-region-config", + type=str, + default=str(REPO_ROOT / "config/training_config_amazon_5p_box.json"), + help="JSON with data_filtering_config.region_boxes[0] = Amazon box.", + ) + p.add_argument( + "--africa-region-config", + type=str, + default=str(REPO_ROOT / "config/training_config_africa_5p_box.json"), + help="JSON with data_filtering_config.region_boxes[0] = Africa box.", + ) + p.add_argument( + "--regions-json", + type=str, + default=None, + help="Optional JSON: {\"amazon\": [lat_min,lat_max,lon_min,lon_max], \"africa\": [...]} overrides config paths.", + ) + p.add_argument( + "--natveg-filter", + action="store_true", + help="Restrict to natveg cells (PCT_NATVEG>0, PCT_NAT_PFT_0<100) via merge with test_static_inverse.csv.", + ) + p.add_argument("--output", type=str, default=None, help="Write text report to this path.") + p.add_argument("--output-csv", type=str, default=None, help="Write long-form stats CSV.") + args = p.parse_args() + + run_dir = Path(args.run_dir).resolve() + if args.inference_dir: + inf_dir = Path(args.inference_dir).resolve() + else: + inf_dir = (run_dir / args.inference_subdir).resolve() + gt_dir = inf_dir / "cnp_predictions" / "soil_2d_ground_truth" + static_path = inf_dir / "cnp_predictions" / "test_static_inverse.csv" + + if args.regions_json: + with open(args.regions_json, "r", encoding="utf-8") as f: + rj = json.load(f) + am = rj["amazon"] + af = rj["africa"] + amazon = RegionBox(float(am[0]), float(am[1]), float(am[2]), float(am[3]), "amazon") + africa = RegionBox(float(af[0]), float(af[1]), float(af[2]), float(af[3]), "africa") + else: + am = _first_region_box_from_config(Path(args.amazon_region_config)) + af = _first_region_box_from_config(Path(args.africa_region_config)) + amazon = RegionBox(am[0], am[1], am[2], am[3], "amazon") + africa = RegionBox(af[0], af[1], af[2], af[3], "africa") + + lines: List[str] = [] + lines.append("5P ground truth: Amazon vs Africa (from inference soil_2d_ground_truth CSVs)") + lines.append(f"Inference dir: {inf_dir}") + lines.append(f"Amazon box [lat_min, lat_max, lon_min, lon_max]: [{amazon.lat_min}, {amazon.lat_max}, {amazon.lon_min}, {amazon.lon_max}]") + lines.append(f"Africa box: [{africa.lat_min}, {africa.lat_max}, {africa.lon_min}, {africa.lon_max}]") + lines.append(f"Natveg filter: {args.natveg_filter}") + lines.append("") + + rows_csv: List[dict] = [] + + for var in FIVE_P_VARS: + path = gt_dir / f"ground_truth_Y_{var}.csv" + if not path.is_file(): + lines.append(f"MISSING {path.name}") + continue + df = pd.read_csv(path) + if "Longitude" not in df.columns or "Latitude" not in df.columns: + print(f"ERROR: {path} missing Longitude/Latitude", file=sys.stderr) + sys.exit(1) + layer_cols = _layer_columns(df, var) + if len(layer_cols) != 10: + print(f"WARNING: {var}: expected 10 layer columns, found {len(layer_cols)}", file=sys.stderr) + + lat = df["Latitude"].to_numpy(dtype=np.float64) + lon = df["Longitude"].to_numpy(dtype=np.float64) + m_am = amazon.mask(lat, lon) + m_af = africa.mask(lat, lon) + base_mask = np.ones(len(df), dtype=bool) + if args.natveg_filter: + nv = _optional_natveg_mask(df, static_path) + if nv is None: + print("WARNING: --natveg-filter set but could not build mask from test_static_inverse; ignoring.", file=sys.stderr) + else: + base_mask = nv + m_am &= base_mask + m_af &= base_mask + + mat = df[layer_cols].to_numpy(dtype=np.float64) if layer_cols else np.zeros((len(df), 0)) + row_sum = np.nansum(mat, axis=1) if mat.size else np.zeros(len(df)) + + lines.append(f"## {var}") + lines.append(f" Cells in Amazon box (after filters): {int(np.sum(m_am))}") + lines.append(f" Cells in Africa box (after filters): {int(np.sum(m_af))}") + + for region_name, mask in (("amazon", m_am), ("africa", m_af)): + pooled = mat[mask].ravel() if mat.size else np.array([]) + st_p = _stats(pooled) + st_sum = _stats(row_sum[mask]) + lines.append(f" [{region_name}] all layer values — n={st_p['n']} mean={st_p['mean']:.6g} std={st_p['std']:.6g} p50={st_p['p50']:.6g} min={st_p['min']:.6g} max={st_p['max']:.6g}") + lines.append(f" [{region_name}] per-cell sum(10 layers) — n={st_sum['n']} mean={st_sum['mean']:.6g} p50={st_sum['p50']:.6g}") + + for kind, st in (("pooled_layers", st_p), ("row_sum_10L", st_sum)): + rows_csv.append({"variable": var, "region": region_name, "stat_kind": kind, **st}) + + # Per-layer mean contrast (Amazon mean - Africa mean) at each layer + if mat.size and np.any(m_am) and np.any(m_af): + prof_am = np.nanmean(mat[m_am], axis=0) + prof_af = np.nanmean(mat[m_af], axis=0) + lines.append(f" Per-layer mean (Amazon - Africa): {np.round(prof_am - prof_af, 6).tolist()}") + lines.append("") + + report = "\n".join(lines) + print(report) + + out_txt = Path(args.output) if args.output else run_dir / "analysis" / "5p_gt_amazon_vs_africa_report.txt" + out_txt.parent.mkdir(parents=True, exist_ok=True) + out_txt.write_text(report, encoding="utf-8") + print(f"Wrote {out_txt}") + + if args.output_csv: + out_c = Path(args.output_csv) + else: + out_c = run_dir / "analysis" / "5p_gt_amazon_vs_africa_stats.csv" + if rows_csv: + out_c.parent.mkdir(parents=True, exist_ok=True) + pd.DataFrame(rows_csv).to_csv(out_c, index=False) + print(f"Wrote {out_c}") + + +if __name__ == "__main__": + main() From c174904e7aa59a3dc9d1add0cac0806a601575e8 Mon Sep 17 00:00:00 2001 From: Dali Wang Date: Thu, 28 May 2026 14:16:51 -0400 Subject: [PATCH 3/5] Support full-domain TESNORTHERA inference and restart export Add batched inference to avoid GPU OOM on ~260K gridcells, restore export-only CSV generation from saved predictions, and switch restart spatial mapping to KD-tree nearest-neighbor lookup to avoid huge distance matrices. --- commands.txt | 34 +- docs/CNP_pipeline_runbook.md | 2 + ..._5P_TWO_REGIONS_BIAS_SCALE_AND_RESTARTS.md | 10 + .../WORKFLOW_PHASE1_PHASE2_PHASE3_RESTARTS.md | 3 + scripts/ai_predictions_to_restart.py | 11 +- scripts/apply_5p_bias_scale_correction.py | 221 +++++++- .../compare_5p_gt_two_regions_inference.py | 501 ++++++++++++++---- scripts/export_saved_predictions_to_csv.py | 191 +++++++ .../merge_5p_bias_corrected_amazon_africa.py | 10 +- scripts/run_inference_all.py | 153 ++++-- 10 files changed, 960 insertions(+), 176 deletions(-) create mode 100644 scripts/export_saved_predictions_to_csv.py diff --git a/commands.txt b/commands.txt index 3a90c8c..bed024b 100644 --- a/commands.txt +++ b/commands.txt @@ -565,4 +565,36 @@ python scripts/ai_predictions_to_netcdf.py \ --site-name amazon303_tworegion_biascorrection \ --output-dir cnp_results/run_20260311_204845_phase3_tworegions/analysis/amazon_5p_restart_comparison_twocorrections - python train_cnp_model.py --training-config-json config/training_config_phase2_tropical_soilp_only.json --variable-list CNP_IO_updated9_dev_dw.txt --tropical-lat-range -30,30 --output-dir-suffix phase2_tropical_soilp_amazon_africa \ No newline at end of file + python train_cnp_model.py --training-config-json config/training_config_phase2_tropical_soilp_only.json --variable-list CNP_IO_updated9_dev_dw.txt --tropical-lat-range -30,30 --output-dir-suffix phase2_tropical_soilp_amazon_africa + + +# --- 05/28/2026: TESNORTHERA full-domain export-only CSV + normal-spinup restart --- +# Uses the completed batched model inference predictions.pkl for the full north domain +# (259,535 samples), so this does not rerun the expensive GPU inference. +cd /mnt/proj-shared/AI4BGC_7xw/AI4BGC + +TESNORTH_RUN=/mnt/proj-shared/AI4BGC_7xw/AI4BGC/cnp_results/run_20260528_0946_tesnorthera_inference_batched_from_tesnorth10pc +TESNORTH_INFER="$TESNORTH_RUN/cnp_inference_entire_domain" +TESNORTH_TRAIN_RUN=/mnt/proj-shared/AI4BGC_7xw/AI4BGC/cnp_results/run_20260319_134216_tesnorth10pc +TESNORTH_VAR_LIST=/mnt/proj-shared/AI4BGC_7xw/AI4BGC/CNP_IO_tesnorth10pct.txt +TESNORTH_RESTART=/mnt/proj-shared/AI4BGC_7xw/AI4BGC/ELM_data/TESNORTHERA5_20adspin_restartfile/uELM_NORTHERA5_ERA5REF_I1850uELMCNPRDCTCBC.elm.r.0021-01-01-00000.nc + +python scripts/export_saved_predictions_to_csv.py \ + --inference-dir "$TESNORTH_INFER" \ + --config "$TESNORTH_TRAIN_RUN/cnp_config.json" \ + --scalers-dir "$TESNORTH_TRAIN_RUN/cnp_predictions/scalers" \ + --derive-np-from-c \ + 2>&1 | tee "$TESNORTH_RUN/export_saved_predictions_to_csv.log" + +python scripts/ai_predictions_to_netcdf.py \ + --ai-predictions "$TESNORTH_INFER/cnp_predictions" \ + --variable-list "$TESNORTH_VAR_LIST" \ + --output "$TESNORTH_INFER/comparison_results/ai_predictions_for_restart.nc" \ + 2>&1 | tee "$TESNORTH_RUN/ai_predictions_to_netcdf.log" + +python -u scripts/ai_predictions_to_restart.py \ + --ai-predictions "$TESNORTH_INFER/comparison_results/ai_predictions_for_restart.nc" \ + --restart-file "$TESNORTH_RESTART" \ + --output "$TESNORTH_RUN/updated_restart_normal_spinup_uELM_NORTHERA5_ERA5REF_I1850uELMCNPRDCTCBC.elm.r.0021-01-01-00000.nc" \ + --variable-list "$TESNORTH_VAR_LIST" \ + 2>&1 | tee "$TESNORTH_RUN/ai_predictions_to_restart.kdtree.log" \ No newline at end of file diff --git a/docs/CNP_pipeline_runbook.md b/docs/CNP_pipeline_runbook.md index 5a58855..1ffa81b 100644 --- a/docs/CNP_pipeline_runbook.md +++ b/docs/CNP_pipeline_runbook.md @@ -161,6 +161,8 @@ python ../../scripts/generate_prediction_quality_report.py cnp_inference_entire_ - Second reads that CSV and writes `cnp_inference_entire_dataset/analysis/` (quality report, plots, `quality_summary_report.txt`, `top_bad_plots/`, etc.). - Replace `cnp_inference_entire_dataset` with your actual inference subdirectory name if different. +**Amazon / Africa five-P only:** To compare predictions to ground truth inside the standard Amazon and Africa boxes (and across multiple runs), use `scripts/compare_5p_gt_two_regions_inference.py` (`pred-eval`). Methodology for Phase3 bias correction of those variables: [REPORT_PHASE3_5P_BIAS_CORRECTION_REVIEW.md](./REPORT_PHASE3_5P_BIAS_CORRECTION_REVIEW.md); workflow: [WORKFLOW_5P_TWO_REGIONS_BIAS_SCALE_AND_RESTARTS.md](./WORKFLOW_5P_TWO_REGIONS_BIAS_SCALE_AND_RESTARTS.md). + ### 4.1) Highlight a site in all scatter plots (optional) Highlight a given site (lon/lat) in every GT vs Pred scatter plot. The site will be marked as a red star. diff --git a/docs/WORKFLOW_5P_TWO_REGIONS_BIAS_SCALE_AND_RESTARTS.md b/docs/WORKFLOW_5P_TWO_REGIONS_BIAS_SCALE_AND_RESTARTS.md index 48580f1..65fc591 100644 --- a/docs/WORKFLOW_5P_TWO_REGIONS_BIAS_SCALE_AND_RESTARTS.md +++ b/docs/WORKFLOW_5P_TWO_REGIONS_BIAS_SCALE_AND_RESTARTS.md @@ -12,6 +12,13 @@ This document records the end-to-end procedure to: The goal is to have a reproducible protocol you can follow again for future runs. +**See also** + +- [REPORT_PHASE3_5P_BIAS_CORRECTION_REVIEW.md](./REPORT_PHASE3_5P_BIAS_CORRECTION_REVIEW.md) — methodology review, reference-site weighting vs full-box fitting, and recommendations for `solutionp_vr` / `occlp_vr`. +- [PHASE3_TWOREGIONS_VS_GOOD_RUN_DIFFERENCE.md](./PHASE3_TWOREGIONS_VS_GOOD_RUN_DIFFERENCE.md) — why Phase3 uses Amazon-only → Africa-only → `merge_5p_bias_corrected_amazon_africa.py` like the “good” Phase2 run. +- [WORKFLOW_PHASE1_PHASE2_PHASE3_RESTARTS.md](./WORKFLOW_PHASE1_PHASE2_PHASE3_RESTARTS.md) — three-phase restart orchestration (`run_phase3_tworegions.sh`, merged 5P CSV layout). +- Regional **prediction vs ground truth** tables: `scripts/compare_5p_gt_two_regions_inference.py` (`pred-eval`); broader pipeline: [CNP_pipeline_runbook.md](./CNP_pipeline_runbook.md). + --- ### 1. Train (or select) the global model: `natveg_improved` @@ -171,7 +178,10 @@ This is the recommended path if you want a **single global NetCDF** where: predictions as the baseline. The fit and corrections are restricted to the two-region boxes (Amazon + Africa); elsewhere the CSVs remain raw Phase2. +For **Phase3-style** workflows, the recommended pattern is **two** passes of `apply_5p_bias_scale_correction.py` (Amazon-only config, then Africa-only config) followed by **`scripts/merge_5p_bias_corrected_amazon_africa.py`**. Rationale and tuning notes: [REPORT_PHASE3_5P_BIAS_CORRECTION_REVIEW.md](./REPORT_PHASE3_5P_BIAS_CORRECTION_REVIEW.md), [PHASE3_TWOREGIONS_VS_GOOD_RUN_DIFFERENCE.md](./PHASE3_TWOREGIONS_VS_GOOD_RUN_DIFFERENCE.md). + - **Script:** `scripts/apply_5p_bias_scale_correction.py` +- **Merge (split regional fits into one CSV set):** `scripts/merge_5p_bias_corrected_amazon_africa.py` - **Key command (example):** ```bash diff --git a/docs/WORKFLOW_PHASE1_PHASE2_PHASE3_RESTARTS.md b/docs/WORKFLOW_PHASE1_PHASE2_PHASE3_RESTARTS.md index 4d52591..92ad0a9 100644 --- a/docs/WORKFLOW_PHASE1_PHASE2_PHASE3_RESTARTS.md +++ b/docs/WORKFLOW_PHASE1_PHASE2_PHASE3_RESTARTS.md @@ -44,6 +44,9 @@ It is based on: - `docs/WORKFLOW_5P_TWO_REGIONS_BIAS_SCALE_AND_RESTARTS.md` - `docs/INSTRUCTIONS_TRENDY_1_AI_RESTART_CREATION.md` +- [REPORT_PHASE3_5P_BIAS_CORRECTION_REVIEW.md](./REPORT_PHASE3_5P_BIAS_CORRECTION_REVIEW.md) — bias/scale methodology, full-region objectives for `solutionp_vr` / `occlp_vr`, and improvement directions. +- [PHASE3_TWOREGIONS_VS_GOOD_RUN_DIFFERENCE.md](./PHASE3_TWOREGIONS_VS_GOOD_RUN_DIFFERENCE.md) — split Amazon/Africa correction + merge vs older single-config Phase3. +- Quantitative **Amazon / Africa 5P pred vs GT**: `scripts/compare_5p_gt_two_regions_inference.py` (`pred-eval`). The recommended way to run all three phases is to use `scripts/run_phase_restarts.sh`, described at the end of this document. diff --git a/scripts/ai_predictions_to_restart.py b/scripts/ai_predictions_to_restart.py index 5ec65af..e3a2de3 100644 --- a/scripts/ai_predictions_to_restart.py +++ b/scripts/ai_predictions_to_restart.py @@ -87,13 +87,16 @@ def create_spatial_mapping(ds_ai: xr.Dataset, ds_model: xr.Dataset) -> tuple[np. print(f" AI coordinates: {len(ai_lon)} gridcells") print(f" Model coordinates: {len(model_lon)} gridcells") - # Create spatial mapping using nearest neighbor (exact same as working script) - from scipy.spatial.distance import cdist + # Create spatial mapping using nearest neighbor without materializing the + # full model x AI distance matrix. The full northern domain has ~260k + # gridcells, so scipy.spatial.distance.cdist would require hundreds of GB. + from scipy.spatial import cKDTree ai_coords = np.column_stack([ai_lon, ai_lat]) model_coords = np.column_stack([model_lon, model_lat]) - distances = cdist(model_coords, ai_coords) - model_to_ai_mapping = np.argmin(distances, axis=1) # 索引是模型格点, 值是最近的 AI 格点 + tree = cKDTree(ai_coords) + distances, model_to_ai_mapping = tree.query(model_coords, k=1) print(f" Spatial mapping created: {len(model_to_ai_mapping)} Model -> {len(set(model_to_ai_mapping))} AI") + print(f" Nearest-neighbor distance max: {float(np.max(distances)):.6g}") # Get grid information from the MODEL file as the master coordinate system n_grid = ds_model.sizes["gridcell"] print(f" Using MODEL gridcell count: {n_grid}") diff --git a/scripts/apply_5p_bias_scale_correction.py b/scripts/apply_5p_bias_scale_correction.py index 816c654..372344e 100644 --- a/scripts/apply_5p_bias_scale_correction.py +++ b/scripts/apply_5p_bias_scale_correction.py @@ -49,7 +49,9 @@ REFERENCE_SITES: List[Tuple[float, float]] = AMAZON_REFERENCE_SITES + AFRICA_REFERENCE_SITES REFERENCE_SITE_EXTRA_WEIGHT_AMAZON: float = 500.0 # extra weight for Amazon reference cells -REFERENCE_SITE_EXTRA_WEIGHT_AFRICA: float = 250.0 # slightly softer extra weight for Africa +REFERENCE_SITE_EXTRA_WEIGHT_AFRICA: float = 250.0 # slightly softer extra weight for Africa (legacy / v2) +# regional_fit_v3: match Amazon ref-site emphasis in Africa for solutionp_vr / occlp_vr only +REFERENCE_SITE_EXTRA_WEIGHT_AFRICA_V3: float = REFERENCE_SITE_EXTRA_WEIGHT_AMAZON REFERENCE_SITE_EXTRA_WEIGHT: float = REFERENCE_SITE_EXTRA_WEIGHT_AMAZON REFERENCE_SITE_ATOL: float = 1e-4 @@ -68,6 +70,20 @@ def contains(self, lat: np.ndarray, lon: np.ndarray) -> np.ndarray: return (lat >= self.lat_min) & (lat <= self.lat_max) & (lon >= self.lon_min) & (lon <= self.lon_max) +def _region_geo_kind(region: RegionBox) -> str: + """Classify box as Amazon vs Africa from lon/lat bounds. + + ``load_region_boxes`` labels the first JSON box ``amazon`` even for Africa-only configs; + geographic bounds disambiguate for ref-site weighting and blend behavior. + """ + lon_mid = 0.5 * (float(region.lon_min) + float(region.lon_max)) + if float(region.lon_min) >= 200.0 or lon_mid >= 200.0: + return "amazon" + if float(region.lon_max) <= 40.0 and float(region.lon_min) >= 0.0: + return "africa" + return "other" + + def load_region_boxes( config_path: str, split_first_region_lat: Optional[float] = None, @@ -135,6 +151,7 @@ def fit_bias_scale( ref_sites: Optional[List[Tuple[float, float]]] = None, ref_site_extra_weight: float = 500.0, ref_site_atol: float = 1e-4, + relative_weight_cap: Optional[float] = None, ) -> Tuple[float, float]: """Fit GT ≈ a * pred + b via least squares (or a * pred only if multiplicative_only). @@ -175,9 +192,12 @@ def fit_bias_scale( if relative_error_weighted: # Minimize sum w_i * (y_i - a*x_i - b)^2 with w_i = 1/(y_i+eps)^2 w = 1.0 / (y.astype(float) + relative_eps) ** 2 + if relative_weight_cap is not None and relative_weight_cap > 0: + w = np.minimum(w, float(relative_weight_cap)) # Boost weight for reference-site cells so fit targets <10% error there if ( - ref_sites + ref_site_extra_weight > 0 + and ref_sites and ref_site_lon is not None and ref_site_lat is not None and ref_site_lon.shape == ref_site_lat.shape @@ -201,6 +221,70 @@ def fit_bias_scale( return float(a), float(b) +def _ref_site_project_minimal_change( + a: float, + b: float, + pred_ref: float, + gt_ref: float, + scale: Optional[float], +) -> Tuple[float, float]: + """One-step projection so a'*pred_ref + b' matches gt_ref (scaled-space if scale set).""" + if scale is not None and scale != 1.0: + pred_s = pred_ref * scale + tgt = gt_ref * scale + lam = 2.0 * (a * pred_s + b - tgt) / (pred_s**2 + 1.0) + a = a - lam * pred_s / 2.0 + b = b - lam / 2.0 + else: + lam = 2.0 * (a * pred_ref + b - gt_ref) / (pred_ref**2 + 1.0) + a = float(a - lam * pred_ref / 2.0) + b = float(b - lam / 2.0) + return a, b + + +def _iterative_ref_site_projection( + a: float, + b: float, + pred_vals: np.ndarray, + gt_vals: np.ndarray, + ref_lon: np.ndarray, + ref_lat: np.ndarray, + region_ref_sites: List[Tuple[float, float]], + scale: Optional[float], + max_sweeps: int = 8, + rel_tol_pct: float = 10.0, +) -> Tuple[float, float]: + """Repeatedly nudge (a,b) so listed reference cells approach GT; multiple sweeps for Africa-scale misfit.""" + for _ in range(max_sweeps): + improved = False + for rlon, rlat in region_ref_sites: + at_site = np.isclose(ref_lon, rlon, atol=REFERENCE_SITE_ATOL) & np.isclose( + ref_lat, rlat, atol=REFERENCE_SITE_ATOL + ) + if not np.any(at_site): + continue + idx = int(np.where(at_site)[0][0]) + pred_ref = float(pred_vals[idx]) + gt_ref = float(gt_vals[idx]) + if scale and scale != 1.0: + pred_s = pred_ref * scale + corrected_ref = (a * pred_s + b) / scale + else: + corrected_ref = a * pred_ref + b + if corrected_ref < 0: + a, b = _ref_site_project_minimal_change(a, b, pred_ref, gt_ref, scale) + improved = True + continue + rel_pct = 100.0 * abs(corrected_ref - gt_ref) / (gt_ref + 1e-12) + if rel_pct <= rel_tol_pct: + continue + a, b = _ref_site_project_minimal_change(a, b, pred_ref, gt_ref, scale) + improved = True + if not improved: + break + return a, b + + def compute_and_apply_corrections( run_dir: str, region_boxes: List[RegionBox], @@ -210,6 +294,9 @@ def compute_and_apply_corrections( scale_factor_vars: Optional[Dict[str, float]] = None, relative_error_weighted_vars: Optional[Set[str]] = None, reference_sites: Optional[List[Tuple[float, float]]] = None, + *, + regional_fit_v2: bool = False, + regional_fit_v3: bool = False, ) -> Dict[str, Dict[str, Dict[str, Dict[str, float]]]]: """Compute per-variable, per-layer, per-region bias/scale corrections and apply them. @@ -224,6 +311,43 @@ def compute_and_apply_corrections( relative_error_weighted_vars = RELATIVE_ERROR_WEIGHTED_VARS if reference_sites is None: reference_sites = REFERENCE_SITES + + # regional_fit_v2: region-wide emphasis for solutionp_vr / occlp_vr (see docs/REPORT_PHASE3_5P_BIAS_CORRECTION_REVIEW.md) + # regional_fit_v3: same Amazon treatment as v2; Africa gets Amazon-strength ref-site WLS + iterative ref projection. + if regional_fit_v2 and regional_fit_v3: + raise ValueError("regional_fit_v2 and regional_fit_v3 are mutually exclusive") + + if regional_fit_v3: + relative_weight_cap_mode: Optional[float] = 1e12 + symmetric_reference_projection = True + africa_correction_alpha = 1.0 + africa_rel_ref_weight = REFERENCE_SITE_EXTRA_WEIGHT_AFRICA_V3 + africa_iterative_ref_projection = True + elif regional_fit_v2: + relative_weight_cap_mode = 1e12 + symmetric_reference_projection = True + africa_correction_alpha = 1.0 + africa_rel_ref_weight = REFERENCE_SITE_EXTRA_WEIGHT_AFRICA + africa_iterative_ref_projection = False + else: + relative_weight_cap_mode = None + symmetric_reference_projection = False + africa_correction_alpha = 0.7 + africa_rel_ref_weight = REFERENCE_SITE_EXTRA_WEIGHT_AFRICA + africa_iterative_ref_projection = False + + if regional_fit_v2: + print( + "regional_fit_v2: cap relative-error weights, disable reference-site WLS boost, " + "symmetric ref-site projection (Amazon + Africa), full Africa correction (no blend)." + ) + if regional_fit_v3: + print( + "regional_fit_v3: Amazon same as v2 (capped WLS, no ref inflation). " + "Africa: capped WLS + Amazon-strength ref-site weights for solutionp_vr/occlp_vr, " + "iterative ref-site projection, full correction (no blend)." + ) + predictions_root = os.path.join( run_dir, "cnp_inference_entire_dataset", @@ -370,10 +494,11 @@ def compute_and_apply_corrections( # Choose region-specific reference sites and extra weights if rel_err: - if region.name.startswith("amazon"): + gk = _region_geo_kind(region) + if gk == "amazon": region_ref_sites = AMAZON_REFERENCE_SITES region_extra_weight = REFERENCE_SITE_EXTRA_WEIGHT_AMAZON - elif region.name.startswith("africa"): + elif gk == "africa": region_ref_sites = AFRICA_REFERENCE_SITES region_extra_weight = REFERENCE_SITE_EXTRA_WEIGHT_AFRICA else: @@ -383,6 +508,16 @@ def compute_and_apply_corrections( region_ref_sites = None region_extra_weight = REFERENCE_SITE_EXTRA_WEIGHT + if rel_err and region_ref_sites: + if regional_fit_v2: + w_extra = 0.0 + elif regional_fit_v3: + w_extra = 0.0 if _region_geo_kind(region) == "amazon" else float(africa_rel_ref_weight) + else: + w_extra = region_extra_weight + else: + w_extra = region_extra_weight + ref_lon = lon[region_mask] if rel_err and region_ref_sites else None ref_lat = lat[region_mask] if rel_err and region_ref_sites else None a, b = fit_bias_scale( @@ -395,19 +530,38 @@ def compute_and_apply_corrections( ref_site_lon=ref_lon, ref_site_lat=ref_lat, ref_sites=region_ref_sites, - ref_site_extra_weight=region_extra_weight, + ref_site_extra_weight=w_extra, ref_site_atol=REFERENCE_SITE_ATOL, + relative_weight_cap=relative_weight_cap_mode, ) - # For relative-error vars: if ref site still has >10% error, nudge (a,b) so ref site is exact (minimal change). - # We only enforce this hard constraint for Amazon regions to avoid overfitting Africa to noisy GT. - if ( + # For relative-error vars: optional minimal-change projection so ref gridcells approach GT. + # Legacy: Amazon only. v2: Amazon + Africa, single nudge. v3: Africa uses iterative sweeps. + gk = _region_geo_kind(region) + apply_ref_projection = ( rel_err and region_ref_sites - and region.name.startswith("amazon") and ref_lon is not None and ref_lat is not None + and (gk == "amazon" or (symmetric_reference_projection and gk == "africa")) + ) + if ( + apply_ref_projection + and regional_fit_v3 + and gk == "africa" + and africa_iterative_ref_projection ): - for rlon, rlat in reference_sites: + a, b = _iterative_ref_site_projection( + a, + b, + pred_vals, + gt_vals, + ref_lon, + ref_lat, + region_ref_sites, + scale, + ) + elif apply_ref_projection: + for rlon, rlat in region_ref_sites: at_site = np.isclose(ref_lon, rlon, atol=REFERENCE_SITE_ATOL) & np.isclose( ref_lat, rlat, atol=REFERENCE_SITE_ATOL ) @@ -424,20 +578,15 @@ def compute_and_apply_corrections( rel_pct = 100 * abs(corrected_ref - gt_ref) / (gt_ref + 1e-12) if corrected_ref >= 0 and rel_pct <= 10: break - # Project (a,b) onto constraint a'*pred_ref + b' = gt_ref (in original space) - # Minimize (a'-a)^2+(b'-b)^2 s.t. a'*pred_ref + b' = gt_ref. Solution: - # lambda = 2*(a*pred_ref + b - gt_ref) / (pred_ref**2 + 1), a' = a - lambda*pred_ref/2, b' = b - lambda/2 if scale and scale != 1.0: pred_orig = pred_ref - # In scaled space: a*pred_s + b = gt_s would give (a*pred_s+b)/scale = gt_orig so a*pred_s+b = gt_orig*scale - # Constraint in scaled space: a'*pred_s + b' = gt_ref*scale pred_s = pred_orig * scale tgt = gt_ref * scale - lam = 2.0 * (a * pred_s + b - tgt) / (pred_s ** 2 + 1.0) + lam = 2.0 * (a * pred_s + b - tgt) / (pred_s**2 + 1.0) a = a - lam * pred_s / 2.0 b = b - lam / 2.0 else: - lam = 2.0 * (a * pred_ref + b - gt_ref) / (pred_ref ** 2 + 1.0) + lam = 2.0 * (a * pred_ref + b - gt_ref) / (pred_ref**2 + 1.0) a = float(a - lam * pred_ref / 2.0) b = float(b - lam / 2.0) break @@ -477,10 +626,9 @@ def compute_and_apply_corrections( else: raw = a * pred_vals + b - # For Africa regions, apply a modest blend between raw prediction and fully corrected value - # to avoid over-correcting in noisier GT regimes. - if region.name.startswith("africa"): - alpha = 0.7 # 70% corrected, 30% original + # For Africa regions, optionally blend toward raw pred (legacy 0.7). regional_fit_v2 uses 1.0 (full correction). + if _region_geo_kind(region) == "africa" and africa_correction_alpha < 1.0: + alpha = float(africa_correction_alpha) raw = alpha * raw + (1.0 - alpha) * pred_vals # Concentrations must be non-negative corrected_df.loc[region_mask, col_name] = np.maximum(raw, 0.0) @@ -571,6 +719,25 @@ def main() -> None: "so the default Amazon site (lat ~ -17) uses a southern fit. Float, or 'none' to disable. Default: -5.0" ), ) + fit_group = parser.add_mutually_exclusive_group() + fit_group.add_argument( + "--regional-fit-v2", + action="store_true", + help=( + "Region-wide emphasis for relative-error vars (solutionp_vr, occlp_vr): cap WLS weights, " + "disable reference-site weight inflation, apply ref-site projection in Amazon and Africa, " + "and use full corrected values in Africa (no 70/30 blend). See docs/REPORT_PHASE3_5P_BIAS_CORRECTION_REVIEW.md." + ), + ) + fit_group.add_argument( + "--regional-fit-v3", + action="store_true", + help=( + "Amazon identical to v2 (capped WLS, no ref inflation). Africa: same caps plus " + "Amazon-strength reference-site WLS weights and iterative ref-site projection for " + "solutionp_vr/occlp_vr; full Africa correction (no blend)." + ), + ) args = parser.parse_args() @@ -630,12 +797,22 @@ def main() -> None: output_subdir=args.output_subdir, multiplicative_only_vars=multiplicative_only_vars, scale_factor_vars=scale_factor_vars, + regional_fit_v2=bool(args.regional_fit_v2), + regional_fit_v3=bool(args.regional_fit_v3), ) # Save parameters to analysis/ for documentation and reuse analysis_dir = os.path.join(run_dir, "analysis") os.makedirs(analysis_dir, exist_ok=True) - params_path = os.path.join(analysis_dir, "bias_scale_params_5P_two_regions.json") + if args.regional_fit_v2 or args.regional_fit_v3: + safe_sub = "".join( + c if c.isalnum() or c in "._-" else "_" for c in args.output_subdir + ) + tag = "regional_fit_v3" if args.regional_fit_v3 else "regional_fit_v2" + params_name = f"bias_scale_params_5P_{tag}_{safe_sub}.json" + else: + params_name = "bias_scale_params_5P_two_regions.json" + params_path = os.path.join(analysis_dir, params_name) with open(params_path, "w", encoding="utf-8") as f: json.dump(params, f, indent=2) print(f"Saved bias/scale parameters to {params_path}") diff --git a/scripts/compare_5p_gt_two_regions_inference.py b/scripts/compare_5p_gt_two_regions_inference.py index 85dcf77..461afca 100644 --- a/scripts/compare_5p_gt_two_regions_inference.py +++ b/scripts/compare_5p_gt_two_regions_inference.py @@ -1,30 +1,51 @@ #!/usr/bin/env python3 """ -Compare ground truth of the five soil P variables between the Amazon and Africa -boxes used in the Phase3 / 5P bias workflow. +Two-region 5P soil tools (Amazon / Africa boxes used in the 5P bias workflow). -Region bounds are read from the same JSON configs as apply_5p_bias_scale_correction -(default: config/training_config_amazon_5p_box.json and -config/training_config_africa_5p_box.json), i.e. the same definitions as -scripts/merge_5p_bias_corrected_amazon_africa.py: +**CLI change:** you must pass a subcommand: ``gt-summary`` or ``pred-eval``. - Amazon: lat [-30, 10], lon [270, 330] (degrees, lon 0–360) - Africa: lat [-15, 15], lon [0, 30] +Subcommands +----------- + +1) **gt-summary** — Ground-truth only: descriptive stats of GT in each region (no pred vs GT). + +2) **pred-eval** — For each inference run (phase1, phase2, phase3, …), compare predictions + to ground truth **at each gridcell** in the **Amazon** box and, separately, in the + **Africa** box. Metrics are **within-region only** (no Amazon-vs-Africa comparison). + The summary CSV includes per-cell aggregates (median/mean/p90 of RMSE, MAE, bias on + column sums, mean relative error per cell) plus **pooled** metrics over all cells×layers + in the region: `pooled_rmse`, `pooled_mae`, `pooled_bias_mean`, `pooled_rel_mae`, + `pooled_r2`, `pooled_pearson_r`, `pooled_n`. + +Region bounds (default, same as ``merge_5p_bias_corrected_amazon_africa.py``): -Ground truth CSVs are read from: + Amazon: lat [-30, 10], lon [270, 330] (lon 0–360 °) + Africa: lat [-15, 15], lon [0, 30] - //cnp_predictions/soil_2d_ground_truth/ground_truth_Y_.csv +Paths (per inference folder, e.g. ``…/cnp_inference_entire_dataset``): -Each file must have Longitude, Latitude and layer columns Y__col1_layer{1..10}. + Ground truth: ``cnp_predictions/soil_2d_ground_truth/ground_truth_Y_.csv`` + Predictions: ``cnp_predictions//predictions_Y_.csv`` + (for merged 5P bias output, ``suffix`` is ``_bias_corrected`` — see ``--prediction-filename-suffix`` or ``runs-json``). -Usage (from repo root): +Examples +-------- - python scripts/compare_5p_gt_two_regions_inference.py \\ + python scripts/compare_5p_gt_two_regions_inference.py gt-summary \\ --run-dir cnp_results/run_20260315_202304_phase3_tworegions - python scripts/compare_5p_gt_two_regions_inference.py \\ - --run-dir cnp_results/run_20260315_113900_phase1_global \\ - --output-csv cnp_results/run_20260315_113900_phase1_global/analysis/5p_gt_amazon_vs_africa.csv + python scripts/compare_5p_gt_two_regions_inference.py pred-eval \\ + --inference-run phase1=/path/.../run_.../cnp_inference_entire_dataset \\ + --inference-run phase2=/path/.../run_.../cnp_inference_entire_dataset \\ + --inference-run phase3_raw=/path/.../run_.../cnp_inference_entire_dataset \\ + --output-summary analysis/5p_pred_eval_summary.csv \\ + --output-per-cell-long analysis/5p_pred_eval_per_cell_long.csv \\ + --output-per-cell-wide analysis/5p_pred_eval_per_cell_wide_rmse.csv + + python scripts/compare_5p_gt_two_regions_inference.py pred-eval \\ + --inference-run phase3_bc=/path/.../cnp_inference_entire_dataset:soil_2d_predictions_5P_bias_corrected_phase2 \\ + --prediction-filename-suffix _bias_corrected \\ + --output-summary analysis/5p_phase3_bias_corrected_summary.csv """ from __future__ import annotations @@ -129,7 +150,6 @@ def _optional_natveg_mask(gt: pd.DataFrame, static_path: Path) -> Optional[np.nd keep = pct_nv > 0 if pct0 is not None: keep &= pct0 < 100.0 - # Align rows: merge on rounded lon/lat (grid snap) gt_k = gt.assign( _lon=np.round(gt["Longitude"].to_numpy(dtype=np.float64), 5), _lat=np.round(gt["Latitude"].to_numpy(dtype=np.float64), 5), @@ -149,53 +169,177 @@ def _optional_natveg_mask(gt: pd.DataFrame, static_path: Path) -> Optional[np.nd return filled -def main() -> None: - p = argparse.ArgumentParser(description="Compare 5P ground truth: Amazon vs Africa (inference CSVs).") - p.add_argument( - "--run-dir", - type=str, - required=True, - help="Run directory containing cnp_inference_entire_dataset (or set --inference-dir).", - ) - p.add_argument( - "--inference-subdir", - type=str, - default="cnp_inference_entire_dataset", - help="Subdirectory under run-dir with cnp_predictions/soil_2d_ground_truth (default: cnp_inference_entire_dataset).", - ) - p.add_argument( - "--inference-dir", - type=str, - default=None, - help="If set, use this path directly as the inference folder (overrides run-dir + inference-subdir).", - ) - p.add_argument( - "--amazon-region-config", - type=str, - default=str(REPO_ROOT / "config/training_config_amazon_5p_box.json"), - help="JSON with data_filtering_config.region_boxes[0] = Amazon box.", - ) - p.add_argument( - "--africa-region-config", - type=str, - default=str(REPO_ROOT / "config/training_config_africa_5p_box.json"), - help="JSON with data_filtering_config.region_boxes[0] = Africa box.", - ) - p.add_argument( - "--regions-json", - type=str, - default=None, - help="Optional JSON: {\"amazon\": [lat_min,lat_max,lon_min,lon_max], \"africa\": [...]} overrides config paths.", +def _load_regions( + regions_json: Optional[str], + amazon_config: str, + africa_config: str, +) -> Tuple[RegionBox, RegionBox]: + if regions_json: + with open(regions_json, "r", encoding="utf-8") as f: + rj = json.load(f) + am = rj["amazon"] + af = rj["africa"] + return ( + RegionBox(float(am[0]), float(am[1]), float(am[2]), float(am[3]), "amazon"), + RegionBox(float(af[0]), float(af[1]), float(af[2]), float(af[3]), "africa"), + ) + am = _first_region_box_from_config(Path(amazon_config)) + af = _first_region_box_from_config(Path(africa_config)) + return RegionBox(am[0], am[1], am[2], am[3], "amazon"), RegionBox(af[0], af[1], af[2], af[3], "africa") + + +def _parse_inference_runs( + runs_json: Optional[str], + inference_run_args: List[str], + default_pred_suffix: str, +) -> List[Tuple[str, Path, str, str]]: + """Return list of (label, inference_dir, predictions_subdir, predictions_filename_suffix).""" + out: List[Tuple[str, Path, str, str]] = [] + if runs_json: + with open(runs_json, "r", encoding="utf-8") as f: + data = json.load(f) + for entry in data: + label = entry["label"] + inf = Path(entry["inference_dir"]).resolve() + sub = entry.get("predictions_subdir", "soil_2d_predictions") + # Per-run suffix in JSON; omit key for default "" (bias-corrected runs set "_bias_corrected"). + sfx = entry.get("predictions_filename_suffix", "") + out.append((label, inf, sub, sfx)) + return out + for s in inference_run_args: + if "=" not in s: + raise ValueError(f"--inference-run must be LABEL=PATH or LABEL=PATH:SUBDIR, got: {s}") + label, rest = s.split("=", 1) + label = label.strip() + rest = rest.strip() + if ":" in rest: + # Last colon separates subdir only if rest looks like path:subdir (heuristic: subdir has no /) + parts = rest.rsplit(":", 1) + if len(parts) == 2 and "/" not in parts[1] and "\\" not in parts[1]: + pth, sub = parts[0], parts[1] + else: + pth, sub = rest, "soil_2d_predictions" + else: + pth, sub = rest, "soil_2d_predictions" + out.append((label, Path(pth).resolve(), sub, default_pred_suffix)) + return out + + +def _add_cell_keys(df: pd.DataFrame, ndigits: int = 5) -> pd.DataFrame: + return df.assign( + lon_k=np.round(df["Longitude"].to_numpy(dtype=np.float64), ndigits), + lat_k=np.round(df["Latitude"].to_numpy(dtype=np.float64), ndigits), ) - p.add_argument( - "--natveg-filter", - action="store_true", - help="Restrict to natveg cells (PCT_NATVEG>0, PCT_NAT_PFT_0<100) via merge with test_static_inverse.csv.", + + +def _load_gt_pred_layers( + inf_dir: Path, + var: str, + predictions_subdir: str, + predictions_filename_suffix: str = "", +) -> Tuple[pd.DataFrame, List[str], List[str]]: + gt_path = inf_dir / "cnp_predictions" / "soil_2d_ground_truth" / f"ground_truth_Y_{var}.csv" + pr_path = ( + inf_dir / "cnp_predictions" / predictions_subdir / f"predictions_Y_{var}{predictions_filename_suffix}.csv" ) - p.add_argument("--output", type=str, default=None, help="Write text report to this path.") - p.add_argument("--output-csv", type=str, default=None, help="Write long-form stats CSV.") - args = p.parse_args() + if not gt_path.is_file(): + raise FileNotFoundError(f"Missing ground truth: {gt_path}") + if not pr_path.is_file(): + raise FileNotFoundError(f"Missing predictions: {pr_path}") + gt_df = pd.read_csv(gt_path) + pr_df = pd.read_csv(pr_path) + gt_cols = _layer_columns(gt_df, var) + pr_cols = _layer_columns(pr_df, var) + if len(gt_cols) != len(pr_cols): + raise ValueError(f"{var}: GT has {len(gt_cols)} layer cols, pred has {len(pr_cols)}") + return gt_df, gt_cols, pr_cols + + +def _align_gt_pred( + gt_df: pd.DataFrame, + pr_df: pd.DataFrame, + gt_cols: List[str], + pr_cols: List[str], +) -> Tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]: + """Return lon, lat, gt_mat (n, L), pred_mat (n, L) aligned on rounded (lon_k, lat_k).""" + g = _add_cell_keys(gt_df[["Longitude", "Latitude"] + gt_cols].copy()) + p = _add_cell_keys(pr_df[["Longitude", "Latitude"] + pr_cols].copy()) + g_idx = g.set_index(["lon_k", "lat_k"]) + p_idx = p.set_index(["lon_k", "lat_k"]) + common = g_idx.index.intersection(p_idx.index) + if len(common) == 0: + raise RuntimeError("GT/Pred merge on (lon_k, lat_k) produced zero rows — check coordinates.") + lon = g_idx.loc[common, "Longitude"].to_numpy(dtype=np.float64) + lat = g_idx.loc[common, "Latitude"].to_numpy(dtype=np.float64) + gmat = g_idx.loc[common, gt_cols].to_numpy(dtype=np.float64) + pmat = p_idx.loc[common, pr_cols].to_numpy(dtype=np.float64) + return lon, lat, gmat, pmat + + +def _pooled_vector_metrics(gt_mat: np.ndarray, pred_mat: np.ndarray, eps: float = 1e-20) -> Dict[str, float]: + """Flatten (n, L) pairs inside region; pooled RMSE, MAE, bias, R², Pearson r (ignoring nan pairs).""" + g = np.asarray(gt_mat, dtype=np.float64).ravel() + p = np.asarray(pred_mat, dtype=np.float64).ravel() + m = np.isfinite(g) & np.isfinite(p) + g, p = g[m], p[m] + if g.size < 4: + return { + "pooled_n": float(g.size), + "pooled_rmse": np.nan, + "pooled_mae": np.nan, + "pooled_bias_mean": np.nan, + "pooled_rel_mae": np.nan, + "pooled_r2": np.nan, + "pooled_pearson_r": np.nan, + } + diff = p - g + rmse = float(np.sqrt(np.mean(diff**2))) + mae = float(np.mean(np.abs(diff))) + bias = float(np.mean(diff)) + rel_mae = float(np.clip(np.mean(np.abs(diff) / (np.abs(g) + eps)), 0.0, 1e6)) + ss_res = float(np.sum(diff**2)) + g_mean = float(np.mean(g)) + ss_tot = float(np.sum((g - g_mean) ** 2)) + pooled_r2 = float(1.0 - ss_res / ss_tot) if ss_tot > eps else np.nan + if np.std(g) > eps and np.std(p) > eps: + pooled_r = float(np.corrcoef(g, p)[0, 1]) + else: + pooled_r = np.nan + return { + "pooled_n": float(g.size), + "pooled_rmse": rmse, + "pooled_mae": mae, + "pooled_bias_mean": bias, + "pooled_rel_mae": rel_mae, + "pooled_r2": pooled_r2, + "pooled_pearson_r": pooled_r, + } + + +def _per_cell_vector_metrics(gt: np.ndarray, pred: np.ndarray, eps: float = 1e-20) -> Dict[str, np.ndarray]: + """gt, pred shape (n, L). Return dict of length-n arrays.""" + diff = pred - gt + rmse_l = np.sqrt(np.nanmean(diff**2, axis=1)) + mae_l = np.nanmean(np.abs(diff), axis=1) + sum_gt = np.nansum(gt, axis=1) + sum_pr = np.nansum(pred, axis=1) + bias_sum = sum_pr - sum_gt + gt_range = np.nanmax(gt, axis=1) - np.nanmin(gt, axis=1) + # Per-layer relative error; clip extreme ratios when |gt|≈0 so regional means stay interpretable + ratio = np.abs(diff) / (np.abs(gt) + eps) + ratio = np.clip(ratio, 0.0, 1e6) + rel_mae_l = np.nanmean(ratio, axis=1) + nrmse_sum = np.clip(np.abs(bias_sum) / (gt_range + eps), 0.0, 1e6) + return { + "rmse_10L": rmse_l, + "mae_10L": mae_l, + "bias_sum": bias_sum, + "rel_mae_mean_layer": rel_mae_l, + "nrmse_sum_vs_gt_range": nrmse_sum, + } + +def cmd_gt_summary(args: argparse.Namespace) -> None: run_dir = Path(args.run_dir).resolve() if args.inference_dir: inf_dir = Path(args.inference_dir).resolve() @@ -204,23 +348,12 @@ def main() -> None: gt_dir = inf_dir / "cnp_predictions" / "soil_2d_ground_truth" static_path = inf_dir / "cnp_predictions" / "test_static_inverse.csv" - if args.regions_json: - with open(args.regions_json, "r", encoding="utf-8") as f: - rj = json.load(f) - am = rj["amazon"] - af = rj["africa"] - amazon = RegionBox(float(am[0]), float(am[1]), float(am[2]), float(am[3]), "amazon") - africa = RegionBox(float(af[0]), float(af[1]), float(af[2]), float(af[3]), "africa") - else: - am = _first_region_box_from_config(Path(args.amazon_region_config)) - af = _first_region_box_from_config(Path(args.africa_region_config)) - amazon = RegionBox(am[0], am[1], am[2], am[3], "amazon") - africa = RegionBox(af[0], af[1], af[2], af[3], "africa") + amazon, africa = _load_regions(args.regions_json, args.amazon_region_config, args.africa_region_config) lines: List[str] = [] - lines.append("5P ground truth: Amazon vs Africa (from inference soil_2d_ground_truth CSVs)") + lines.append("5P ground truth: Amazon vs Africa (distribution summary only)") lines.append(f"Inference dir: {inf_dir}") - lines.append(f"Amazon box [lat_min, lat_max, lon_min, lon_max]: [{amazon.lat_min}, {amazon.lat_max}, {amazon.lon_min}, {amazon.lon_max}]") + lines.append(f"Amazon box: [{amazon.lat_min}, {amazon.lat_max}, {amazon.lon_min}, {amazon.lon_max}]") lines.append(f"Africa box: [{africa.lat_min}, {africa.lat_max}, {africa.lon_min}, {africa.lon_max}]") lines.append(f"Natveg filter: {args.natveg_filter}") lines.append("") @@ -237,9 +370,6 @@ def main() -> None: print(f"ERROR: {path} missing Longitude/Latitude", file=sys.stderr) sys.exit(1) layer_cols = _layer_columns(df, var) - if len(layer_cols) != 10: - print(f"WARNING: {var}: expected 10 layer columns, found {len(layer_cols)}", file=sys.stderr) - lat = df["Latitude"].to_numpy(dtype=np.float64) lon = df["Longitude"].to_numpy(dtype=np.float64) m_am = amazon.mask(lat, lon) @@ -248,7 +378,7 @@ def main() -> None: if args.natveg_filter: nv = _optional_natveg_mask(df, static_path) if nv is None: - print("WARNING: --natveg-filter set but could not build mask from test_static_inverse; ignoring.", file=sys.stderr) + print("WARNING: --natveg-filter set but could not build mask; ignoring.", file=sys.stderr) else: base_mask = nv m_am &= base_mask @@ -258,24 +388,17 @@ def main() -> None: row_sum = np.nansum(mat, axis=1) if mat.size else np.zeros(len(df)) lines.append(f"## {var}") - lines.append(f" Cells in Amazon box (after filters): {int(np.sum(m_am))}") - lines.append(f" Cells in Africa box (after filters): {int(np.sum(m_af))}") + lines.append(f" Cells in Amazon box: {int(np.sum(m_am))}") + lines.append(f" Cells in Africa box: {int(np.sum(m_af))}") for region_name, mask in (("amazon", m_am), ("africa", m_af)): pooled = mat[mask].ravel() if mat.size else np.array([]) st_p = _stats(pooled) st_sum = _stats(row_sum[mask]) - lines.append(f" [{region_name}] all layer values — n={st_p['n']} mean={st_p['mean']:.6g} std={st_p['std']:.6g} p50={st_p['p50']:.6g} min={st_p['min']:.6g} max={st_p['max']:.6g}") - lines.append(f" [{region_name}] per-cell sum(10 layers) — n={st_sum['n']} mean={st_sum['mean']:.6g} p50={st_sum['p50']:.6g}") - + lines.append(f" [{region_name}] pooled layers — n={st_p['n']} mean={st_p['mean']:.6g} p50={st_p['p50']:.6g}") + lines.append(f" [{region_name}] per-cell sum(10L) — n={st_sum['n']} mean={st_sum['mean']:.6g} p50={st_sum['p50']:.6g}") for kind, st in (("pooled_layers", st_p), ("row_sum_10L", st_sum)): rows_csv.append({"variable": var, "region": region_name, "stat_kind": kind, **st}) - - # Per-layer mean contrast (Amazon mean - Africa mean) at each layer - if mat.size and np.any(m_am) and np.any(m_af): - prof_am = np.nanmean(mat[m_am], axis=0) - prof_af = np.nanmean(mat[m_af], axis=0) - lines.append(f" Per-layer mean (Amazon - Africa): {np.round(prof_am - prof_af, 6).tolist()}") lines.append("") report = "\n".join(lines) @@ -286,15 +409,205 @@ def main() -> None: out_txt.write_text(report, encoding="utf-8") print(f"Wrote {out_txt}") - if args.output_csv: - out_c = Path(args.output_csv) - else: - out_c = run_dir / "analysis" / "5p_gt_amazon_vs_africa_stats.csv" + out_c = Path(args.output_csv) if args.output_csv else run_dir / "analysis" / "5p_gt_amazon_vs_africa_stats.csv" if rows_csv: out_c.parent.mkdir(parents=True, exist_ok=True) pd.DataFrame(rows_csv).to_csv(out_c, index=False) print(f"Wrote {out_c}") +def cmd_pred_eval(args: argparse.Namespace) -> None: + from functools import reduce + + runs = _parse_inference_runs(args.runs_json, args.inference_run, args.prediction_filename_suffix or "") + if len(runs) < 1: + print("ERROR: provide --runs-json or at least one --inference-run LABEL=DIR", file=sys.stderr) + sys.exit(1) + + amazon, africa = _load_regions(args.regions_json, args.amazon_region_config, args.africa_region_config) + + summary_rows: List[dict] = [] + long_rows: List[dict] = [] + wide_fragments_by_label: Dict[str, List[pd.DataFrame]] = {label: [] for label, _, _, _ in runs} + + print("pred-eval: per-cell prediction vs ground truth (Amazon and Africa separately; not cross-region)\n") + print(f"Runs: {[r[0] for r in runs]}") + print("Prediction CSV stem: predictions_Y_.csv (suffix empty unless --prediction-filename-suffix or runs-json).\n") + print(f"Amazon box: [{amazon.lat_min}, {amazon.lat_max}, {amazon.lon_min}, {amazon.lon_max}]") + print(f"Africa box: [{africa.lat_min}, {africa.lat_max}, {africa.lon_min}, {africa.lon_max}]\n") + + for label, inf_dir, pred_subdir, pred_sfx in runs: + static_path = inf_dir / "cnp_predictions" / "test_static_inverse.csv" + for var in FIVE_P_VARS: + gt_df, gt_cols, pr_cols = _load_gt_pred_layers(inf_dir, var, pred_subdir, pred_sfx) + pr_df = pd.read_csv( + inf_dir / "cnp_predictions" / pred_subdir / f"predictions_Y_{var}{pred_sfx}.csv" + ) + lon, lat, gmat, pmat = _align_gt_pred(gt_df, pr_df, gt_cols, pr_cols) + + base_mask = np.ones(len(lon), dtype=bool) + if args.natveg_filter: + coord_df = pd.DataFrame({"Longitude": lon, "Latitude": lat}) + nv = _optional_natveg_mask(coord_df, static_path) + if nv is not None and len(nv) == len(lon): + base_mask = nv + else: + print(f"WARNING [{label} {var}]: natveg mask not applied.", file=sys.stderr) + + mets = _per_cell_vector_metrics(gmat, pmat) + + for region_box, region_name in ((amazon, "amazon"), (africa, "africa")): + rmask = region_box.mask(lat, lon) & base_mask + if not np.any(rmask): + print(f"WARNING [{label} {var} {region_name}]: zero cells", file=sys.stderr) + continue + + row = { + "run_label": label, + "predictions_subdir": pred_subdir, + "predictions_filename_suffix": pred_sfx, + "region": region_name, + "variable": var, + "n_cells": int(np.sum(rmask)), + } + for mk, arr in mets.items(): + sub = arr[rmask] + row["median_" + mk] = float(np.nanmedian(sub)) + row["mean_" + mk] = float(np.nanmean(sub)) + row["p90_" + mk] = float(np.nanpercentile(sub, 90)) + pooled = _pooled_vector_metrics(gmat[rmask], pmat[rmask]) + row.update(pooled) + summary_rows.append(row) + + lon_r = lon[rmask] + lat_r = lat[rmask] + lk = np.round(lon_r, 5) + lak = np.round(lat_r, 5) + for i in range(len(lon_r)): + long_rows.append( + { + "run_label": label, + "predictions_subdir": pred_subdir, + "predictions_filename_suffix": pred_sfx, + "region": region_name, + "variable": var, + "Longitude": float(lon_r[i]), + "Latitude": float(lat_r[i]), + "lon_k": float(lk[i]), + "lat_k": float(lak[i]), + "rmse_10L": float(mets["rmse_10L"][rmask][i]), + "mae_10L": float(mets["mae_10L"][rmask][i]), + "bias_sum": float(mets["bias_sum"][rmask][i]), + "rel_mae_mean_layer": float(mets["rel_mae_mean_layer"][rmask][i]), + "nrmse_sum_vs_gt_range": float(mets["nrmse_sum_vs_gt_range"][rmask][i]), + } + ) + + wide_fragments_by_label[label].append( + pd.DataFrame( + { + "region": region_name, + "variable": var, + "lon_k": lk, + "lat_k": lak, + "rmse_10L": mets["rmse_10L"][rmask], + } + ) + ) + + sum_df = pd.DataFrame(summary_rows) + + if args.output_summary: + outp = Path(args.output_summary) + outp.parent.mkdir(parents=True, exist_ok=True) + sum_df.to_csv(outp, index=False) + print(f"Wrote summary ({len(sum_df)} rows): {outp}") + else: + with pd.option_context("display.max_rows", 200, "display.width", 200): + print(sum_df.to_string(index=False)) + + if args.output_per_cell_long and long_rows: + outp = Path(args.output_per_cell_long) + outp.parent.mkdir(parents=True, exist_ok=True) + pd.DataFrame(long_rows).to_csv(outp, index=False) + print(f"Wrote per-cell long ({len(long_rows)} rows): {outp}") + + if args.output_per_cell_wide: + on = ["region", "variable", "lon_k", "lat_k"] + wide_dfs: List[pd.DataFrame] = [] + for label, _, _, _ in runs: + frags = wide_fragments_by_label.get(label, []) + if not frags: + continue + w = pd.concat(frags, ignore_index=True) + w = w.rename(columns={"rmse_10L": f"rmse_10L__{label}"}) + wide_dfs.append(w) + if wide_dfs: + wide_df = reduce(lambda a, b: a.merge(b, on=on, how="outer"), wide_dfs) + outp = Path(args.output_per_cell_wide) + outp.parent.mkdir(parents=True, exist_ok=True) + wide_df.to_csv(outp, index=False) + print(f"Wrote per-cell wide rmse_10L ({len(wide_df)} rows): {outp}") + else: + print("No wide table written (no data).", file=sys.stderr) + + +def main() -> None: + parser = argparse.ArgumentParser(description="Two-region 5P tools: GT summary or pred vs GT per cell.") + sub = parser.add_subparsers(dest="command", required=True) + + p_gt = sub.add_parser("gt-summary", help="GT-only distribution stats per region (legacy).") + p_gt.add_argument("--run-dir", type=str, required=True) + p_gt.add_argument("--inference-subdir", type=str, default="cnp_inference_entire_dataset") + p_gt.add_argument("--inference-dir", type=str, default=None) + p_gt.add_argument("--amazon-region-config", type=str, default=str(REPO_ROOT / "config/training_config_amazon_5p_box.json")) + p_gt.add_argument("--africa-region-config", type=str, default=str(REPO_ROOT / "config/training_config_africa_5p_box.json")) + p_gt.add_argument("--regions-json", type=str, default=None) + p_gt.add_argument("--natveg-filter", action="store_true") + p_gt.add_argument("--output", type=str, default=None) + p_gt.add_argument("--output-csv", type=str, default=None) + p_gt.set_defaults(func=cmd_gt_summary) + + p_ev = sub.add_parser( + "pred-eval", + help="Per-cell pred vs GT in Amazon and Africa; repeat runs for phase1/2/3 comparison.", + ) + p_ev.add_argument( + "--runs-json", + type=str, + default=None, + help='JSON list of {"label","inference_dir","predictions_subdir?","predictions_filename_suffix?"}.', + ) + p_ev.add_argument( + "--prediction-filename-suffix", + type=str, + default="", + help='Appended to predictions_Y_.csv for every --inference-run (e.g. _bias_corrected). Per-run: use runs-json.', + ) + p_ev.add_argument( + "--inference-run", + action="append", + default=[], + metavar="LABEL=PATH[:SUBDIR]", + help="Repeat. PATH = cnp_inference_entire_dataset. Optional :predictions_subdir after last path segment.", + ) + p_ev.add_argument("--amazon-region-config", type=str, default=str(REPO_ROOT / "config/training_config_amazon_5p_box.json")) + p_ev.add_argument("--africa-region-config", type=str, default=str(REPO_ROOT / "config/training_config_africa_5p_box.json")) + p_ev.add_argument("--regions-json", type=str, default=None) + p_ev.add_argument("--natveg-filter", action="store_true") + p_ev.add_argument("--output-summary", type=str, default=None) + p_ev.add_argument("--output-per-cell-long", type=str, default=None, help="Long CSV: one row per cell per run per variable.") + p_ev.add_argument( + "--output-per-cell-wide", + type=str, + default=None, + help="Wide CSV: keys region,variable,lon_k,lat_k and rmse_10L__