Fix RampInterpolator double-counting and calc_knownRp pressure conversion

lyopronto/calc_knownRp.py

@@ -63,7 +63,7 @@ def dry(vial,product,ht,Pchamber,Tshelf,dt):
     if Pch_t.max_setpt() > functions.Vapor_pressure(Tsh_t.max_setpt()):
         warn("Chamber pressure setpoint exceeds vapor pressure at shelf temperature " +\
         "setpoint(s). Drying cannot proceed.")
-        return np.array([[0.0, Tsh_t(0), Tsh_t(0), Tsh_t(0), Pch_t(0), 0.0, 0.0]])
+        return np.array([[0.0, Tsh_t(0), Tsh_t(0), Tsh_t(0), Pch_t(0) * constant.Torr_to_mTorr, 0.0, 0.0]])
 
     inputs = (vial, product, ht, Pch_t, Tsh_t, dt, Lpr0)
 

lyopronto/functions.py

@@ -287,14 +287,21 @@ def __init__(self, rampspec, count_ramp_against_dt=True):
             times = np.array([0.0, self.dt_setpt[0] / constant.hr_To_min])
         # Older logic: setpoint_dt includes the ramp time.
         # Kept for backward compatibility, but add a check if insufficient time allowed for ramp
-        if count_ramp_against_dt: 
+        if count_ramp_against_dt:
+            # If there is no "init" value, dt_setpt[0] was already consumed.
+            has_init = "init" in rampspec
             for i in range(1, len(self.setpt)):
-                # If less dt_setpt than setpt provided, repeat the last dt
-                totaltime = self.dt_setpt[min(len(self.dt_setpt)-1, i-1)] / constant.hr_To_min
+                # If fewer dt_setpt than setpt provided, repeat the last dt
+                dt_idx = i - 1 if has_init else i
+                totaltime = self.dt_setpt[min(len(self.dt_setpt) - 1, dt_idx)] / constant.hr_To_min
                 ramptime = abs((self.setpt[i] - self.setpt[i-1]) / self.ramp_rate) / constant.hr_To_min
                 holdtime = totaltime - ramptime
-                if ramptime > holdtime:
-                    warn(f"Ramp time from {self.setpt[i-1]:.2e} to {self.setpt[i]:.2e} exceeds total time for setpoint change, {totaltime}.")
+                if holdtime < 0:
+                    warn(f"Ramp time ({ramptime * constant.hr_To_min:.1f} min) from "
+                         f"{self.setpt[i-1]:.2e} to {self.setpt[i]:.2e} exceeds "
+                         f"total stage time ({totaltime * constant.hr_To_min:.1f} min). "
+                         f"Clamping hold time to 0.")
+                    holdtime = 0.0
                 times = np.append(times, [ramptime, holdtime])
         else:
         # Newer logic: setpoint_dt applies *after* the ramp is complete.

lyopronto/opt_Pch_Tsh.py

@@ -38,7 +38,7 @@ def dry(vial,product,ht,Pchamber,Tshelf,dt,eq_cap,nVial):
     percent_dried = Lck/Lpr0*100.0        # Percent dried
 
     # Initial chamber pressure
-    P0 = 0.1    # Initial guess for chamber pressure [Torr]
+    P0 = (Pchamber['min'] + Pchamber.get('max', Pchamber['min']*3))/2.0
 
     # Initial product and shelf temperatures
     T0=product['T_pr_crit']   # [degC]
@@ -64,7 +64,7 @@ def fun(x):
             {'type':'ineq','fun':lambda x: functions.Ineq_Constraints(x[0],x[1],product['T_pr_crit'],x[2],eq_cap['a'],eq_cap['b'],nVial)[0]},  # equipment capability inequlity
             {'type':'ineq','fun':lambda x: functions.Ineq_Constraints(x[0],x[1],product['T_pr_crit'],x[2],eq_cap['a'],eq_cap['b'],nVial)[1]})  # maximum product temperature inequality
         # Bounds for the unknowns
-        bnds = ((Pchamber['min'],None),(None,None),(None,None),(Tshelf['min'],Tshelf['max']),(None,None),(None,None),(None,None))
+        bnds = ((Pchamber['min'],Pchamber.get('max', None)),(None,None),(None,None),(Tshelf['min'],Tshelf['max']),(None,None),(None,None),(None,None))
         # Minimize the objective function i.e. maximize the sublimation rate
         res = sp.minimize(fun,x0,bounds = bnds, constraints = cons)
         [Pch,dmdt,Tbot,Tsh,Psub,Tsub,Kv] = res['x']    # Results [Torr], [kg/hr], [degC], [degC], [Torr], [degC], [cal/s/K/cm^2]

