Skip to content

feat: add enhanced probabilistic forecasting visualizations - #612

Open
Sir-Sloth-The-Lazy wants to merge 29 commits into
mllam:mainfrom
Sir-Sloth-The-Lazy:feat/enhanced-visualization-proposal
Open

feat: add enhanced probabilistic forecasting visualizations#612
Sir-Sloth-The-Lazy wants to merge 29 commits into
mllam:mainfrom
Sir-Sloth-The-Lazy:feat/enhanced-visualization-proposal

Conversation

@Sir-Sloth-The-Lazy

@Sir-Sloth-The-Lazy Sir-Sloth-The-Lazy commented Apr 20, 2026

Copy link
Copy Markdown
Contributor

Describe your changes

Ports the ensemble / latent-space plotting functions from the prob_model_lam branch into the refactored ForecasterModule hierarchy and adds a suite of probabilistic forecast-verification diagnostics.

To view only this PR's additions, compare against the post-refactor branch:

👉 View the clean diff on my fork

The goal is to land the visualization layer ahead of the full Graph-EFM port so that when probabilistic models arrive they already have the calibration, reliability, and uncertainty-quality plots that ensemble forecasting requires. Every function is self-contained and accepts pre-computed tensors/arrays, so model wiring can be added incrementally.

Fixes to existing plots

  • plot_error_map_absolute — unnormalized companion heatmap with a shared colour scale so errors can be compared directly across variables
  • get_var_cmap(var_name) helper — routes anomaly/error fields to RdBu_r, uncertainty/spread fields to YlOrRd, and raw physical fields to viridis
  • Ensemble-std panel in plot_ensemble_prediction uses YlOrRd (previously plasma)

New plotting functions

Function Purpose
plot_spread_skill Scatter of ensemble spread vs RMSE with 1:1 diagonal
plot_rank_histogram Talagrand bar chart with flat reference line
plot_reliability_diagram Calibration curve against 45° diagonal
plot_crps_leadtime CRPS-vs-lead-time line plot, one line per variable
plot_spaghetti Contour-bundle plot for ensemble spatial coherence
plot_plume Per-location fan diagram with 10/25/50/75/90th percentile bands
plot_exceedance_prob Probability-of-exceedance spatial map
save_forecast_animation Animated GIF from a list of figures (Pillow-based)
plot_power_spectrum Radially-averaged 2-D FFT comparison (truth vs mean vs members)
plot_latent_pca NumPy-SVD PCA scatter of prior vs variational latent samples

Ported from prob_model_lam

  • plot_ensemble_prediction — adapted to use datastore + pcolormesh via plot_on_axis (no more imshow + GRID_LIMITS)
  • plot_latent_samples — ported largely unchanged, with numpy-safe tensor handling

Wiring into ForecasterModule

  • create_metric_log_dict now logs the absolute error map alongside the normalized one for every metric
  • When the metric name contains "crps", the lead-time line plot is logged too

Dependencies: no new dependencies. Uses only existing stack (cartopy, matplotlib, numpy, torch, xarray, Pillow).

Issue Link

  • Partially addresses Merge Graph-EFM model from prob_model_lam branch #62 — ports the prob_model_lam "Add plotting functions for ensemble forecasts" item and goes further with calibration/reliability diagnostics that were not present in either branch.
  • Partially addresses Refactor model class hierarchy #49 — extends the new class hierarchy (ForecasterModule) with visualization hooks that are ready for probabilistic/ensemble Forecaster subclasses.

Type of change

  • 🐛 Bug fix (non-breaking change that fixes an issue)
  • ✨ New feature (non-breaking change that adds functionality)
  • 💥 Breaking change (fix or feature that would cause existing functionality to not work as expected)
  • 📖 Documentation (Addition or improvements to documentation)

Checklist before requesting a review

  • My branch is up-to-date with the target branch - if not update your fork with the changes from the target branch (use pull with --rebase option if possible).
  • I have performed a self-review of my code
  • For any new/modified functions/classes I have added docstrings that clearly describe its purpose, expected inputs and returned values
  • I have placed in-line comments to clarify the intent of any hard-to-understand passages of my code
  • I have updated the README to cover introduced code changes
  • I have added tests that prove my fix is effective or that my feature works
  • I have given the PR a name that clearly describes the change, written in imperative form (context).
  • I have requested a reviewer and an assignee (assignee is responsible for merging). This applies only if you have write access to the repo, otherwise feel free to tag a maintainer to add a reviewer and assignee.

