// 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.
// The hysteresis band and its rationale moved to `MeasureOptions::noise_pct`: it is a per-measurement
// setting now, and a second copy of the number here would be a second thing to keep in step.
constexpr double kPi = 3.14159265358979323846;

/// The flat level's band, as a half-width in shares of the window's peak-to-peak. 5 % is *under*
/// the overshoot a band-limited edge rings to -- Gibbs, about 9 % of the swing at the harmonic
/// counts the mock's `square` uses -- so the ringing's excursion falls outside the densest band
/// while a plateau's own ripple, a few percent, stays inside it. Measured on the mock's 30-harmonic
/// square at 1.3 MS/s: the flat top reads 1.000000 with this band, 1.003573 at 3 %, 1.003495 at
/// 10 % -- the wider bands reach the ringing again. A share and not an absolute width: the band is
/// the same quantity `threshold_pct` and the hysteresis are measured in.
constexpr double kFlatBandShare = 0.05;

/// Which bin of a side is the level: the *outermost* one holding at least this share of the
/// densest bin's buckets. One rule for both kinds of waveform. Where a plateau exists it is by
/// far the densest band and the rule cannot reach past it. Where none exists -- a sine, a ramp --
/// the signal dwells in every band it visits, the whole side qualifies, and the outermost is the
/// honest answer: the level is then the top of what the window reaches, near the peak, rather
/// than a bin of an arbitrary mid-range density that a shifted window would move.
constexpr double kFlatTieShare = 0.5;

}  // namespace

void measure_set_level(MeasureLevel& acc, const MeasureStats& s, double pct, double noise_frac) {
  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);
  const double want = noise_frac > 0.0 ? noise_frac * span : 0.0;
  acc.hysteresis = std::max(0.0, std::min(want, 0.5 * room));
}

void measure_filter_set(MeasureFilter& f, double fc, double dt) {
  f.reset();
  f.a = 1.0;
  if (fc > 0.0 && dt > 0.0) {
    const double nyquist = 0.5 / dt;
    const double use = fc < nyquist * 0.99 ? fc : nyquist * 0.99;
    f.a = 1.0 - std::exp(-2.0 * kPi * use * dt);
  }
}

void measure_scan(const MeasureChunk& chunk, MeasureStats& s, MeasureFilter& filter) {
  // 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 = filter.step(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,
                          bool flat_hist) {
  double sum = 0.0;
  std::size_t finite = 0;
  // The *buckets'* own range, which is not the window's: the window's partial edge buckets are
  // read raw and are not in this list. Only the non-empty buckets: an all-non-finite bucket is
  // stored as an empty aggregate and its `lo`/`hi` say nothing.
  double blo = 0.0;
  double bhi = 0.0;
  bool any = false;
  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;
    if (!any) {
      blo = k.lo;
      bhi = k.hi;
      any = true;
    } else {
      if (k.lo < blo) blo = k.lo;
      if (k.hi > bhi) bhi = k.hi;
    }
    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;
  if (!flat_hist || !any || !(bhi > blo)) return;
  // One observation per bucket, at its midpoint, carrying that bucket's *exact* sum and finite
  // count -- so the level the caller finally reports is a mean of real samples and not a bin
  // centre. Weighting by the bucket (and not by its 64 samples) is what keeps the flat regions
  // and the ringing on the same footing: the question is how much of the *waveform* dwells in a
  // band, and a bucket is a span of it. A range that is degenerate (every bucket the same value)
  // has no band to pick and leaves `folded` at 0, which is the caller's "no histogram".
  s.flat.lo = blo;
  s.flat.hi = bhi;
  const double scale = static_cast<double>(MeasureFlatHist::kBins) / (bhi - blo);
  for (std::size_t b = 0; b < count; ++b) {
    const MeasureBucket& k = buckets[b];
    if (k.count == 0) continue;
    const double mid = 0.5 * (k.lo + k.hi);
    std::size_t bin = mid <= blo ? 0u : static_cast<std::size_t>((mid - blo) * scale);
    if (bin >= MeasureFlatHist::kBins) bin = MeasureFlatHist::kBins - 1;
    s.flat.sum[bin] += k.sum;
    s.flat.count[bin] += k.count;
    ++s.flat.buckets[bin];
    ++s.flat.folded;
  }
}

void measure_accumulate(const MeasureChunk& chunk,
                        const std::function<double(std::size_t)>& time_of, MeasureLevel* levels,
                        std::size_t level_count, MeasureFilter& filter) {
  // Samples outer, levels inner -- the reverse of what reads more naturally, and deliberately. The filter
  // holds one state for the whole walk, so a levels-outer loop would step it once per level per sample and
  // every level after the first would read a signal filtered `k` times over. Each level's own state is
  // independent of the others, so the swap costs nothing else.
  for (std::size_t i = 0; i < chunk.count; ++i) {
    const double v = filter.step(static_cast<double>(chunk.value[i]));
    if (!std::isfinite(v)) continue;
    for (std::size_t k = 0; k < level_count; ++k) {
      MeasureLevel& acc = levels[k];
      // 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;
    }
  }
}