tests/test_calc_knownRp.py

@@ -410,3 +410,27 @@ def test_flux_profile_non_monotonic(self, reference_case):
         late_stage = flux[int(len(flux) * 0.8) :]
         assert np.all(np.diff(late_stage) <= 0.0), "Flux should decrease in late stage"
 
+    def test_early_return_pressure_in_mtorr(self, knownRp_standard_setup):
+        """Regression test for the bug where Pch_t(0) was returned without the
+        * constant.Torr_to_mTorr conversion, making column 4 three orders of
+        magnitude too small compared to normal output.
+        """
+        vial, product, ht, _, _, dt = knownRp_standard_setup
+
+        # Set chamber pressure very high so it exceeds vapor pressure at Tsh
+        Pchamber_high = {"setpt": [10.0], "dt_setpt": [1800.0], "ramp_rate": 0.5}
+        Tshelf = {"init": -40.0, "setpt": [-35.0], "dt_setpt": [1800.0], "ramp_rate": 1.0}
+
+        with pytest.warns(UserWarning, match="Chamber pressure"):
+            output = calc_knownRp.dry(vial, product, ht, Pchamber_high, Tshelf, dt)
+
+        assert output.shape == (1, 7), "Early return should produce exactly one row"
+
+        # Column 4 is chamber pressure in mTorr
+        Pch_mTorr = output[0, 4]
+        Pch_Torr = 10.0  # The setpoint we passed in
+
+        assert Pch_mTorr == Pch_Torr * constant.Torr_to_mTorr, (
+            f"Early-return pressure should be in mTorr ({Pch_Torr * constant.Torr_to_mTorr}), "
+            f"got {Pch_mTorr}"
+        )

tests/test_freezing.py

@@ -119,7 +119,7 @@ def test_freezing_basics(self, freezing_params):
         check_max_time(results, Tshelf, dt)
         assert results[-1, 1] == pytest.approx(Tshelf["setpt"][-1])
         # Since default setup has long hold, product should approach shelf