Checklist for reviewers

Each PR comes with its own improvements and flaws. The reviewer should check the following:

  • the code is readable
  • the code is well tested
  • the code is documented (including return types and parameters)
  • the code is easy to maintain

Author checklist after completed review

  • I have added a line to the CHANGELOG describing this change, in a section
    reflecting type of change (add section where missing):
    • added: when you have added new functionality
    • changed: when default behaviour of the code has been changed
    • fixes: when your contribution fixes a bug
    • maintenance: when your contribution is relates to repo maintenance, e.g. CI/CD or documentation

Checklist for assignee

  • PR is up to date with the base branch
  • the tests pass
  • (if the PR is not just maintenance/bugfix) the PR is assigned to the next milestone. If it is not, propose it for a future milestone.
  • author has added an entry to the changelog (and designated the change as added, changed, fixed or maintenance)
  • Once the PR is ready to be merged, squash commits and merge the PR.

Sir-Sloth-The-Lazy and others added 29 commits February 21, 2026 17:42
- Update test_datasets.py to use ForecasterModule instead of GraphLAM
- Update test_plotting.py to use ForecasterModule instead of GraphLAM
- Fix interior_mask_bool property shape (1,) -> (N,) for correct loss masking
- Fix all_gather_cat to handle single-device runs without incorrect dim collapse
…r hierarchy

- Replace opaque argparse.Namespace with explicit keyword arguments in
  StepPredictor, BaseGraphModel, BaseHiGraphModel, GraphLAM, HiLAM,
  and HiLAMParallel __init__ methods
- Reorder methods in step_predictor.py: forward/expand_to_batch now
  appear before clamping methods
- Update all instantiation sites (train_model.py, test_training.py,
  test_prediction_model_classes.py) to pass explicit kwargs
- HiLAM helper methods (make_same/up/down_gnns) now use self.hidden_dim
  and self.hidden_layers instead of args parameter

Addresses review comments on PR mllam#208.
- Rename border to boundary in Forecaster
- Pass Forecaster object to ForecasterModule init instead of Predictor
- Remove inline imports in ForecasterModule
- Move loss-related pred_std logic fully into ForecasterModule
- Delete obsolete test_refactored_hierarchy.py
Co-authored-by: Joel Oskarsson <joel.oskarsson@outlook.com>
- Add predicts_std property to StepPredictor, Forecaster and ARForecaster
  so ForecasterModule can query the forecaster instead of taking output_std
  as a separate constructor argument
- Remove output_std parameter from ForecasterModule; use
  self._forecaster.predicts_std throughout
- Move fallback per_var_std logic out of forecast_for_batch into each
  step method so pred_std is None before fallback, enabling direct None
  checks instead of hparam checks
- Replace len(datastore.boundary_mask) with datastore.num_grid_points in
  StepPredictor to avoid relying on boundary_mask
- Move get_state_feature_weighting and ARForecaster inline imports to
  module-level imports in forecaster_module.py and train_model.py
- Fix statement ordering in StepPredictor.__init__ so register_buffer for
  grid_static_features appears directly after building the tensor
- Replace dict+loop pattern for registering state_mean/state_std buffers
  with two direct register_buffer calls
- Remove all internal Item N checklist references from comments
- Remove TORCH_FORCE_NO_WEIGHTS_ONLY_LOAD env var hack; pass
  weights_only=False explicitly to load_from_checkpoint calls and
  weights_only=True to torch.load in test_graph_creation.py
- Add test_step_predictor_no_static_features to verify models initialise
  and run correctly when the datastore returns None for static features
- Fix graph= -> graph_name= and model.forecaster -> model._forecaster in
  tests to match current API
…r_batch

Makes the forecasting path tolerant to batch-folded execution so that
future ensemble generation can fold (S, B) into (S*B) before calling
ARForecaster, without any changes to ARForecaster or StepPredictor.