MeasureFlatLevels measure_flat_levels(const MeasureStats& b, double pct) {
  MeasureFlatLevels out;
  const MeasureFlatHist& h = b.flat;
  const double span = b.hi - b.lo;
  if (h.folded == 0 || !(h.hi > h.lo) || !(span > 0.0)) return out;
  // The split, sized exactly as `measure_set_level()` sizes the level the crossing types read, so
  // "the top" and "the level a period is counted about" are placed by one rule.
  const double level = b.lo + std::clamp(pct, 0.0, 1.0) * span;
  const double width = (h.hi - h.lo) / static_cast<double>(MeasureFlatHist::kBins);
  const double band = kFlatBandShare * span;
  for (int upper = 1; upper >= 0; --upper) {
    const auto in_side = [&](std::size_t i) {
      const double centre = h.lo + (static_cast<double>(i) + 0.5) * width;
      return (centre > level) == (upper != 0);
    };
    // One scan for the densest bin of this side, one for the outermost of the equally dense
    // ones, one to sum the band about it: 192 iterations of an array of 64, per flat level.
    std::size_t densest = 0;
    for (std::size_t i = 0; i < MeasureFlatHist::kBins; ++i) {
      if (in_side(i) && h.buckets[i] > densest) densest = h.buckets[i];
    }
    // A side of the split with no bucket in it has no level: `top` on a window that never reaches
    // above the split reads `--`, which is exactly what it read before.
    if (densest == 0) continue;
    const double floor_buckets = kFlatTieShare * static_cast<double>(densest);
    std::size_t anchor = 0;
    bool have_anchor = false;
    for (std::size_t i = 0; i < MeasureFlatHist::kBins; ++i) {
      if (!in_side(i) || static_cast<double>(h.buckets[i]) < floor_buckets) continue;
      // The upper side wants the *last* qualifier and the lower one the first, so a bottom-up scan
      // is "always overwrite" for one and "only ever write once" for the other.
      if (upper != 0 || !have_anchor) anchor = i;
      have_anchor = true;
    }
    if (!have_anchor) continue;
    const double centre = h.lo + (static_cast<double>(anchor) + 0.5) * width;
    double sum = 0.0;
    std::size_t count = 0;
    for (std::size_t i = 0; i < MeasureFlatHist::kBins; ++i) {
      const double offset = h.lo + (static_cast<double>(i) + 0.5) * width - centre;
      if (offset < -band || offset > band) continue;
      sum += h.sum[i];
      count += h.count[i];
    }
    if (count == 0) continue;
    if (upper != 0) {
      out.top = sum / static_cast<double>(count);
      out.top_ok = true;
    } else {
      out.bottom = sum / static_cast<double>(count);
      out.bottom_ok = true;
    }
  }
  return out;
}

int measure_level_pcts(MeasureType type, const MeasureOptions& options, double pct[2]) {
  const auto share = [](double v) { return std::clamp(v, 0.0, 100.0) / 100.0; };
  switch (type) {
    case MeasureType::Top:
    case MeasureType::Bottom:
    case MeasureType::Above:
    case MeasureType::Below:
      pct[0] = share(options.threshold_pct);
      return 1;
    case MeasureType::Between: {
      // The band's two edges, clamped and *ordered*: the shares name a lower and an upper edge, so
      // "between 80 % and 20 %" is the band between 20 % and 80 %. Normalizing rather than refusing
      // is the honest reading of a band -- the set of samples does not depend on which edge the user
      // typed first -- and the two accumulators the pair is answered from are keyed on these ordered
      // values, so the same band is one walk however it was written.
      const double a = share(options.threshold_pct);
      const double b = share(options.threshold_pct2);
      pct[0] = std::min(a, b);
      pct[1] = std::max(a, b);
      return 2;
    }
    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] = 0.5;
      return 1;
    default:
      return 0;
  }
}

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";
    case MeasureType::Above: return "above";
    case MeasureType::Below: return "below";
    case MeasureType::Between: return "between";
  }
  return "?";
}

int measure_type_count() { return 13; }

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;
    case 9: return MeasureType::Period;
    // Appended rather than filed beside `Top`/`Bottom`, where they read best in the dropdown:
    // `layout_store` writes *this* index into the session file, so renumbering would turn a saved
    // `max` into an `above` on the next start -- silently, and only for users who had one.
    case 10: return MeasureType::Above;
    case 11: return MeasureType::Below;
    default: return MeasureType::Between;
  }
}

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

bool measure_uses_band(MeasureType type) { return type == MeasureType::Between; }

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

bool measure_uses_noise(MeasureType type) {
  return type == MeasureType::PulseDuration || type == MeasureType::Period ||
         type == MeasureType::Frequency;
}

bool measure_can_show_on_plot(MeasureType) { return true; }

