diff --git a/hplc/quant.py b/hplc/quant.py index 17c8e95..67fddf2 100644 --- a/hplc/quant.py +++ b/hplc/quant.py @@ -261,6 +261,12 @@ def _assign_windows( # Determine if peaks should be added. if len(known_peaks) > 0: + # Anchor the index math to the first time value of the (possibly + # cropped) chromatogram. `find_peaks` returns positional indices into + # this array, so converting a user-specified time to an index must be + # measured relative to the start time rather than assuming t = 0. + t0 = self.df[self.time_col].values[0] + # Get the enforced peak positions if type(known_peaks) == dict: _known_peaks = list(known_peaks.keys()) @@ -268,13 +274,13 @@ def _assign_windows( _known_peaks = known_peaks # Find the nearest location in the time array given the user-specified time - enforced_location_inds = ( - np.int_(np.array(_known_peaks) / self._dt) - self._crop_offset + enforced_location_inds = np.int_( + (np.array(_known_peaks) - t0) / self._dt ) # Update the user specified times with the nearest location updated_loc = np.round( - self._dt * (enforced_location_inds + self._crop_offset), + t0 + self._dt * enforced_location_inds, decimals=self._time_precision, ) if type(known_peaks) == dict: @@ -309,7 +315,7 @@ def _assign_windows( self._peak_indices = np.append(self._peak_indices, loc) if self._added_peaks is None: self._added_peaks = [] - self._added_peaks.append((loc + self._crop_offset) * self._dt) + self._added_peaks.append(t0 + loc * self._dt) if type(known_peaks) == dict: _sel_loc = updated_known_peaks[_known_peaks[i]] if "width" in _sel_loc.keys(): @@ -373,9 +379,12 @@ def _assign_windows( "window_type", ] = "interpeak" - # If more than one split ind, set up all ranges. - elif split_inds[0] != 0: - split_inds += 1 + # Otherwise, build the boundaries for every background segment. This + # must run for all multi-segment cases; a gap at the first sample + # (split_inds[0] == 0) was previously skipped, which dropped the final + # background segment and mislabeled interpeak windows. + else: + split_inds = split_inds + 1 split_inds = np.insert(split_inds, 0, 0) split_inds = np.append(split_inds, len(tidx)) @@ -693,11 +702,43 @@ def deconvolve_peaks( raise ValueError( f"Could not adjust bounds for peak at {v['location'][i]} because bound keys do not contain at least one of the following: `location`, `amplitude`, `scale`, `skew`. " ) + # Clamp the location initial guess into its bounds. The guess is + # the peak time rounded to `_time_precision`, while the bounds are + # the window's raw (unrounded) time range; rounding can push the + # guess just outside that range, which scipy rejects as + # "`x0` is infeasible". The true (unrounded) peak time always lies + # inside the window, so clamping only undoes the rounding excursion. + _loc_lo, _loc_hi = _param_bounds["location"] + p0[paridx["location"]] = min( + max(p0[paridx["location"]], _loc_lo), _loc_hi + ) + for _, val in _param_bounds.items(): bounds[0].append(val[0]) bounds[1].append(val[1]) self._param_bounds.append(_param_bounds) + # Safety net: ensure every parameter's bounds are valid and contain + # the initial guess before calling the optimizer, so a degenerate or + # inverted bound surfaces as an actionable message instead of scipy's + # opaque "`x0` is infeasible" / "lower bound ... strictly less ...". + for _j, (_g, _lo, _hi) in enumerate(zip(p0, bounds[0], bounds[1])): + _name = parorder[_j % 4] + _pk = v["location"][_j // 4] + if _lo >= _hi: + raise ValueError( + f"Invalid bounds for '{_name}' of the peak near retention " + f"time {_pk}: lower bound ({_lo:.4g}) is not less than upper " + f"bound ({_hi:.4g}). Try adjusting `param_bounds` or cropping " + f"the chromatogram to exclude edge artifacts." + ) + if not (_lo <= _g <= _hi): + raise ValueError( + f"Initial guess for '{_name}' of the peak near retention " + f"time {_pk} ({_g:.4g}) lies outside its bounds " + f"[{_lo:.4g}, {_hi:.4g}]. Try adjusting `param_bounds`." + ) + # Perform the inference popt, _ = scipy.optimize.curve_fit( self._fit_skewnorms, @@ -875,16 +916,18 @@ def fit_peaks( peak_df["peak_id"] = peak_df["peak_id"].astype(int) self.peaks = peak_df - # Compute the mixture + # Compute the mixture. Build the columns directly from the sorted peak + # table so that column `i` always corresponds to `peak_id == i + 1`. The + # previous implementation filled columns in window/detection order, which + # could disagree with the retention-time-sorted `peak_id` (e.g. when a + # large skew shifts a fitted retention time across a neighbor) and + # mislabel per-peak traces in `show`. time = self.df[self.time_col].values out = np.zeros((len(time), len(peak_df))) - iter = 0 - for _, _v in self._peak_props.items(): - for _, v in _v.items(): - params = [v["amplitude"], v["retention_time"], - v["scale"], v["alpha"]] - out[:, iter] = self._compute_skewnorm(time, *params) - iter += 1 + for i, (_, row) in enumerate(peak_df.iterrows()): + params = [row["amplitude"], row["retention_time"], + row["scale"], row["skew"]] + out[:, i] = self._compute_skewnorm(time, *params) self.unmixed_chromatograms = np.round(out, decimals=precision) if return_peaks: return peak_df diff --git a/tests/test_chromatogram.py b/tests/test_chromatogram.py index b8e8061..5f09078 100644 --- a/tests/test_chromatogram.py +++ b/tests/test_chromatogram.py @@ -539,6 +539,142 @@ def test_generic_param_bounding(): assert True +def _skewnorm_signal(t, params): + """Build a signal as a sum of amplitude-weighted (skew)normal peaks.""" + sig = np.zeros_like(t, dtype=float) + for amp, loc, scale, alpha in params: + sig += amp * scipy.stats.skewnorm(alpha, loc, scale).pdf(t) + return sig + + +def test_known_peaks_nonzero_start_time(): + """ + Regression test for bug E: enforced (`known_peaks`) locations are mapped to + array indices relative to the chromatogram's start time, not by assuming the + time axis begins at t = 0. On a chromatogram whose time starts at t != 0, a + shallow peak that is not auto-detected must still be enforced at the correct + retention time. + """ + t = np.arange(10, 30, 0.01) + # A tall auto-detected peak and a shallow peak that prominence filtering skips. + sig = _skewnorm_signal(t, [(1000, 16.0, 0.3, 0), (80, 22.0, 0.3, 0)]) + df = pd.DataFrame({'time': t, 'signal': sig}) + chrom = hplc.quant.Chromatogram(df) + peaks = chrom.fit_peaks(known_peaks={22.0: {'width': 1}}, + prominence=0.5, correct_baseline=False, verbose=False) + # The enforced peak is placed at ~22 (would land out-of-range on the t0=0 bug). + assert np.any(np.abs(peaks['retention_time'].values - 22.0) < 0.5) + assert np.any(np.abs(peaks['retention_time'].values - 16.0) < 0.5) + + +def test_unmixed_columns_match_peak_id(): + """ + Regression test for bug F: `unmixed_chromatograms[:, peak_id - 1]` must hold + the trace for that `peak_id`. When detection order differs from + retention-time order (here forced by enforcing an early peak in a window that + also contains a later auto-detected peak), the columns were previously left + in detection order and disagreed with the sorted `peak_id`. + """ + t = np.arange(0, 30, 0.01) + # Auto-detected tall peak at 16.0; shallow enforced peak earlier at 15.0. + sig = _skewnorm_signal(t, [(1000, 16.0, 0.2, 0), (90, 15.0, 0.2, 0)]) + df = pd.DataFrame({'time': t, 'signal': sig}) + chrom = hplc.quant.Chromatogram(df) + peaks = chrom.fit_peaks(known_peaks={15.0: {'width': 1}}, prominence=0.5, + correct_baseline=False, buffer=200, verbose=False) + + time = chrom.df['time'].values + for _, row in peaks.iterrows(): + col = chrom.unmixed_chromatograms[:, int(row['peak_id']) - 1] + t_at_max = time[np.argmax(col)] + assert np.abs(t_at_max - row['retention_time']) < 0.5, ( + f"column for peak_id {row['peak_id']} peaks at {t_at_max}, " + f"expected near {row['retention_time']}") + + +def test_infeasible_location_guess_does_not_crash(): + """ + Regression test for bug H: the location initial guess is the peak time + *rounded* to `_time_precision`, while its bounds are the window's *raw* + (unrounded) time range. When `_time_precision = |ceil(log10(dt))|` is coarser + than the actual sample spacing, rounding can push the guess just past the + window's max, which scipy rejects as "`x0` is infeasible". The clamp must pull + the guess back inside its bounds so fitting proceeds. + + This is purely a function of `dt` vs `_time_precision`, so it is reproduced + synthetically rather than from recorded data: a step of `dt = 0.011` gives + `_time_precision = |ceil(log10(0.011))| = 1`, so a peak near the end of the + record (here ~2.985) rounds *up* to 3.0 — past its window's raw max (~2.98). + On the unfixed code this raises before any optimizer iteration; with the clamp + the fit completes. (Verified: removing the clamp makes this signal raise + "Initial guess for 'location' ... lies outside its bounds".) + """ + t = np.arange(0, 3.0, 0.011) + sig = _skewnorm_signal(t, [(6000, 1.0, 0.03, 0), (6000, 2.985, 0.03, 0)]) + df = pd.DataFrame({'time': t, 'signal': sig}) + assert int(np.abs(np.ceil(np.log10(np.mean(np.diff(t)))))) == 1, ( + 'fixture must keep `_time_precision == 1` for the rounding excursion') + + chrom = hplc.quant.Chromatogram(df) + try: + peaks = chrom.fit_peaks(correct_baseline=False, verbose=False, + max_iter=5000) + except ValueError as e: + # The clamp removes the bounds/feasibility crash entirely; if a ValueError + # still surfaces it must not be the infeasible-bounds one under test. + assert 'infeasible' not in str(e).lower() + assert 'outside its bounds' not in str(e).lower() + assert 'lower bound' not in str(e).lower() + raise + # The fit proceeds past the bounds check and recovers both peaks. + assert len(peaks) == 2 + + +def test_peak_adjacent_to_start_assigns_interpeak(): + """ + Smoke test for bug G: a chromatogram with a peak right at the start and a + long trailing background must run end-to-end and assign interpeak windows + over the background. The G fix corrects background-window splitting when a + gap lands at the first sample (`split_inds[0] == 0`). That exact branch is + hard to trigger deterministically without depending on `scipy.peak_widths` + internals, so this guards the surrounding behavior rather than isolating the + branch. + """ + t = np.arange(0, 30, 0.01) + sig = _skewnorm_signal(t, [(400, 0.4, 0.15, 0), (600, 18.0, 0.4, 0)]) + df = pd.DataFrame({'time': t, 'signal': sig}) + chrom = hplc.quant.Chromatogram(df) + chrom.fit_peaks(correct_baseline=False, verbose=False, max_iter=5000) + scores = chrom.assess_fit(verbose=False) + # The trailing background between the two peaks is captured as interpeak. + assert (chrom.window_df['window_type'] == 'interpeak').any() + assert (scores['window_type'] == 'interpeak').any() + + +def test_degenerate_bounds_raise_clear_error(): + """ + The pre-fit bounds validation (the safety net added for bug H) must surface + degenerate or infeasible parameter bounds as a clear, actionable ValueError + instead of scipy's opaque message. Issue #22 is exactly the degenerate + amplitude-bound case. + """ + t = np.arange(0, 20, 0.01) + sig = _skewnorm_signal(t, [(1000, 10.0, 0.3, 0)]) + df = pd.DataFrame({'time': t, 'signal': sig}) + + # Degenerate amplitude bound [0, 0]: lower bound not strictly less than upper. + chrom = hplc.quant.Chromatogram(df) + with pytest.raises(ValueError, match='Invalid bounds'): + chrom.fit_peaks(param_bounds={'amplitude': [0, 0]}, + correct_baseline=False, verbose=False) + + # Amplitude bound [2x, 3x] the peak value excludes the initial guess (1x). + chrom = hplc.quant.Chromatogram(df) + with pytest.raises(ValueError, match='lies outside its bounds'): + chrom.fit_peaks(param_bounds={'amplitude': [2, 3]}, + correct_baseline=False, verbose=False) + + def test_multipeak_param_bounds_validated_per_peak(): """ Regression test for bug B: when a global `param_bounds` is applied to a