from dataclasses import dataclass import numpy as np TRIGGER_THRESHOLD = 30.0 PULSES_PER_REVOLUTION = 8 @dataclass(frozen=True) class DetectedCycle: number: int start_offset: int end_offset: int angle: np.ndarray signal: np.ndarray trigger_offsets: tuple[int, ...] def _trigger_runs(trigger: np.ndarray) -> list[tuple[int, int, int]]: """Return start offset, end offset and length for each high trigger run.""" values = np.asarray(trigger, dtype=float) high = np.nan_to_num(values, nan=-np.inf) >= TRIGGER_THRESHOLD changes = np.diff(np.r_[False, high, False].astype(np.int8)) starts = np.flatnonzero(changes == 1) ends = np.flatnonzero(changes == -1) - 1 return [ (int(start), int(end), int(end - start + 1)) for start, end in zip(starts, ends) ] def _zero_marker_threshold(lengths: np.ndarray) -> float: """Adaptively separate ordinary keyphasor pulses from the wider 0-degree marker. The wide marker's extra sensing time varies by machine (measured ratios from ~1.4x to ~1.7x of the ordinary pulse), so a fixed margin misses machines whose zero marker is only slightly wider. Otsu-style thresholding splits the run lengths into two clusters regardless of the exact ratio; a near-single cluster (ratio < 1.15) means no reliable zero marker exists and the returned threshold rejects every run. """ lengths = np.asarray(lengths, dtype=float) if len(lengths) < 3: return float("inf") uniq = np.unique(lengths) if len(uniq) < 2: return float("inf") best_variance = -1.0 best_threshold = float("inf") best_means = (float(lengths.min()), float(lengths.max())) for t in uniq[:-1]: low = lengths[lengths <= t] high = lengths[lengths > t] if len(low) == 0 or len(high) == 0: continue weight_low = len(low) / len(lengths) weight_high = len(high) / len(lengths) mean_low = float(low.mean()) mean_high = float(high.mean()) between = weight_low * weight_high * (mean_low - mean_high) ** 2 if between > best_variance: best_variance = between best_threshold = float(t) best_means = (mean_low, mean_high) if best_means[1] < best_means[0] * 1.15: return float("inf") return best_threshold def detect_cycles(samples: np.ndarray) -> tuple[list[DetectedCycle], dict[str, float | int]]: """Detect complete crankshaft cycles from a wave_sample array. The trigger contract intentionally mirrors showPV: second_value >= 30 is high, the longer high run marks 0 degrees, and eight pulses make one turn. Offsets are array offsets rather than sample_index values so the function also works when a file has a non-zero or sparse sample_index column. """ values = np.asarray(samples, dtype=float) if values.ndim != 2 or values.shape[1] < 3 or not len(values): return [], { "triggerRunCount": 0, "zeroMarkerCount": 0, "completeCycleCount": 0, } runs = _trigger_runs(values[:, 2]) diagnostics: dict[str, float | int] = { "triggerRunCount": len(runs), "zeroMarkerCount": 0, "completeCycleCount": 0, } if len(runs) < PULSES_PER_REVOLUTION + 1: return [], diagnostics lengths = np.asarray([run[2] for run in runs], dtype=float) ordinary_length = float(np.median(lengths)) long_limit = _zero_marker_threshold(lengths) zero_runs = [run for run in runs if run[2] > long_limit] diagnostics.update( ordinaryTriggerLength=ordinary_length, zeroMarkerThreshold=long_limit, zeroMarkerCount=len(zero_runs), ) if len(zero_runs) < 2: return [], diagnostics cycles: list[DetectedCycle] = [] signal = values[:, 1] for cycle_number, (zero_run, next_zero_run) in enumerate( zip(zero_runs, zero_runs[1:]), start=1, ): start = zero_run[0] end = next_zero_run[0] cycle_runs = [run for run in runs if start <= run[0] <= end] if len(cycle_runs) != PULSES_PER_REVOLUTION + 1 or end <= start: continue angle = np.full(end - start, np.nan, dtype=float) for pulse_index in range(PULSES_PER_REVOLUTION): segment_start = cycle_runs[pulse_index][0] segment_end = cycle_runs[pulse_index + 1][0] local_start = segment_start - start local_end = segment_end - start if local_end <= local_start: continue angle[local_start:local_end] = np.linspace( pulse_index * 45.0, (pulse_index + 1) * 45.0, local_end - local_start, endpoint=False, ) cycles.append( DetectedCycle( number=cycle_number, start_offset=start, end_offset=end, angle=angle, signal=signal[start:end].copy(), trigger_offsets=tuple(run[0] for run in cycle_runs[:-1]), ) ) diagnostics["completeCycleCount"] = len(cycles) return cycles, diagnostics def build_angle_vector(sample_count: int, cycles: list[DetectedCycle]) -> np.ndarray: """Build a sparse full-file crank-angle vector (0-180 degrees). Mirrors showPV's display angle: the 0-360 keyphasor angle is folded so that the piston stroke position reads 0 -> 180 -> 0 across a full revolution. """ angle = np.full(sample_count, np.nan, dtype=float) for cycle in cycles: full = cycle.angle angle[cycle.start_offset : cycle.end_offset] = np.where( full <= 180.0, full, 360.0 - full, ) return angle def downsample_indices( values: np.ndarray, target_count: int, required_indices: set[int] | None = None, ) -> np.ndarray: """Min/max downsample while retaining requested period boundaries. Each bucket contributes its first/last/min/max points. Required offsets are added afterwards, so trigger boundaries survive even when the overview has far fewer pixels than the raw waveform. """ values = np.asarray(values, dtype=float) count = len(values) if count == 0: return np.empty(0, dtype=np.int64) target_count = max(int(target_count), 16) required = { int(index) for index in (required_indices or set()) if 0 <= int(index) < count } if count <= target_count: return np.arange(count, dtype=np.int64) bucket_count = max(1, target_count // 4) edges = np.linspace(0, count, bucket_count + 1, dtype=np.int64) selected = set(required) for bucket_index in range(bucket_count): start = int(edges[bucket_index]) end = int(edges[bucket_index + 1]) if end <= start: continue bucket = values[start:end] finite = np.isfinite(bucket) candidates = {start, end - 1} if finite.any(): finite_values = bucket.copy() finite_values[~finite] = np.nan candidates.add(start + int(np.nanargmin(finite_values))) candidates.add(start + int(np.nanargmax(finite_values))) selected.update(candidates) if len(selected) > target_count: required_sorted = sorted(required) remaining = sorted(selected.difference(required)) slots = max(0, target_count - len(required_sorted)) if slots: positions = np.linspace(0, len(remaining) - 1, slots, dtype=np.int64) selected = set(required_sorted).union(remaining[int(position)] for position in positions) else: # Keep every period boundary even if it exceeds the requested budget. selected = set(required_sorted) elif len(selected) < target_count: missing = np.setdiff1d( np.arange(count, dtype=np.int64), np.asarray(sorted(selected), dtype=np.int64), assume_unique=True, ) slots = min(target_count - len(selected), len(missing)) if slots: positions = np.linspace(0, len(missing) - 1, slots, dtype=np.int64) selected.update(int(missing[int(position)]) for position in positions) return np.asarray(sorted(selected), dtype=np.int64)