MeasureResult compute_measure(const MeasureStats& b, const MeasureLevel* const levels[2],
                              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: {
      // Straight out of the window's own extremes, which every type already needs: this one *is*
      // the extremes, ripple and overshoot included -- that is the difference between it and
      // `Amplitude`, which reads the two levels the signal dwells at instead.
      out.primary = 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::Amplitude: {
      // The delta between the two flat levels, so a ringing edge's overshoot is not part of it.
      // The split is the mid level, fixed rather than an option -- the dialog offers no threshold
      // for this type -- and the reading is insensitive to it anyway: the band a level is read from
      // is a *value* band, so any split between the two plateaus picks the same two.
      const MeasureFlatLevels flat = measure_flat_levels(b, 0.5);
      if (flat.top_ok && flat.bottom_ok) {
        out.primary = flat.top - flat.bottom;
        out.marks.push_back({true, flat.top});
        out.marks.push_back({true, flat.bottom});
        out.valid = true;
        return out;
      }
      // No density estimate to separate a plateau from its ringing: a window shorter than a
      // bucket, or an HF-rejected walk, which cannot use the pyramid at all. The delta then
      // degenerates to the window's own peak to peak -- the same number when there is nothing to
      // exclude, and the only honest one when the ripple cannot be told from the level.
      out.primary = b.hi - b.lo;
      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 flat level of one side of the split: the densest band there, whose samples' exact mean
      // is the value. The band is a *value* band, so the reading does not follow the window's
      // alignment, and it excludes both a plateau's ripple and an edge's overshoot -- which is the
      // whole point, because those are exactly what a mean over a side reports instead.
      double pct[2] = {0.0, 0.0};
      if (measure_level_pcts(type, options, pct) == 0) return out;
      const bool want_top = type == MeasureType::Top;
      const MeasureFlatLevels flat = measure_flat_levels(b, pct[0]);
      if (want_top ? flat.top_ok : flat.bottom_ok) {
        out.primary = want_top ? flat.top : flat.bottom;
        out.marks.push_back({true, out.primary});
        out.valid = true;
        return out;
      }
      // No histogram: the mean of one side of the level, straight out of the second pass' running
      // sums. A window the pyramid cannot summarize is one bucket or two of a zoomed-in view, or
      // an HF-rejected walk, and on those the split's own means are the honest reading -- and they
      // are what `Above`/`Below` report on *every* window.
      const MeasureLevel* lv = levels != nullptr ? levels[0] : nullptr;
      const std::size_t s = want_top ? 1 : 0;
      if (lv == nullptr || lv->count[s] == 0) return out;
      out.primary = lv->sum[s] / static_cast<double>(lv->count[s]);
      out.marks.push_back({true, out.primary});
      out.valid = true;
      return out;
    }
    case MeasureType::Above:
    case MeasureType::Below: {
      // The plain mean of the samples on one side of the share -- what `Top`/`Bottom` reported
      // before they became flat levels, and a different question from them: a mean over a side
      // carries the ripple and the edge's overshoot (0.986496 on the rippled square where the
      // plateau is 1), which is exactly why the flat-level types exist beside these.
      const MeasureLevel* lv = levels != nullptr ? levels[0] : nullptr;
      const std::size_t s = type == MeasureType::Above ? 1 : 0;
      if (lv == nullptr || lv->count[s] == 0) return out;
      out.primary = lv->sum[s] / static_cast<double>(lv->count[s]);
      out.marks.push_back({true, out.primary});
      out.valid = true;
      return out;
    }
    case MeasureType::Between: {
      // The band's interior: what lies above the lower edge and *not* above the upper one. Both
      // accumulators split the same samples the same way -- pass two walks one bucket list and one
      // filter for the whole window -- so the band is the *difference* of the two upper sides, and
      // it is exact rather than approximate. It needs no third walk of its own, and the bucket
      // fast path serves it unchanged: a bucket wholly inside the band is strictly above the lower
      // level and strictly below the upper, so both accumulators fold it whole.
      const MeasureLevel* lower = levels != nullptr ? levels[0] : nullptr;
      const MeasureLevel* upper = levels != nullptr ? levels[1] : nullptr;
      if (lower == nullptr || upper == nullptr) return out;
      // Equal counts mean an empty interior: nothing lies in the band, so there is nothing to
      // average -- the pane reads `--`, as it does for any other type with no samples on the side
      // it was asked about. `lower` may not count *fewer*, which would say the two accumulators
      // saw different samples.
      if (lower->count[1] <= upper->count[1]) return out;
      out.primary = (lower->sum[1] - upper->sum[1]) /
                    static_cast<double>(lower->count[1] - upper->count[1]);
      out.marks.push_back({true, out.primary});
      out.valid = true;
      return out;
    }
    case MeasureType::PulseDuration:
    case MeasureType::Period:
    case MeasureType::Frequency: {
      if (levels == nullptr || levels[0] == nullptr) return out;
      const MeasureLevel& level = *levels[0];
      // 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