Compressor/backend/app/algorithms/cycles.py

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 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 = ordinary_length + max(5.0, ordinary_length * 0.45)
    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)