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 _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)