fix: guard against empty peak array in get_ignition_delay

examples/h2o2-shocktube/h2o2-results.yaml

@@ -0,0 +1,56 @@
+average deviation function: 7.6446209680420125
+average error function: 58.91886752410112
+datasets:
+- absolute deviation: 7.6446209680420125
+  datapoints:
+  - composition:
+    - {InChI: 1S/H2/h1H, amount: '0.00444', species-name: H2}
+    - {InChI: 1S/O2/c1-2, amount: '0.00556', species-name: O2}
+    - {InChI: 1S/Ar, amount: '0.99', species-name: Ar}
+    composition type: mole fraction
+    experimental ignition delay: 0.00047154 second
+    pressure: 220000.0 pascal
+    simulated ignition delay: 0.0010007455162881925 second
+    temperature: 1164.48 kelvin
+  - composition:
+    - {InChI: 1S/H2/h1H, amount: '0.00444', species-name: H2}
+    - {InChI: 1S/O2/c1-2, amount: '0.00556', species-name: O2}
+    - {InChI: 1S/Ar, amount: '0.99', species-name: Ar}
+    composition type: mole fraction
+    experimental ignition delay: 0.00044803 second
+    pressure: 220000.0 pascal
+    simulated ignition delay: 0.0009972083218715213 second
+    temperature: 1164.97 kelvin
+  - composition:
+    - {InChI: 1S/H2/h1H, amount: '0.00444', species-name: H2}
+    - {InChI: 1S/O2/c1-2, amount: '0.00556', species-name: O2}
+    - {InChI: 1S/Ar, amount: '0.99', species-name: Ar}
+    composition type: mole fraction
+    experimental ignition delay: 0.00029157 second
+    pressure: 220000.0 pascal
+    simulated ignition delay: 0.000574010524129157 second
+    temperature: 1264.2 kelvin
+  - composition:
+    - {InChI: 1S/H2/h1H, amount: '0.00444', species-name: H2}
+    - {InChI: 1S/O2/c1-2, amount: '0.00556', species-name: O2}
+    - {InChI: 1S/Ar, amount: '0.99', species-name: Ar}
+    composition type: mole fraction
+    experimental ignition delay: 0.00020593 second
+    pressure: 220000.0 pascal
+    simulated ignition delay: 0.00042129329629920055 second
+    temperature: 1332.57 kelvin
+  - composition:
+    - {InChI: 1S/H2/h1H, amount: '0.00444', species-name: H2}
+    - {InChI: 1S/O2/c1-2, amount: '0.00556', species-name: O2}
+    - {InChI: 1S/Ar, amount: '0.99', species-name: Ar}
+    composition type: mole fraction
+    experimental ignition delay: 8.811e-05 second
+    pressure: 220000.0 pascal
+    simulated ignition delay: 0.00021169358253286696 second
+    temperature: 1519.18 kelvin
+  dataset: H2O2Ar-shocktube-example.yaml
+  dataset_id: 0
+  error function: 58.91886752410112
+  standard deviation: 0.1
+error function standard deviation: 0.0
+model: h2o2.cti

examples/h2o2-shocktube/main.py

@@ -5,10 +5,12 @@
 # Local imports
 from pyteck.eval_model import evaluate_model
 
-evaluate_model(model_name='h2o2.cti',
-               spec_keys_file='spec_keys.yaml',
-               dataset_file='dataset_file.txt',
-               data_path='data',
-               model_path='models',
-               results_path='results',
-)
+if __name__ == '__main__':
+    evaluate_model(model_name='h2o2.cti',
+                   spec_keys_file='spec_keys.yaml',
+                   dataset_file='dataset_file.txt',
+                   data_path='data',
+                   model_path='models',
+                   results_path='results',
+                   skip_validation=True,
+    )

pyteck/simulation.py

@@ -100,35 +100,40 @@ def get_ignition_delay(time, target, target_name, ignition_type):
         # Get indices of peaks
         peak_inds = detect_peaks(target, edge=None, mph=1.e-9*np.max(target))
 
-        # Get index of largest peak (overall ignition delay)
-        max_ind = peak_inds[np.argmax(target[peak_inds])]
-
-        #ign_delays = time[peak_inds[np.where((time[peak_inds[peak_inds <= max_ind]]) > 0.0)]]
-
-        ign_delays = time[peak_inds[peak_inds <= max_ind]]
+        if peak_inds.size == 0:
+            ign_delays = np.array([])
+        else:
+            # Get index of largest peak (overall ignition delay)
+            max_ind = peak_inds[np.argmax(target[peak_inds])]
+            ign_delays = time[peak_inds[peak_inds <= max_ind]]
 
     elif ignition_type == 'd/dt max':
         target = first_derivative(time, target)
         # Get indices of peaks. Set a minimum peak height of 1e-7% of the
         # maximum value to avoid noise peaks.
         peak_inds = detect_peaks(target, edge=None, mph=1.e-9*np.max(target))
 
-        # Get index of largest peak (overall ignition delay)
-        max_ind = peak_inds[np.argmax(target[peak_inds])]
-
-        ign_delays = time[peak_inds[np.where((time[peak_inds[peak_inds <= max_ind]]) > 0.0)]]
+        if peak_inds.size == 0:
+            ign_delays = np.array([])
+        else:
+            # Get index of largest peak (overall ignition delay)
+            max_ind = peak_inds[np.argmax(target[peak_inds])]
+            ign_delays = time[peak_inds[np.where((time[peak_inds[peak_inds <= max_ind]]) > 0.0)]]
 
     elif ignition_type == '1/2 max':
         # maximum value, and associated index
         max_val = np.max(target)
         peak_inds = detect_peaks(target, edge=None, mph=1.e-9*np.max(target))