Prediction is kept folded through the existing deterministic logging and
aggregation paths so all dim assumptions in training_step, validation_step,
and test_step remain correct. Unfolding to (*leading, T, N, F) is deferred
to ensemble-specific subclasses (e.g. EnsForecasterModule).

Adds test_fold_unfold_equivalence to confirm ARForecaster's rollout is
rank-transparent under a pre-entry fold.
Remove the thin forecast_for_batch wrapper, inlining the batch
unpacking and forecaster call into each step method. Also fix
_forecaster -> forecaster attribute references in checkpoint
loading and tests to match the actual attribute name.
…sue-49

Resolved conflicts:
- Removed ar_model.py (replaced by ForecasterModule/ARForecaster hierarchy)
- Updated test_training.py: use ForecasterModule imports, dynamic device
  allocation, and wandb.init(mode="disabled")
- Updated test_plotting.py: replace ARModel.all_gather_cat with
  ForecasterModule.all_gather_cat, fix model.args.create_gif reference

Ported features from ar_model.py:
- GIF export (create_gif param + plot_examples logic) into ForecasterModule
- Added common_step to ForecasterModule; refactored training/validation/
  test steps to delegate to it, eliminating duplicated forecaster calls
Ports ensemble/latent plotting functions from prob_model_lam into the
new ForecasterModule hierarchy and adds a suite of probabilistic
forecast-verification diagnostics described in docs/visualization_proposal.md.

Fixes to existing plots:
- plot_error_map_absolute: unnormalized companion heatmap with a shared
  colour scale so errors can be compared across variables
- get_var_cmap helper: picks RdBu_r for anomaly/error fields, YlOrRd for
  uncertainty/spread, viridis for raw physical fields
- plot_ensemble_prediction std panel uses YlOrRd

New plotting functions:
- plot_spread_skill: ensemble spread vs RMSE with 1:1 diagonal
- plot_rank_histogram: Talagrand diagram with flat reference
- plot_reliability_diagram: calibration curve vs 45 degree diagonal
- plot_crps_leadtime: CRPS vs lead-time line plot
- plot_spaghetti: contour-bundle for ensemble spatial coherence
- plot_plume: per-location fan diagram with 10/25/50/75/90 pct. bands
- plot_exceedance_prob: probability-of-exceedance spatial map
- save_forecast_animation: GIF from a list of figures
- plot_power_spectrum: radially-averaged 2D FFT comparison
- plot_latent_pca: NumPy-SVD PCA scatter of prior vs VI samples

Ported from prob_model_lam:
- plot_ensemble_prediction adapted to datastore + pcolormesh
- plot_latent_samples

Wiring:
- ForecasterModule.create_metric_log_dict now logs the absolute error
  map alongside the normalized one, and adds a CRPS-vs-lead-time line
  plot whenever the metric name contains "crps".
Wrap three overlong docstring/comment lines and remove dead
thresholds_np assignment left over from reliability-diagram fix.
@sadamov
sadamov self-requested a review April 20, 2026 14:14
@sadamov sadamov added the enhancement New feature or request label Apr 20, 2026
@Sir-Sloth-The-Lazy

Copy link
Copy Markdown
Contributor Author

@joeloskarsson @sadamov looking forward for your dig on this PR whenever you have time. ;D

@AswaniSahoo

Copy link
Copy Markdown

Nice work on this ! The ForecasterModule / ARForecaster / StepPredictor split looks clean, and the ensemble vis toolbox is a solid addition. A couple of things I noticed while going through the diff:

  1. plot_spread_skill - shared axis scaling across variables

The current implementation computes a single global_max across all variables:

global_max = max(spread_np.max(), rmse_np.max())
ax.set_xlim(0, global_max * 1.05)
ax.set_ylim(0, global_max * 1.05)

If one variable dominates (e.g. geopotential height having much larger absolute errors than 2m temperature), the scatter for smaller-error variables gets compressed into a corner. Per-subplot scaling (local_max = max(s.max(), r.max())) might be more informative here, since each subplot already shows a single variable anyway. The 1:1 reference line would still work the same way.

  1. plot_ensemble_prediction - silent member cap

