Decomposition terms update

.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

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:

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",

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=[
             [], [], [],
         ],

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

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, Tuple, Union
+
+import numpy as np
+
+
+_SSTUnits = str
+
+
+def _maybe_scalar(value: Union[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,
+    }

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,
+    )