-        max_ind = peak_inds[np.argmax(target[peak_inds])]
 
-        # TODO: interpolate for actual half-max value
-        # Find index associated with the 1/2 max value, but only consider
-        # points before the peak
-        half_idx = (np.abs(target[0:max_ind] - 0.5 * max_val)).argmin()
-        ign_delays = np.array([time[half_idx]])
+        if peak_inds.size == 0:
+            ign_delays = np.array([])
+        else:
+            max_ind = peak_inds[np.argmax(target[peak_inds])]
+            # TODO: interpolate for actual half-max value
+            # Find index associated with the 1/2 max value, but only consider
+            # points before the peak
+            half_idx = (np.abs(target[0:max_ind] - 0.5 * max_val)).argmin()
+            ign_delays = np.array([time[half_idx]])
 
     elif ignition_type == 'd/dt max extrapolated':
         # First need to evaluate derivative of the target
@@ -137,10 +142,13 @@ def get_ignition_delay(time, target, target_name, ignition_type):
         # Get indices of peaks, and index of largest peak, which corresponds to
         # the point of maximum deriative
         peak_inds = detect_peaks(target_derivative, edge=None, mph=1.e-9*np.max(target))
-        max_ind = peak_inds[np.argmax(target_derivative[peak_inds])]
 
-        # use slope to extrapolate to intercept with baseline value (0 by default)
-        ign_delays = np.array([time[max_ind] - (target[max_ind] / target_derivative[max_ind])])
+        if peak_inds.size == 0:
+            ign_delays = np.array([])
+        else:
+            max_ind = peak_inds[np.argmax(target_derivative[peak_inds])]
+            # use slope to extrapolate to intercept with baseline value (0 by default)
+            ign_delays = np.array([time[max_ind] - (target[max_ind] / target_derivative[max_ind])])
 
         # TODO: handle target with nonzero baseline?
     else:

pyteck_issue_get_ignition_delay_crash.md

@@ -0,0 +1,70 @@
+## Bug: `get_ignition_delay` crashes on empty peak array
+
+### Description
+
+`get_ignition_delay()` in `pyteck/simulation.py` raises `ValueError: attempt to get argmax of an empty sequence` when `detect_peaks` returns no peaks. This happens when the input signal is flat or nearly flat (e.g., from a non-igniting reactor simulation).
+
+All four ignition types are affected: `max`, `d/dt max`, `1/2 max`, and `d/dt max extrapolated`.
+
+The function already has a safety check at the bottom:
+
+```python
+if ign_delays.size == 0:
+    # ... warning + return np.array([0.0])
+```
+
+But execution never reaches it because `peak_inds[np.argmax(target[peak_inds])]` crashes first when `peak_inds` is empty.
+
+### Minimal Reproducer
+
+```python
+import numpy as np
+from pyteck.simulation import get_ignition_delay
+
+# Flat signal — no ignition occurred
+time = np.linspace(0, 0.05, 100)
+temperature = np.full_like(time, 500.0)
+
+# Crashes with: ValueError: attempt to get argmax of an empty sequence
+ign = get_ignition_delay(time, temperature, 'temperature', 'd/dt max')
+```
+
+### Expected Behavior
+
+The function should return `np.array([0.0])` (the existing "no ignition detected" sentinel), not crash.
+
+### Root Cause
+
+In each branch (`max`, `d/dt max`, etc.), `detect_peaks` may return an empty array when the signal has no meaningful peaks. The line:
+
+```python
+max_ind = peak_inds[np.argmax(target[peak_inds])]
+```
+
+then calls `np.argmax` on an empty slice, which raises `ValueError`.
+
+### Suggested Fix
+
+Add an empty-array guard after each `detect_peaks` call so that execution falls through to the existing safety check:
+
+```python
+peak_inds = detect_peaks(target, edge=None, mph=1.e-9*np.max(target))
+
+if peak_inds.size == 0:
+    ign_delays = np.array([])
+else:
+    max_ind = peak_inds[np.argmax(target[peak_inds])]
+    ign_delays = time[peak_inds[peak_inds <= max_ind]]
+```
+
+This needs to be applied in all four branches (`max`, `d/dt max`, `1/2 max`, `d/dt max extrapolated`).
+
+### Context
+
+I'm using `get_ignition_delay` as a standalone function to post-process Cantera 0-D reactor output across a large parameter grid (temperature × pressure × equivalence ratio × multiple kinetic models). The grid intentionally includes conditions where ignition does not occur, which triggers this crash.
+
+### Environment
+
+- PyTeCK development branch (commit with `get_ignition_delay` as a standalone function)
+- NumPy 1.x
+- Python 3.9

setup.py

@@ -26,9 +26,9 @@
 long_description = readme + '\n\n' + changelog + '\n\n' + citation
 
 install_requires = [
-    'pyyaml>=3.12,<4.0',
-    'pint>=0.7.2,<0.9',
-    'numpy>=1.13.0,<2.0',
+    'pyyaml>=3.12',
+    'pint>=0.7.2',
+    'numpy>=1.13.0',
     'tables',
     'pyked>=0.4.1',
     'scipy>=1.0.0',
@@ -60,6 +60,7 @@
     package_data={'pyteck': ['tests/*.xml', 'tests/*.yaml', 'tests/dataset_file.txt', 'tests/*.cti']},
     install_requires=install_requires,
     zip_safe=False,
+    python_requires='>=3.7',
 
     license='MIT License',