n_members = min(samples_np.shape[0], 6) caps the displayed members at 6, but the docstring just says samples : (S, N_grid) without mentioning this. Might be worth a brief note.

  1. CRPS lead-time integration

The auto-detection in create_metric_log_dict:

if "crps" in metric_name:
    crps_fig = vis.plot_crps_leadtime(...)

is nice forward-looking design for when CRPS gets added to the metrics module.

@joeloskarsson

Copy link
Copy Markdown
Collaborator

With more things merged, could you get this up to date with main so it is easier to discuss the diff :)

@joeloskarsson joeloskarsson self-assigned this May 20, 2026
@joeloskarsson
joeloskarsson self-requested a review May 20, 2026 14:05
Comment thread neural_lam/vis.py
return s


def get_var_cmap(var_name: str) -> str:

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is this used? Don't fully get when it would be used. If I know I am plotting an error, I could just as well hard-code the colormap as call this with a hard-coded string? It would however be nice to use different cmaps different colors for different atmospheric variables, but this does not seem to be the case now?

Comment thread neural_lam/vis.py
else:
std_vmin, std_vmax = vrange_std

n_members = min(samples_np.shape[0], 6)

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Are we sure we want to limit this to 6? I guess this will need to be configured somewhere, but could be left as responsibility of caller (through number of entries in samples).

Comment thread neural_lam/vis.py
Comment on lines +596 to +600
img_side_size = int(np.sqrt(num_mesh_nodes))
assert img_side_size**2 == num_mesh_nodes, (
"Number of mesh nodes is not a square number, "
"cannot plot latent samples as images"
)

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think we should figure out a better way to do these plots, rather than using imshow for a square. In a model like Graph-EFM, the latent vectors actually have a spatial position associated with them. A better approach here is to actually use those positions. Maybe an option would be to give this method also the coordinates associated with the N_mesh dimension?

Comment thread neural_lam/vis.py


@matplotlib.rc_context(utils.fractional_plot_bundle(1))
def plot_spread_skill(spread, rmse, var_names, title=None):

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is this intended for one forecast? Or the average spread and skill across the entire dataset? How do we show which lead time each point corresponds to?

@joeloskarsson

Copy link
Copy Markdown
Collaborator

It would here be very useful if you would put in here also one example plot from each new method, as we then make sure we reason about the same visual throughout the review of this :)

@sadamov

sadamov commented May 21, 2026

Copy link
Copy Markdown
Collaborator

While I am a big fan of verification&vizualization, and thanks for this PR @Sir-Sloth-The-Lazy, we also need to be aware how costly probabilistic plots can be. There is a real risk that validation steps take longer for plotting than actual inference. So setting sensible defaults and also moving some plots to post-training evaluation with a separate library, might be reasonable here. Just wanted to drop this note and maybe ask for some performance benchmarks once the details are settled.

@Sir-Sloth-The-Lazy

Copy link
Copy Markdown
Contributor Author

@sadamov , Thank you for the review. I would be reading more about how to make plotting efficient, duely noted :). @joeloskarsson I would be now attending to this PR and completing this ! I will need sometime to attend to all of your concerns. Sorry for the delay 😅

@joeloskarsson

joeloskarsson commented May 21, 2026

Copy link
Copy Markdown
Collaborator

This is a very good point @sadamov, we should probably not implement a bunch of ensemble-related plots just because we can. This is something that would be good to get more community input on, to see what people want/need. @Sir-Sloth-The-Lazy a suggestion would be to 1. Mark in your list which plots are mainly for validation (during training) and which are for evaluation runs from a saved checkpoint. 2. Also add a figure of each plot, as the look now. 3. Either post this here, and direct people from slack to have a look, or post this directly as a thread on slack (in #neural-lam).
That way we can potentially get some more inputs on this and align with community interests.

@Sir-Sloth-The-Lazy

Copy link
Copy Markdown
Contributor Author

@joeloskarsson I have started the work keeping your and @sadamov guidance in mind :) I would pushing the changes and updation soon.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

enhancement New feature or request

Projects

None yet

Development

Successfully merging this pull request may close these issues.

4 participants