From cfe64690045f89e12db281929423340747a078d2 Mon Sep 17 00:00:00 2001 From: dgilford Date: Tue, 13 Jan 2026 15:29:58 +0000 Subject: [PATCH 1/3] Add log decomposition API --- notebooks/sample_output_analyses.ipynb | 4 +- run_sample.py | 8 +- src/tcpyPI/__init__.py | 3 + src/tcpyPI/analyses.py | 165 +++++++++++++++++++++++++ tests/test_log_decomposition_api.py | 97 +++++++++++++++ 5 files changed, 271 insertions(+), 6 deletions(-) create mode 100644 src/tcpyPI/analyses.py create mode 100644 tests/test_log_decomposition_api.py diff --git a/notebooks/sample_output_analyses.ipynb b/notebooks/sample_output_analyses.ipynb index 1d06a1e..b4ac4f7 100644 --- a/notebooks/sample_output_analyses.ipynb +++ b/notebooks/sample_output_analyses.ipynb @@ -28,7 +28,7 @@ "import matplotlib.pyplot as plt\n", "import time\n", "# load PI calculation module\n", - "from tcpyPI import pi\n", + "from tcpyPI import pi, log_decompose_pi\n", "# imports utility functions\n", "from tcpyPI.utilities import *" ] @@ -310,7 +310,7 @@ " )\n", " \n", " result1 = xr.apply_ufunc(\n", - " decompose_pi,\n", + " log_decompose_pi,\n", " pi_ds['vmax'], sst, pi_ds['t0'],\n", " kwargs=dict(CKCD=CKCD),\n", " input_core_dims=[\n", diff --git a/run_sample.py b/run_sample.py index f3f5c71..1fce2b3 100644 --- a/run_sample.py +++ b/run_sample.py @@ -11,8 +11,8 @@ import xarray as xr # load in pyPI modules -from tcpyPI import pi -from tcpyPI.utilities import * +from tcpyPI import pi, log_decompose_pi +from tcpyPI.utilities import pi_efficiency, pi_diseq_resid # define the sample data locations datdir='./data/' @@ -94,9 +94,9 @@ def run_sample_analyses(ds,CKCD=0.9): ) result = xr.apply_ufunc( - decompose_pi, + log_decompose_pi, ds['vmax'], ds['sst']+273.15, ds['t0'], - kwargs=dict(CKCD=CKCD), + kwargs=dict(CKCD=CKCD, sst_units="K"), input_core_dims=[ [], [], [], ], diff --git a/src/tcpyPI/__init__.py b/src/tcpyPI/__init__.py index 4a5fb1b..8d038b9 100644 --- a/src/tcpyPI/__init__.py +++ b/src/tcpyPI/__init__.py @@ -10,3 +10,6 @@ from .pi import * from .constants import * from .utilities import * + +# User-facing analyses API +from .analyses import log_decompose_pi, pi_log_decomposition diff --git a/src/tcpyPI/analyses.py b/src/tcpyPI/analyses.py new file mode 100644 index 0000000..2261b7d --- /dev/null +++ b/src/tcpyPI/analyses.py @@ -0,0 +1,165 @@ +"""Analysis helpers for tcpyPI. + +These functions provide a small, user-facing API for decomposing potential intensity +into logarithmic components following Wing et al. (2015), Eq. 2. + +Notes +----- +- The decomposition implemented here matches `tcpyPI.utilities.decompose_pi`: + `lnpi` is `ln(V^2)` (i.e., `2 * ln(V)`), not `ln(V)`. +- Inputs can be scalars or NumPy arrays; invalid physical states yield `nan`. +""" + +from __future__ import annotations + +from typing import Any, Dict, Literal, Tuple + +import numpy as np + + +_SSTUnits = Literal["C", "K"] + + +def _maybe_scalar(value: np.ndarray | float) -> Any: + return float(value) if np.ndim(value) == 0 else value + + +def _sst_to_kelvin(sst: Any, sst_units: _SSTUnits) -> Any: + units = sst_units.upper() + if units == "K": + return sst + if units == "C": + return np.asarray(sst) + 273.15 + raise ValueError(f"Unsupported sst_units={sst_units!r}; expected 'C' or 'K'.") + + +def log_decompose_pi( + pi: Any, + sst: Any, + t0: Any, + CKCD: float = 0.9, + *, + sst_units: _SSTUnits = "K", +) -> Tuple[Any, Any, Any, float]: + """Log-decompose potential intensity into efficiency, disequilibrium, and Ck/Cd. + + Parameters + ---------- + pi + Potential intensity wind speed (m/s). + sst + Sea surface temperature, in units given by `sst_units`. + t0 + Outflow temperature (K). + CKCD + Ratio of exchange coefficients (Ck/Cd). + sst_units + Units for `sst`: `'K'` (default) or `'C'`. + + Returns + ------- + (lnpi, lneff, lndiseq, lnCKCD) + `lnpi` is `ln(pi^2)`. + + Notes + ----- + This matches the behavior of `tcpyPI.utilities.decompose_pi`, including edge cases: + - If efficiency <= 0: all terms are `nan` except `lnCKCD`. + - If efficiency > 0 and pi <= 0: `lneff` is returned, but `lnpi` and `lndiseq` are `nan`. + """ + + pi_arr = np.asarray(pi, dtype=float) + t0_arr = np.asarray(t0, dtype=float) + sst_k = np.asarray(_sst_to_kelvin(sst, sst_units), dtype=float) + + # Fast path: preserve the original Numba-compiled scalar behavior. + if pi_arr.ndim == 0 and t0_arr.ndim == 0 and sst_k.ndim == 0: + from .utilities import decompose_pi as _decompose_pi + + return _decompose_pi(float(pi_arr), float(sst_k), float(t0_arr), CKCD=CKCD) + + pi_arr, sst_k, t0_arr = np.broadcast_arrays(pi_arr, sst_k, t0_arr) + + lnCKCD = float(np.log(CKCD)) + + efficiency = (sst_k - t0_arr) / t0_arr + valid_eff = efficiency > 0 + valid_pi = pi_arr > 0 + + lneff = np.full(efficiency.shape, np.nan, dtype=float) + lneff[valid_eff] = np.log(efficiency[valid_eff]) + + lnpi = np.full(pi_arr.shape, np.nan, dtype=float) + lnpi[valid_pi] = 2.0 * np.log(pi_arr[valid_pi]) + + lndiseq = np.full(pi_arr.shape, np.nan, dtype=float) + valid = valid_eff & valid_pi + lndiseq[valid] = lnpi[valid] - lneff[valid] - lnCKCD + + return ( + _maybe_scalar(lnpi), + _maybe_scalar(lneff), + _maybe_scalar(lndiseq), + lnCKCD, + ) + + +def pi_log_decomposition( + SSTC: Any, + MSL: Any, + P: Any, + TC: Any, + R: Any, + CKCD: float = 0.9, + *, + ascent_flag: int = 0, + diss_flag: int = 1, + V_reduc: float = 0.8, + ptop: float = 50, + miss_handle: int = 1, +) -> Dict[str, Any]: + """Run `tcpyPI.pi` and return outputs plus the Wing et al. (2015) log decomposition. + + Parameters match `tcpyPI.pi` (SST in Celsius). + + Returns + ------- + dict + Keys: `vmax`, `pmin`, `ifl`, `t0`, `otl`, `lnpi`, `lneff`, `lndiseq`, `lnCKCD`. + """ + + from .pi import pi as _pi + + vmax, pmin, ifl, t0, otl = _pi( + SSTC, + MSL, + P, + TC, + R, + CKCD=CKCD, + ascent_flag=ascent_flag, + diss_flag=diss_flag, + V_reduc=V_reduc, + ptop=ptop, + miss_handle=miss_handle, + ) + + lnpi, lneff, lndiseq, lnCKCD = log_decompose_pi( + vmax, + SSTC, + t0, + CKCD=CKCD, + sst_units="C", + ) + + return { + "vmax": vmax, + "pmin": pmin, + "ifl": ifl, + "t0": t0, + "otl": otl, + "lnpi": lnpi, + "lneff": lneff, + "lndiseq": lndiseq, + "lnCKCD": lnCKCD, + } diff --git a/tests/test_log_decomposition_api.py b/tests/test_log_decomposition_api.py new file mode 100644 index 0000000..f187ed2 --- /dev/null +++ b/tests/test_log_decomposition_api.py @@ -0,0 +1,97 @@ +import numpy as np + +import tcpyPI +from tcpyPI import utilities + + +def test_log_decompose_pi_matches_utilities_scalar(): + out = tcpyPI.log_decompose_pi(70, 300, 200, CKCD=0.9, sst_units="K") + expected = utilities.decompose_pi(70, 300, 200, CKCD=0.9) + np.testing.assert_allclose(out, expected, rtol=0, atol=1e-12) + + +def test_log_decompose_pi_sst_units_celsius(): + out_k = tcpyPI.log_decompose_pi(70, 300, 200, CKCD=0.9, sst_units="K") + out_c = tcpyPI.log_decompose_pi(70, 26.85, 200, CKCD=0.9, sst_units="C") + np.testing.assert_allclose(out_k, out_c, rtol=0, atol=1e-12) + + +def test_log_decompose_pi_vectorized_inputs(): + pi = np.array([70.0, 0.0]) + lnpi, lneff, lndiseq, lnCKCD = tcpyPI.log_decompose_pi( + pi, 300.0, 200.0, CKCD=0.9, sst_units="K" + ) + + assert np.isfinite(lnpi[0]) + assert np.isfinite(lneff[0]) + assert np.isfinite(lndiseq[0]) + + assert np.isnan(lnpi[1]) + assert np.isfinite(lneff[1]) + assert np.isnan(lndiseq[1]) + + assert np.isfinite(lnCKCD) + + +def test_pi_log_decomposition_smoke_and_consistency(): + # Example from `tcpyPI.pi` docstring + SSTC = 30 + MSL = 1010 + level_data = np.array( + [ + [1000, 28, 18], + [975, 25, 18], + [950, 24, 16], + [925, 23, 13], + [900, 22, 12], + [875, 20, 11], + [850, 19, 10], + [825, 18, 10], + [800, 16, 9], + [775, 15, 8], + [750, 13, 7], + [700, 11, 4], + [650, 8, 3], + [600, 5, 1.7], + [550, 2, 1.2], + [500, -2, 1.7], + [450, -6, 0.7], + [400, -11, 0.2], + [350, -18, 0.15], + [300, -27, 0.10], + [250, -37, 0.11], + [225, -43, 0.08], + [200, -49, 0.05], + [175, -57, 0.03], + [150, -65, 0.014], + [125, -73, 0.005], + [100, -79, 0.003], + [70, -73, 0.002], + [50, -64, 0.002], + ] + ) + P = level_data[:, 0] + TC = level_data[:, 1] + R = level_data[:, 2] + + out = tcpyPI.pi_log_decomposition(SSTC, MSL, P, TC, R, CKCD=0.9) + + for key in ["vmax", "pmin", "ifl", "t0", "otl", "lnpi", "lneff", "lndiseq", "lnCKCD"]: + assert key in out + + vmax, pmin, ifl, t0, otl = tcpyPI.pi(SSTC, MSL, P, TC, R, CKCD=0.9) + assert out["vmax"] == vmax + assert out["pmin"] == pmin + assert out["ifl"] == ifl + assert out["t0"] == t0 + assert out["otl"] == otl + + expected_lnpi, expected_lneff, expected_lndiseq, expected_lnCKCD = tcpyPI.log_decompose_pi( + vmax, SSTC, t0, CKCD=0.9, sst_units="C" + ) + np.testing.assert_allclose( + [out["lnpi"], out["lneff"], out["lndiseq"], out["lnCKCD"]], + [expected_lnpi, expected_lneff, expected_lndiseq, expected_lnCKCD], + rtol=0, + atol=1e-12, + ) From 9ca8b8466ca4cd33cf4928be19a5ec56e76c8ac6 Mon Sep 17 00:00:00 2001 From: dgilford Date: Tue, 13 Jan 2026 15:30:12 +0000 Subject: [PATCH 2/3] Docs: document log decomposition --- .gitignore | 5 +++++ README.md | 11 +++++++++++ 2 files changed, 16 insertions(+) diff --git a/.gitignore b/.gitignore index ec5a025..78bf118 100644 --- a/.gitignore +++ b/.gitignore @@ -7,6 +7,7 @@ dist/ # Jupyter .ipynb_checkpoints +.pytest_cache/ # Environment .env @@ -26,3 +27,7 @@ Thumbs.db # Data *.nc !data/*.nc + +# Local agent/dev artifacts (not tracked) +AGENTS.md +notebooks/recreate_numerical_instability.ipynb diff --git a/README.md b/README.md index e87c98b..0e36686 100644 --- a/README.md +++ b/README.md @@ -69,6 +69,17 @@ Users who want to apply the PI calculation to a set of local environmental condi (VMAX,PMIN,IFL,TO,LNB)=pi(SST,MSL,P,T,R) ``` +### Log Decomposition (Wing et al. 2015) + +tcpyPI provides a simple API for log-decomposing PI into efficiency, disequilibrium, and `Ck/Cd` terms (where `lnpi = ln(V^2)`): + +```python +from tcpyPI import log_decompose_pi, pi_log_decomposition + +out = pi_log_decomposition(SSTC, MSL, P, TC, R, CKCD=0.9) # SSTC/TC in Celsius +lnpi, lneff, lndiseq, lnCKCD = log_decompose_pi(out["vmax"], SSTC, out["t0"], CKCD=0.9, sst_units="C") +``` + ### Running a pyPI Sample Included in the pyPI release is a sample script [run_sample.py](run_sample.py) which runs global sample data from MERRA2 (in 2004) through pi.py, vectorizes the output, and performs several simple analyses. To run, simply: From 508fee2cb7aadfb8f06b7fcb631fc4ad1bb339ec Mon Sep 17 00:00:00 2001 From: dgilford Date: Tue, 13 Jan 2026 15:45:07 +0000 Subject: [PATCH 3/3] Fix Python 3.7 typing in analyses --- src/tcpyPI/analyses.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/tcpyPI/analyses.py b/src/tcpyPI/analyses.py index 2261b7d..290bca0 100644 --- a/src/tcpyPI/analyses.py +++ b/src/tcpyPI/analyses.py @@ -12,15 +12,15 @@ from __future__ import annotations -from typing import Any, Dict, Literal, Tuple +from typing import Any, Dict, Tuple, Union import numpy as np -_SSTUnits = Literal["C", "K"] +_SSTUnits = str -def _maybe_scalar(value: np.ndarray | float) -> Any: +def _maybe_scalar(value: Union[np.ndarray, float]) -> Any: return float(value) if np.ndim(value) == 0 else value