// SPDX-License-Identifier: MIT
// SPDX-FileCopyrightText: 2026 Martino Ferrari

#include "measure.hpp"

#include <algorithm>
#include <cmath>

namespace udpscope {

namespace {

/// One crossing *into* side `e` at time `t`, which is also what closes the interval the walk was
/// in. Shared by the raw walk and the bucket walk so the two cannot disagree about what an edge
/// records -- an interval's width is the distance from the edge that entered it to the edge that
/// left it, which is the pairing the pulse duration averages.
void enter_side(MeasureLevel& acc, std::size_t e, double t) {
  const std::size_t other = 1 - e;
  if (acc.enters[e] == 0) acc.first_enter[e] = t;
  acc.last_enter[e] = t;
  ++acc.enters[e];
  if (acc.open[other]) {
    acc.width_sum[other] += t - acc.open_at[other];
    ++acc.widths[other];
    acc.open[other] = false;
  }
  acc.open_at[e] = t;
  acc.open[e] = true;
}

/// The hysteresis band's half-width, as a share of the window's peak-to-peak. Measured rather than
/// chosen: 5 % reads the mock's 0.1 interferer correctly and fails at 0.2 (`period=0.075` for a
/// true 0.2), while 10 % holds to an interferer amplitude of 0.3 -- a third of the waveform -- and
/// to white noise at sigma=0.05. Above that the band starts eating a real waveform's own extremes
/// on anything less square than a sine.
constexpr double kHysteresisFrac = 0.10;

}  // namespace

void measure_set_level(MeasureLevel& acc, const MeasureStats& s, double pct) {
  const double span = s.hi - s.lo;
  acc.level = s.lo + pct * span;
  // Half the room to the nearer extreme is the clamp. A threshold at 90 % has a tenth of the span
  // above it, so a fixed 10 % band would put `level + h` *above every sample*: the walk then never
  // decides the upper side and the measurement reads nothing at all, which is strictly worse than
  // a noisy reading. Measured on a 3.3 V logic signal with ringing at thr=90 %: `enters=0` with a
  // fixed band, `enters=4 period=0.200000` with the clamp.
  const double room = std::min(s.hi - acc.level, acc.level - s.lo);
  acc.hysteresis = std::max(0.0, std::min(kHysteresisFrac * span, 0.5 * room));
}

void measure_scan(const MeasureChunk& chunk, MeasureStats& s) {
  // The chunk's own extremes first, in a loop that touches nothing but the values: *where* an
  // extreme is costs an index, and the caller turns the two that survive into stamps once at the
  // end. The sum is accumulated per chunk and folded in once, which is also kinder numerically
  // than one running total over a few million terms.
  double sum = 0.0;
  std::size_t count = 0;
  double lo = 0.0;
  double hi = 0.0;
  std::size_t at_lo = 0;
  std::size_t at_hi = 0;
  bool any = false;
  for (std::size_t i = 0; i < chunk.count; ++i) {
    const double v = static_cast<double>(chunk.value[i]);
    if (!std::isfinite(v)) continue;
    if (!any) {
      lo = hi = v;
      at_lo = at_hi = i;
      any = true;
    } else {
      if (v < lo) {
        lo = v;
        at_lo = i;
      }
      if (v > hi) {
        hi = v;
        at_hi = i;
      }
    }
    sum += v;
    ++count;
  }
  if (!any) return;
  s.sum += sum;
  s.count += count;
  if (!s.valid) {
    s.lo = lo;
    s.hi = hi;
    s.lo_i = chunk.offset + at_lo;
    s.hi_i = chunk.offset + at_hi;
    s.lo_search = 0;
    s.hi_search = 0;
    s.valid = true;
    return;
  }
  // Strictly better only: the *first* sample to reach an extreme is the one a reader's eye lands
  // on, so a later chunk that merely ties it does not move the mark.
  if (lo < s.lo) {
    s.lo = lo;
    s.lo_i = chunk.offset + at_lo;
    s.lo_search = 0;
  }
  if (hi > s.hi) {
    s.hi = hi;
    s.hi_i = chunk.offset + at_hi;
    s.hi_search = 0;
  }
}

void measure_scan_buckets(const MeasureBucket* buckets, std::size_t count, MeasureStats& s) {
  double sum = 0.0;
  std::size_t finite = 0;
  for (std::size_t b = 0; b < count; ++b) {
    const MeasureBucket& k = buckets[b];
    // A bucket of nothing (every sample non-finite) is stored as an empty aggregate, and the raw
    // walk skips exactly those samples -- so it contributes nothing here either.
    if (k.count == 0) continue;
    sum += k.sum;
    finite += k.count;
    if (!s.valid) {
      s.lo = k.lo;
      s.hi = k.hi;
      s.lo_i = s.hi_i = k.offset;
      s.lo_search = s.hi_search = k.span;
      s.valid = true;
      continue;
    }
    if (k.lo < s.lo) {
      s.lo = k.lo;
      s.lo_i = k.offset;
      s.lo_search = k.span;
    }
    if (k.hi > s.hi) {
      s.hi = k.hi;
      s.hi_i = k.offset;
      s.hi_search = k.span;
    }
  }
  s.sum += sum;
  s.count += finite;
}

void measure_accumulate(const MeasureChunk& chunk,
                        const std::function<double(std::size_t)>& time_of, MeasureLevel* levels,
                        std::size_t level_count) {
  for (std::size_t k = 0; k < level_count; ++k) {
    MeasureLevel& acc = levels[k];
    for (std::size_t i = 0; i < chunk.count; ++i) {
      const double v = static_cast<double>(chunk.value[i]);
      if (!std::isfinite(v)) continue;
      // The *sums* split at the bare level: `Top` and `Bottom` are defined as the mean above and
      // below the threshold, and a band here would move a documented exact reading.
      const int s = v > acc.level ? 1 : 0;
      acc.sum[static_cast<std::size_t>(s)] += v;
      ++acc.count[static_cast<std::size_t>(s)];
      // The *edge* decision reads the band instead. A sample inside it is on neither side, so it
      // decides no crossing and the walk keeps the side it last decided -- which is what stops a
      // signal that merely wobbles about the level from counting every wobble as an edge. With
      // `hysteresis == 0` this is the single-level test, `v == level` skip included.
      int d;
      if (v > acc.level + acc.hysteresis) {
        d = 1;
      } else if (v < acc.level - acc.hysteresis) {
        d = 0;
      } else {
        continue;
      }
      if (acc.have_side && d != (acc.side ? 1 : 0)) {
        enter_side(acc, static_cast<std::size_t>(d), time_of(chunk.offset + i));
      }
      acc.side = d != 0;
      acc.have_side = true;
    }
  }
}

bool measure_accumulate_bucket(const MeasureBucket& k,
                               const std::function<double(std::size_t)>& time_of,
                               MeasureLevel& acc) {
  if (k.count == 0) return true;
  // A non-finite sample in the bucket would make the edge's position a guess: the raw walk takes
  // the first *finite* sample of the new side, and an aggregate cannot say which that is.
  if (k.count != k.span) return false;
  // The *sums* need only the bare level: strictly on one side of it, so every sample decides the
  // same bin and none is *at* it, which is what makes the fold exact rather than an approximation.
  const bool above = k.lo > acc.level;
  const bool below = k.hi < acc.level;
  if (!above && !below) return false;
  const std::size_t e = above ? 1u : 0u;
  // The band is deliberately *not* folded into the test above. Tightening the whole condition to
  // `level +/- h` refuses every bucket lying inside the band -- a tenth of the amplitude either side
  // of the mid level is where a waveform moves fastest, so those are exactly the buckets a crossing
  // passes through -- which hands the raw walk a share of the window rather than a bucket or two per
  // crossing, and breaks the "`raw` stays at a fraction of a percent of `covered`" invariant.
  // Instead: a bucket on one side of the level can always be *binned*, and only a *side change*
  // whose position is inside the bucket is refused.
  const double top = acc.level + acc.hysteresis;
  const double bot = acc.level - acc.hysteresis;
  // Does the bucket reach past the band at all? One wholly inside it decides no side, exactly as a
  // raw sample inside it decides none, so the walk state is left alone.
  const bool reaches = above ? (k.hi > top) : (k.lo < bot);
  const bool change = reaches && acc.have_side && e != (acc.side ? 1u : 0u);
  // A bucket that *enters* a side while straddling its own threshold holds the crossing somewhere
  // inside itself, and only the raw walk can say where. One wholly past the threshold holds it at
  // the bucket's first sample, exactly, which is where the raw walk would have found it.
  const bool fully_outside = above ? (k.lo > top) : (k.hi < bot);
  // Nothing may be folded before this point: a refusal hands the bucket to `measure_accumulate()`,
  // which adds the same samples again. Sums added above the refusal double-counted every opened
  // bucket, which read as `top` 0.614304 on a unit sine against the 2/pi it must be.
  if (change && !fully_outside) return false;
  acc.sum[e] += k.sum;
  acc.count[e] += k.count;
  if (change) enter_side(acc, e, time_of(k.offset));
  if (reaches) {
    acc.side = e != 0;
    acc.have_side = true;
  }
  return true;
}

bool measure_level_pct(MeasureType type, const MeasureOptions& options, double& pct) {
  switch (type) {
    case MeasureType::Top:
    case MeasureType::Bottom:
      pct = std::clamp(options.threshold_pct, 0.0, 100.0) / 100.0;
      return true;
    case MeasureType::PulseDuration:
    case MeasureType::Period:
    case MeasureType::Frequency:
      // The *mid* level, for all three. No option here on purpose: a period is the spacing of
      // one crossing per cycle, and a threshold near either extreme counts noise as crossings.
      pct = 0.5;
      return true;
    default:
      return false;
  }
}


const char* measure_type_name(MeasureType type) {
  switch (type) {
    case MeasureType::Average: return "average";
    case MeasureType::Min: return "min";
    case MeasureType::Max: return "max";
    case MeasureType::PeakToPeak: return "peak to peak";
    case MeasureType::Amplitude: return "amplitude";
    case MeasureType::Top: return "top";
    case MeasureType::Bottom: return "bottom";
    case MeasureType::PulseDuration: return "pulse duration";
    case MeasureType::Period: return "period";
    case MeasureType::Frequency: return "frequency";
  }
  return "?";
}

int measure_type_count() { return 10; }

MeasureType measure_type_at(int index) {
  switch (index) {
    case 0: return MeasureType::Average;
    case 1: return MeasureType::Top;
    case 2: return MeasureType::Bottom;
    case 3: return MeasureType::Max;
    case 4: return MeasureType::Min;
    case 5: return MeasureType::PeakToPeak;
    case 6: return MeasureType::Amplitude;
    case 7: return MeasureType::PulseDuration;
    case 8: return MeasureType::Frequency;
    default: return MeasureType::Period;
  }
}

bool measure_uses_threshold(MeasureType type) {
  return type == MeasureType::Top || type == MeasureType::Bottom ||
         type == MeasureType::PulseDuration;
}

bool measure_uses_side(MeasureType type) { return type == MeasureType::PulseDuration; }

bool measure_can_show_on_plot(MeasureType) { return true; }

MeasureResult compute_measure(const MeasureStats& b, const MeasureLevel* level, MeasureType type,
                              const MeasureOptions& options) {
  MeasureResult out;
  if (!b.valid) return out;

  switch (type) {
    case MeasureType::Average: {
      if (b.count == 0) return out;
      out.primary = b.sum / static_cast<double>(b.count);
      // A mean is a level like any other, so it can be shown as one ruler.
      out.marks.push_back({true, out.primary});
      out.valid = true;
      return out;
    }
    case MeasureType::Min:
    case MeasureType::Max: {
      // A ruler at the extreme, and a marker at the time it happened: "the minimum" is a level
      // and a place. The place comes from the shared scan, which takes the *first* sample to
      // reach the extreme -- the one a reader's eye lands on when the extreme is a plateau.
      const bool want_max = type == MeasureType::Max;
      out.primary = want_max ? b.hi : b.lo;
      out.marks.push_back({true, out.primary});
      out.marks.push_back({false, want_max ? b.at_hi : b.at_lo});
      out.valid = true;
      return out;
    }
    case MeasureType::PeakToPeak:
    case MeasureType::Amplitude: {
      // Both come straight out of the window's own extremes, which every type already needs:
      // the two rulers are the same pair, and the numbers differ by a factor of two.
      out.primary = type == MeasureType::PeakToPeak ? b.hi - b.lo : 0.5 * (b.hi - b.lo);
      // A ruler at each extreme: the *span* is the measurement, and two levels are how a
      // plot states a span without a band between them.
      out.marks.push_back({true, b.lo});
      out.marks.push_back({true, b.hi});
      out.valid = true;
      return out;
    }
    case MeasureType::Top:
    case MeasureType::Bottom: {
      // The mean of one side of the level, straight out of the second pass' running sums. A flat
      // window has no range at all, so every sample lands on the "below" side -- which is why
      // `top` on a signal that never moves reads `--` rather than that one value.
      if (level == nullptr) return out;
      const std::size_t s = type == MeasureType::Top ? 1 : 0;
      if (level->count[s] == 0) return out;
      out.primary = level->sum[s] / static_cast<double>(level->count[s]);
      out.marks.push_back({true, out.primary});
      out.valid = true;
      return out;
    }
    case MeasureType::PulseDuration:
    case MeasureType::Period:
    case MeasureType::Frequency: {
      if (level == nullptr) return out;
      // The pulse is one side of the mid level: its widths run from the edge that enters it to
      // the edge that leaves it, and a period is the spacing of the edges that *enter* it -- one
      // per cycle, whatever else the waveform does in between.
      const std::size_t s = options.above ? 1 : 0;
      const double enter_span = level->last_enter[s] - level->first_enter[s];
      const double enter_period = level->enters[s] >= 2
                                      ? enter_span / static_cast<double>(level->enters[s] - 1)
                                      : 0.0;
      if (type == MeasureType::PulseDuration) {
        if (level->widths[s] == 0) return out;
        out.primary = level->width_sum[s] / static_cast<double>(level->widths[s]);
        // Duty as width/period, not as pulse-time over the window: a window that starts or
        // ends mid-pulse holds partial cycles, and counting only the complete ones against
        // the whole span would report a duty that depends on where the window happened to
        // begin. The period is the spacing of the pulse's own enter-edges, and it needs two of
        // them -- a window holding one pulse reports the width alone.
        out.secondary = enter_period > 0.0 ? out.primary / enter_period : 0.0;
        out.has_secondary = enter_period > 0.0;
        out.marks.push_back({false, level->first_enter[s]});
        out.marks.push_back({false, level->first_enter[s] + out.primary});
        out.valid = true;
        return out;
      }
      if (!(enter_period > 0.0)) return out;
      out.primary = type == MeasureType::Period ? enter_period : 1.0 / enter_period;
      out.marks.push_back({false, level->first_enter[s]});
      out.marks.push_back({false, level->first_enter[s] + enter_period});
      out.valid = true;
      return out;
    }
  }
  return out;
}


}  // namespace udpscope