Add option to utilize GPR model to create zero eccentricity initial orbital parameters - #12
Conversation
…rbital parameters
There was a problem hiding this comment.
Pull request overview
This PR extends SimulationSupport.EccentricityControl.InitialOrbitalParameters to optionally compute zero-eccentricity initial orbital parameters using a PN baseline plus a Gaussian Process Regression (GPR) correction, and adds a Click-based CLI wrapper for the same functionality.
Changes:
- Add
method=("PN"|"GPR")andgpr_checkpointsto select PN-only vs. PN+GPR-corrected parameter estimation. - Refactor PN computation into a dedicated helper and add GPR helpers for feature vector assembly and checkpoint-based corrections.
- Add a CLI command (
initial-orbital-parameters) with text/JSON output.
Suppressed comments (2)
src/SimulationSupport/EccentricityControl/InitialOrbitalParameters.py:429
NotImplementedErroris raised with an empty message. Since this is a user-facing code path (selected bymethod == "GPR"and nonzero eccentricity), it should explain the limitation (e.g. GPR currently supports only zero eccentricity).
raise NotImplementedError("")
src/SimulationSupport/EccentricityControl/InitialOrbitalParameters.py:161
- Spelling typo in the comment (
supperted->supported).
# Only zero eccentricity is supperted here, since it utilizes ZeroEccParamsFromPN
💡 Add a code-review agent skill or configure MCP servers for context-aware, tailored reviews. Learn more in the docs.
| """Estimate the initial orbital parameters for a BBH evolution. | ||
|
|
||
| Estimates the initial coordinate separation, D_0, orbital angular velocity, Omega_0, and | ||
| radial expansion velocity, adot_0, from a Post-Newtonian approximation, optionally | ||
| corrected by a trained Gaussian Process Regression (GPR) model. | ||
|
|
||
| Specify the target eccentricity and either '--separation', '--orbital-angular-velocity', | ||
| '--num-orbits', or '---time-to-merger'. | ||
| """ | ||
| orbital_params_specified = [ | ||
| param is not None | ||
| for param in [ | ||
| separation, | ||
| orbital_angular_velocity, | ||
| num_orbits, | ||
| time_to_merger, | ||
| ] | ||
| ] | ||
| if sum(orbital_params_specified) != 1: | ||
| raise click.UsageError( | ||
| "Specify either '--separation', '--orbital-angular-velocity'," | ||
| "'--num-orbits', or '---time-to-merger'." | ||
| ) |
| if gpr_adot_checkpoint: | ||
| gpr_checkpoints["adot"] = gpr_adot_checkpoint | ||
| if not gpr_checkpoints: | ||
| raise click.UsagError( |
| import click | ||
| import numpy as np | ||
| import rich |
There was a problem hiding this comment.
Please add these to the requirements
| from SimulationSupport.gpr import ( | ||
| load_gpr_checkpoint, | ||
| predict_with_gpr_model, | ||
| ) |
| f"GPR checkpoint expects input feature {missing_feature}," | ||
| "which is not available. Available features:" | ||
| f"{sorted(available_values.keys())}. Update the" | ||
| "'available_values' mapping in this module to match" | ||
| "the checkpoint's input_features metadata." |
| corrected = {"omega": pn_omega, "adot": pn_adot} | ||
| for quantity_name, checkpoint_path in gpr_checkpoints.items(): | ||
| if quantity_name not in corrected: | ||
| raise ValueError( | ||
| f"Unknown quantity `{quantity_name}` in `gpr_checkpoints`." | ||
| "Expected one of `D_0`, `Omega_0`, or `adot_0`." | ||
| ) |
| separation: Optional[float] = None, | ||
| orbital_angular_velocity: Optional[float] = None, | ||
| radial_expansion_velocity: Optional[float] = None, | ||
| method: str = "PN", # accpets both "PN" and "GPR" |
There was a problem hiding this comment.
You don't need this comment, it's enough to have it in the docs below
| help=( | ||
| "Dimensionless spin vector of the smaller black hole, for example," | ||
| "written as '--spin-b 0.0, 0.0, 0.1'." | ||
| ), |
nilsvu
left a comment
There was a problem hiding this comment.
Please add a test that runs the added code. Please also look at the Copilot review and consider the suggestions.
| ) -> Tuple[float, float, float]: | ||
| r"""Estimate initial orbital parameters from a Post-Newtonian approximation. | ||
| r"""Estimate initial orbital parameters from a Post-Newtonian approximation | ||
| or a Gaussian Process Regression (GPR) correction to the PN approximation. |
There was a problem hiding this comment.
Please make the first line of the docs fit on one line
| separation: Optional[float] = None, | ||
| orbital_angular_velocity: Optional[float] = None, | ||
| radial_expansion_velocity: Optional[float] = None, | ||
| method: str = "PN", # accpets both "PN" and "GPR" |
There was a problem hiding this comment.
You don't need this comment, it's enough to have it in the docs below
| import click | ||
| import numpy as np | ||
| import rich |
There was a problem hiding this comment.
Please add these to the requirements
| gpr_checkpoints=gpr_checkpoints, | ||
| ) | ||
| else: | ||
| return _initial_orbital_parameters_gpr_eccentric( |
There was a problem hiding this comment.
Delete this for now. Just have _initial_orbital_parameters_gpr and assert ecc==0 in it.
|
|
||
| # Compute the initial orbital parameters from the Post-Newtonian approximation. | ||
| # Only zero eccentricity is supperted here, since it utilizes ZeroEccParamsFromPN | ||
| assert eccentricity == 0.0, ( |
There was a problem hiding this comment.
Move this assert into the _initial_orbital_parameters_pn function
| return pn_separation, orbital_angular_velocity, radial_expansion_velocity | ||
|
|
||
|
|
||
| def _initial_orbital_parameters_gpr_eccentric( |
|
|
||
| # CLI | ||
| # The function can be imported and called from Python directly, or it can be called with the CLI. | ||
| def get_initial_orbital_parameters( |
There was a problem hiding this comment.
You don't need the extra function, move this code into the CLI endpoint function below
| if sum(orbital_params_specified) != 1: | ||
| raise click.UsageError( | ||
| "Specify either '--separation', '--orbital-angular-velocity'," | ||
| "'--num-orbits', or '---time-to-merger'." |
There was a problem hiding this comment.
A lot of this checking is already done in initial_orbital_param. Can you remove some of the redundant checking? Just build the target_params and pass to initial_orbital_param, so the CLI doesn't add much logic (it should be a very thin wrapper around the Python function).
| return separation, orbital_angular_velocity, radial_expansion_velocity | ||
|
|
||
|
|
||
| def _build_gpr_feature_vector(input_features, available_values): |
There was a problem hiding this comment.
This function is only used in 1 place, can you just put the code there?
Expanded file to allow option of utilizing either PN approximation or a pre-trained GPR model to create zero eccentricity initial orbital parameters. Original PN logic remains unchanged.
method = GPR option: first computes the standard PN baseline for D_0, Omega_0, and adot_0, then applies a correction from a trained GPR model on top of the PN baseline.
currently only works for zero eccentricity GPR; eccentric GPR option is left for a future PR.
runnable both in Python and as a CLI.