Add temperature estimation #422

add temperature estimation

make estimation function independently importable

Author: peterhollender

src/openlifu/plan/solution.py

@@ -22,6 +22,7 @@
     find_centroid,
     get_beamwidth,
     get_mask,
+    model_tx_temperature_rise,
 )
 from openlifu.sim import SimSetup, run_simulation
 from openlifu.util.annotations import OpenLIFUFieldData
@@ -30,6 +31,7 @@
 from openlifu.util.units import getunitconversion, rescale_coords, rescale_data_arr
 from openlifu.xdc import Transducer
 
+logger = logging.getLogger(__name__)
 
 def _construct_nc_filepath_from_json_filepath(json_filepath:Path) -> Path:
     """Construct a default filepath to netCDF file given filepath to associated solution json file."""
@@ -342,10 +344,45 @@ def analyze(self,
         solution_analysis.TIC = np.mean(TIC)
         solution_analysis.voltage_V = self.voltage
         solution_analysis.power_W = np.mean(power_W)
-
+        solution_analysis.estimated_tx_temperature_rise_C = self.estimate_tx_temperature_rise(
+            t_sec=solution_analysis.sequence_duration_s,
+        )
         solution_analysis.param_constraints = param_constraints
         return solution_analysis
 
+    def estimate_tx_temperature_rise(self,
+                                  t_sec: float,
+                                  T0_degC: float = 30.0):
+
+        """
+        Standalone temperature prediction function for thermal modeling.
+
+        Based on physics-based power decay gradient model fitted from experimental data.
+        Model: dT/dt = (t + t_shift)^(-n) + C
+
+        Parameters:
+        -----------
+        self: Solution
+            The treatment solution containing voltage and pulse information.
+        t_sec: float
+            Time in seconds for which to predict temperature rise.
+        T0_degC: float
+            Initial temperature in Celsius. Default is 30.0°C.
+
+        Returns:
+        --------
+        float
+            Predicted temperature rise in Celsius
+        """
+        return model_tx_temperature_rise(
+            voltage=self.voltage,
+            t_sec=t_sec,
+            duty_cycle=self.get_sequence_dutycycle(),
+            apodization_fraction=np.mean(self.apodizations),
+            frequency_kHz=self.pulse.frequency / 1e3,
+            T0_degC=T0_degC,
+        )
+
     def compute_scaling_factors(
             self,
             focal_pattern: FocalPattern,

src/openlifu/plan/solution_analysis.py

@@ -1,6 +1,7 @@
 from __future__ import annotations
 
 import json
+import logging
 from dataclasses import dataclass, field
 from typing import Annotated, Any, Dict, Tuple, Type
 
@@ -13,6 +14,8 @@
 from openlifu.util.dict_conversion import DictMixin
 from openlifu.util.units import getunitconversion, getunittype
 
+logger = logging.getLogger(__name__)
+
 DEFAULT_ORIGIN = np.zeros(3)
 
 PARAM_FORMATS = {
@@ -46,7 +49,8 @@
     "power_W": [None, "0.2f", "W", "Emitted Power"],
     "duty_cycle_pulse_train_pct": [None, "0.1f", "%", "Pulse Train Duty Cycle"],
     "duty_cycle_sequence_pct": [None, "0.1f", "%", "Sequence Duty Cycle"],
-    "sequence_duration_s": [None, "0.0f", "s", "Sequence Duration"],}
+    "sequence_duration_s": [None, "0.0f", "s", "Sequence Duration"],
+    "estimated_tx_temperature_rise_C": [None, "0.2f", "°C", "Estimated TX Temperature Rise"]}
 
 @dataclass
 class SolutionAnalysis(DictMixin):
@@ -143,6 +147,8 @@ class SolutionAnalysis(DictMixin):
     sequence_duration_s: Annotated[float | None, OpenLIFUFieldData("Sequence duration (s)", "Total duration of the sequence (s)")] = None
     """Total duration of the sequence (s)"""
 
+    estimated_tx_temperature_rise_C: Annotated[float | None, OpenLIFUFieldData("Estimated TX temperature rise (°C)", "Estimated temperature rise (assume 30°C baseline)")] = None
+    """Estimated temperature rise in Celsius (assume 30°C baseline)"""
 
     param_constraints: Annotated[Dict[str, ParameterConstraint], OpenLIFUFieldData("Parameter constraints", None)] = field(default_factory=dict)
     """TODO: Add description"""
@@ -576,3 +582,69 @@ def get_beamwidth(
     negoff, posoff = get_beam_bounds(da, focus=focus, dim=dim, cutoff=float(cutoff), origin=origin, min_offset=min_offset, max_offset=max_offset)
     bw = posoff - negoff
     return bw
+
+def model_tx_temperature_rise(voltage: float,
+                              t_sec: float,
+                              duty_cycle: float=1.0,
+                              apodization_fraction: float=1.0,
+                              frequency_kHz: float=400.0,
+                              T0_degC: float = 30.0):
+    """
+    Temperature prediction function for thermal modeling.
+
+    Based on physics-based power decay gradient model fitted from experimental data.
+    Model: dT/dt = (t + t_shift)^(-n) + C
+
+    Parameters:
+    -----------
+    voltage: float
+        Voltage to use when running sonication.
+    duty_cycle: float
+        Duty cycle of the sonication.
+    apodization_fraction: float
+        Fraction of apodization applied.
+    frequency_kHz: float
+        Frequency in kHz.
+    t_sec: float
+        Time in seconds for which to predict temperature rise.
+    T0_degC: float
+        Initial temperature in Celsius. Default is 30.0°C.
+
+    Returns:
+    --------
+    float
+        Predicted temperature rise in Celsius
+    """
+    P = voltage**2 * duty_cycle
+    P = P * apodization_fraction
+    t = t_sec
+    T0 = T0_degC
+
+    if T0 < 20 or T0 > 40:
+        logger.warning("Initial temperature T0 must be between 20 and 40 degrees Celsius for the model to be valid.")
+
+    if P < 50 or P > 500:
+        logger.warning("Squared Voltage must be between 50 and 500 V^2 for the model to be valid.")
+
+    if t < 1 or t > 600:
+        logger.warning("Time t must be between 1 and 600 seconds for the model to be valid.")
+
+    if frequency_kHz < 380 or frequency_kHz > 420:
+        logger.warning("Frequency must be between 380 and 420 kHz for the model to be valid.")
+
+    # Predict power law parameters using polynomial regression (degree 2)
+    n = (2.131832 + -0.003475*P + -0.044916*T0 +
+        0.000002*P*P + 0.000049*P*T0 + 0.000474*T0*T0)
+
+    t_shift = (1.074525 + 0.101532*P + -0.097741*T0 +
+            -0.000080*P*P + -0.001891*P*T0 + 0.005975*T0*T0)
+
+    C = (0.051572 + -0.000072*P + -0.001668*T0 +
+        -0.000000*P*P + 0.000004*P*T0 + 0.000007*T0*T0)
+
+    # Apply integrated power decay gradient model
+    t_adj = t + t_shift
+    integral_term = (t_adj**(1-n) - t_shift**(1-n)) / (1 - n)
+    temperature_rise_C = integral_term + C * t
+
+    return temperature_rise_C

tests/test_solution_analysis.py

@@ -1,9 +1,15 @@
 from __future__ import annotations
 
+import logging
+
 import pytest
 
 from openlifu.plan.param_constraint import ParameterConstraint
-from openlifu.plan.solution_analysis import SolutionAnalysis, SolutionAnalysisOptions
+from openlifu.plan.solution_analysis import (
+    SolutionAnalysis,
+    SolutionAnalysisOptions,
+    model_tx_temperature_rise,
+)
 
 # ---- Tests for SolutionAnalysis ----
 
@@ -80,3 +86,148 @@ def test_to_dict_from_dict_solution_analysis_options(example_solution_analysis_o
     options_dict = example_solution_analysis_options.to_dict()
     new_options = SolutionAnalysisOptions.from_dict(options_dict)
     assert new_options == example_solution_analysis_options
+
+
+# ---- Tests for model_tx_temperature_rise ----
+
+# Valid mid-range parameters used as a baseline throughout the tests.
+# P = voltage^2 * duty_cycle = 14^2 * 1.0 = 196 V^2 (valid range: 50-500)
+_BASE_VOLTAGE = 14.0       # V
+_LOW_VOLTAGE = 8.0         # V  → P = 64 V^2 (near low end of valid range)
+_HIGH_VOLTAGE = 21.0       # V  → P = 441 V^2 (near high end of valid range)
+_BASE_T0 = 30.0            # °C (mid-range of valid 20-40 °C)
+_BASE_FREQ = 400.0         # kHz (centre of valid 380-420 kHz)
+
+
+def test_low_voltage_less_heating_than_high_voltage():
+    """Higher voltage (more power) should produce a greater temperature rise."""
+    t = 120.0  # mid-range time, well within valid 1-600 s
+    rise_low = model_tx_temperature_rise(_LOW_VOLTAGE, t, T0_degC=_BASE_T0, frequency_kHz=_BASE_FREQ)
+    rise_high = model_tx_temperature_rise(_HIGH_VOLTAGE, t, T0_degC=_BASE_T0, frequency_kHz=_BASE_FREQ)
+    assert rise_low < rise_high, (
+        f"Expected lower voltage ({_LOW_VOLTAGE} V, rise={rise_low:.4f} °C) to produce "
+        f"less heating than higher voltage ({_HIGH_VOLTAGE} V, rise={rise_high:.4f} °C)"
+    )
+
+
+def test_little_heating_right_after_start():
+    """Temperature rise at t=1 s (just started) should be much less than at t=300 s."""
+    rise_early = model_tx_temperature_rise(_BASE_VOLTAGE, t_sec=1.0, T0_degC=_BASE_T0)
+    rise_later = model_tx_temperature_rise(_BASE_VOLTAGE, t_sec=300.0, T0_degC=_BASE_T0)
+    assert rise_early < rise_later, (
+        f"Expected rise at t=1 s ({rise_early:.4f} °C) to be less than rise at t=300 s ({rise_later:.4f} °C)"
+    )
+    # Additionally confirm the early rise is small in absolute terms
+    assert rise_early < 5.0, (
+        f"Expected temperature rise at t=1 s to be < 5 °C, got {rise_early:.4f} °C"
+    )
+
+
+def test_temperature_rises_monotonically_with_time():
+    """Temperature rise must be strictly increasing across a span of time points."""
+    times = [1.0, 10.0, 60.0, 180.0, 360.0, 600.0]
+    rises = [model_tx_temperature_rise(_BASE_VOLTAGE, t, T0_degC=_BASE_T0) for t in times]
+    for i in range(len(rises) - 1):
+        assert rises[i] < rises[i + 1], (
+            f"Temperature rise not monotonically increasing: "
+            f"rise[{times[i]} s]={rises[i]:.4f} >= rise[{times[i+1]} s]={rises[i+1]:.4f}"
+        )
+
+
+def test_temperature_rises_monotonically_with_voltage():
+    """Temperature rise must increase with voltage (at fixed time and other params)."""
+    # Voltages chosen so that P = V^2 stays within the valid 50-500 V^2 range
+    voltages = [8.0, 11.0, 14.0, 17.0, 21.0]
+    t = 120.0
+    rises = [model_tx_temperature_rise(v, t, T0_degC=_BASE_T0) for v in voltages]
+    for i in range(len(rises) - 1):
+        assert rises[i] < rises[i + 1], (
+            f"Temperature rise not monotonically increasing with voltage: "
+            f"rise[{voltages[i]} V]={rises[i]:.4f} >= rise[{voltages[i+1]} V]={rises[i+1]:.4f}"
+        )
+
+
+def test_lower_duty_cycle_produces_less_heating():
+    """Reducing duty cycle reduces effective power and therefore temperature rise."""
+    t = 120.0
+    voltage = _BASE_VOLTAGE
+    rise_full = model_tx_temperature_rise(voltage, t, duty_cycle=1.0, T0_degC=_BASE_T0)
+    rise_half = model_tx_temperature_rise(voltage, t, duty_cycle=0.5, T0_degC=_BASE_T0)
+    assert rise_half < rise_full, (
+        f"Expected 50 % duty cycle ({rise_half:.4f} °C) to produce less heating "
+        f"than 100 % duty cycle ({rise_full:.4f} °C)"
+    )
+
+
+def test_lower_apodization_produces_less_heating():
+    """Partial apodization reduces effective power and therefore temperature rise."""
+    t = 120.0
+    rise_full = model_tx_temperature_rise(_BASE_VOLTAGE, t, apodization_fraction=1.0, T0_degC=_BASE_T0)
+    rise_half = model_tx_temperature_rise(_BASE_VOLTAGE, t, apodization_fraction=0.5, T0_degC=_BASE_T0)
+    assert rise_half < rise_full, (
+        f"Expected apodization=0.5 ({rise_half:.4f} °C) to produce less heating "
+        f"than apodization=1.0 ({rise_full:.4f} °C)"
+    )
+
+
+@pytest.mark.parametrize("bad_T0", [19.9, 40.1])
+def test_warning_emitted_for_T0_out_of_range(bad_T0, caplog):
+    """A warning must be logged when T0 is outside the valid 20-40 °C range."""
+    with caplog.at_level(logging.WARNING, logger="openlifu.plan.solution_analysis"):
+        model_tx_temperature_rise(_BASE_VOLTAGE, t_sec=60.0, T0_degC=bad_T0)
+    assert any("T0" in record.message or "temperature" in record.message.lower()
+               for record in caplog.records), (
+        f"Expected a warning about T0 out of range for T0={bad_T0} °C"
+    )
+
+
+@pytest.mark.parametrize(("bad_voltage","bad_duty_cycle"), [
+    (6.0, 1.0),   # P = 36 < 50
+    (25.0, 1.0),  # P = 625 > 500
+])
+def test_warning_emitted_for_power_out_of_range(bad_voltage, bad_duty_cycle, caplog):
+    """A warning must be logged when the squared-voltage power is outside 50-500 V^2."""
+    with caplog.at_level(logging.WARNING, logger="openlifu.plan.solution_analysis"):
+        model_tx_temperature_rise(bad_voltage, t_sec=60.0, duty_cycle=bad_duty_cycle, T0_degC=_BASE_T0)
+    assert any("voltage" in record.message.lower() or "squared" in record.message.lower() or "v^2" in record.message.lower()
+               for record in caplog.records), (
+        f"Expected a warning about power out of range for voltage={bad_voltage} V"
+    )
+
+
+@pytest.mark.parametrize("bad_time", [0.5, 601.0])
+def test_warning_emitted_for_time_out_of_range(bad_time, caplog):
+    """A warning must be logged when t is outside the valid 1-600 s range."""
+    with caplog.at_level(logging.WARNING, logger="openlifu.plan.solution_analysis"):
+        model_tx_temperature_rise(_BASE_VOLTAGE, t_sec=bad_time, T0_degC=_BASE_T0)
+    assert any("time" in record.message.lower() or "seconds" in record.message.lower()
+               for record in caplog.records), (
+        f"Expected a warning about time out of range for t={bad_time} s"
+    )
+
+
+@pytest.mark.parametrize("bad_freq", [379.9, 420.1])
+def test_warning_emitted_for_frequency_out_of_range(bad_freq, caplog):
+    """A warning must be logged when frequency is outside the valid 380-420 kHz range."""
+    with caplog.at_level(logging.WARNING, logger="openlifu.plan.solution_analysis"):
+        model_tx_temperature_rise(_BASE_VOLTAGE, t_sec=60.0, frequency_kHz=bad_freq, T0_degC=_BASE_T0)
+    assert any("frequency" in record.message.lower() or "khz" in record.message.lower()
+               for record in caplog.records), (
+        f"Expected a warning about frequency out of range for freq={bad_freq} kHz"
+    )
+
+
+def test_no_warnings_for_valid_inputs(caplog):
+    """No warnings should be emitted when all inputs are within their valid ranges."""
+    with caplog.at_level(logging.WARNING, logger="openlifu.plan.solution_analysis"):
+        model_tx_temperature_rise(
+            voltage=_BASE_VOLTAGE,
+            t_sec=60.0,
+            duty_cycle=1.0,
+            apodization_fraction=1.0,
+            frequency_kHz=_BASE_FREQ,
+            T0_degC=_BASE_T0,
+        )
+    assert len(caplog.records) == 0, (
+        f"Unexpected warnings for valid inputs: {[r.message for r in caplog.records]}"
+    )