-        assert results[-1, 2] == pytest.approx(results[-1, 2], abs=0.1)
+        assert results[-1, 2] == pytest.approx(Tshelf["setpt"][-1], abs=0.1)
 
     def test_freezing_cin(self, freezing_params):
         """Test that freezing with imitated controlled ice nucleation is physically

tests/test_functions.py

@@ -643,6 +643,57 @@ def test_ramp_interpolator_insufficient_dt(self):
         with pytest.warns(UserWarning, match="Ramp"):
             functions.RampInterpolator(Tshelf, count_ramp_against_dt=True)
 
+    def test_ramp_interpolator_noinit_multisetpt_no_double_count(self):
+        """Test that the no-init path does not double-count dt_setpt[0].
+
+        Regression test for https://github.com/LyoHUB/LyoPRONTO/issues/17.
+        Without the fix, the first call to dt_setpt consumed index 0 during
+        initialization (times starts with [0, dt_setpt[0]]) and in the loop
+        at i=1.
+        """
+        Pchamber = {
+            "setpt": [0.1, 0.2],
+            "dt_setpt": [60, 120],   # 1 hr for first stage, 2 hr for second
+            "ramp_rate": 1.0,
+        }
+        ramp = functions.RampInterpolator(Pchamber, count_ramp_against_dt=True)
+
+        # First stage: hold at 0.1 for 1 hr (from initialization)
+        # Second stage: ramp from 0.1->0.2 (ramptime = 0.1/1.0 / 60 hr),
+        #               then hold for remainder of dt_setpt[1]=120 min = 2 hr
+        ramptime_2 = abs(0.2 - 0.1) / 1.0 / constant.hr_To_min
+        holdtime_2 = 120 / constant.hr_To_min - ramptime_2
+
+        expected_times = np.cumsum([0.0, 60 / constant.hr_To_min, ramptime_2, holdtime_2])
+        np.testing.assert_allclose(ramp.times, expected_times)
+
+        # Times should be monotonically non-decreasing (holdtime can be 0)
+        assert np.all(np.diff(ramp.times) >= 0), "Times must be monotonically non-decreasing"
+
+        # Values at boundaries
+        assert ramp(0.0) == pytest.approx(0.1)
+        assert ramp(ramp.times[-1]) == pytest.approx(0.2)
+        assert ramp(100.0) == pytest.approx(0.2)
+
+    def test_ramp_interpolator_holdtime_clamped_to_zero(self):
+        """Test that negative holdtime is clamped to 0 with a warning.
+
+        Regression test for https://github.com/LyoHUB/LyoPRONTO/issues/17.
+        When ramp time exceeds the total stage time, holdtime goes negative
+        which would produce non-monotonic cumulative times.
+        """
+        Tshelf = {
+            "init": 0.0,
+            "setpt": [100.0],          # Need to ramp 100 deg
+            "dt_setpt": [1],           # Only 1 minute allowed
+            "ramp_rate": 0.1,          # 0.1 deg/min -> takes 1000 min
+        }
+        with pytest.warns(UserWarning, match="Ramp time"):
+            ramp = functions.RampInterpolator(Tshelf, count_ramp_against_dt=True)
+
+        # holdtime should be clamped to 0, not negative
+        assert np.all(np.diff(ramp.times) >= 0), "Times must be monotonically non-decreasing"
+
     def test_ramp_interpolator_out_of_bounds(self):
         """Test RampInterpolator behavior outside defined time range."""
         Tshelf = {

tests/test_main.py

@@ -1,10 +1,23 @@
 """Tests for the high-level API (formerly main.py)."""
 
-from contextlib import chdir
+import os
+from contextlib import contextmanager
 import pytest
 import numpy as np
 from lyopronto import *
 
+
+@contextmanager
+def chdir(path):
+    """contextlib.chdir was added in Python 3.11; this provides equivalent
+    functionality for older supported versions."""
+    old = os.getcwd()
+    os.chdir(path)
+    try:
+        yield
+    finally:
+        os.chdir(old)
+
 class TestHighLevelAPI:
     """Tests for the high-level API functions in lyopronto.high_level."""
 

tests/test_opt_Pch.py

@@ -46,7 +46,7 @@ def opt_pch_consistency(output, setup):
     assert np.all(Pch_values >= Pchamber["min"] * constant.Torr_to_mTorr), (
         "Pressure should be >= min bound"
     )
-    if hasattr(Pchamber, "max"):
+    if "max" in Pchamber:
         assert np.all(Pch_values <= Pchamber["max"] * constant.Torr_to_mTorr), (
             "Pressure should be <= max bound"
         )
@@ -355,12 +355,13 @@ def test_opt_pch_reference(self, repo_root, opt_pch_reference_inputs):
         # Instead, check that output is reasonable and matches or exceeds the performance.
         opt_pch_consistency(output, opt_pch_reference_inputs)
         assert_complete_drying(output)
-        # Drying time should be equal to or better than reference
+        # Drying time should be equal to or better than reference (with small tolerance
+        # for floating-point differences across Python versions)
         drying_time_ref = output_ref[-1, 0]
         drying_time = output[-1, 0]
-        assert drying_time <= drying_time_ref, (
-            f"Drying time {drying_time:.2f} hr should be <= reference "
-            + f"{drying_time_ref:.2f} hr"
+        assert drying_time <= drying_time_ref + 1e-3, (
+            f"Drying time {drying_time:.6f} hr should be <= reference "
+            + f"{drying_time_ref:.6f} hr"
         )
         # array_compare = np.isclose(output, output_ref, atol=1e-3)
         # assert array_compare.all(), (

tests/test_opt_Pch_Tsh.py

@@ -102,7 +102,7 @@ def opt_both_consistency(output, setup):
     assert np.all(Pch_values >= Pchamber["min"] * constant.Torr_to_mTorr), (
         "Pressure should be >= min bound"
     )
-    if hasattr(Pchamber, "max"):
+    if "max" in Pchamber:
         assert np.all(Pch_values <= Pchamber["max"] * constant.Torr_to_mTorr), (
             "Pressure should be <= max bound"
         )

tests/test_opt_Tsh.py

@@ -33,7 +33,7 @@ def opt_tsh_consistency(output, setup):
     # Chamber pressure should start at first setpoint
     # Note: May not reach final setpoint if drying completes first
     Pch_values = output[:, 4]
-    Pch_check = functions.RampInterpolator(Pchamber)(output[:, 0])
+    Pch_check = functions.RampInterpolator(Pchamber)(output[:, 0]) * constant.Torr_to_mTorr
     np.testing.assert_allclose(Pch_values, Pch_check, atol=0.1)
 
     # Shelf temperature (column 3) should start at init
@@ -49,7 +49,7 @@ def opt_tsh_consistency(output, setup):
     assert np.all(Tsh_values >= Tshelf["min"]), (
         "Shelf temperature should be >= min bound"
     )
-    if hasattr(Tshelf, "max"):
+    if "max" in Tshelf:
         assert np.all(Tsh_values <= Tshelf["max"]), (
             "Shelf temperature should be <= max bound"
         )