// analysis_rates.mbt — population-level firing rate helpers (v0.41.3).
//
// Complements `vecplot.mbt::Monitor::firing_rate` (which operates on
// the rising-edge buffer) by providing direct spike-train utilities
// that accept arrays of spike times — useful for offline analysis
// and for cross-validation against Julia.
//
// Conventions:
//   - spike times are in milliseconds (Float)
//   - rates are returned in Hz (spikes per second)
//   - time windows are inclusive on both ends
//
// API:
//   neuron_firing_rate(spike_times, t_start, t_end) → Float (Hz)
//   population_firing_rate(per_neuron_spikes, t_start, t_end) → Float (Hz)
//     (mean over active neurons)
//   population_total_rate(per_neuron_spikes, t_start, t_end) → Float (Hz)
//     (sum over all neurons)
//   firing_rate_dynamics(per_neuron_spikes, window, dt) → Array[Float]
//     (sliding-window mean rate at each timestep)
//   fano_factor(per_neuron_counts) → Float
//     (variance / mean over per-neuron spike counts)

///|
/// Count spikes in `spike_times` that fall within `[t_start, t_end]`.
/// Both bounds are inclusive.
fn count_in_window(spike_times : Array[Float], t_start : Float, t_end : Float) -> Int {
  let mut count = 0
  for i in 0..= t_start && t <= t_end {
      count = count + 1
    }
  }
  count
}

///|
/// Firing rate (Hz) of a single neuron over `[t_start, t_end]`.
/// Returns 0.0F if the window duration is zero.
pub fn neuron_firing_rate(
  spike_times : Array[Float],
  t_start : Float,
  t_end : Float,
) -> Float {
  let duration_ms = t_end - t_start
  if duration_ms <= 0.0F {
    return 0.0F
  }
  let count = count_in_window(spike_times, t_start, t_end)
  Float::from_int(count) * 1000.0F / duration_ms
}

///|
/// Mean firing rate (Hz) across a population (one element per neuron).
/// Silent neurons contribute 0 Hz. Window is `[t_start, t_end]`.
pub fn population_firing_rate(
  per_neuron_spikes : Array[Array[Float]],
  t_start : Float,
  t_end : Float,
) -> Float {
  let n_pop = per_neuron_spikes.length()
  if n_pop == 0 {
    return 0.0F
  }
  let mut total = 0.0F
  for i in 0.. Float {
  let n_pop = per_neuron_spikes.length()
  if n_pop == 0 {
    return 0.0F
  }
  let mut total = 0.0F
  for i in 0.. 0.0F {
      let count = count_in_window(per_neuron_spikes[i], t_start, t_end)
      total = total + Float::from_int(count) * 1000.0F / duration_ms
    }
  }
  total
}

///|
/// Sliding-window mean firing rate (Hz). Returns an array of length
/// `n_steps` where each entry is the mean rate over the window
/// `[i*dt, i*dt + window]`. `window` is in ms.
pub fn firing_rate_dynamics(
  per_neuron_spikes : Array[Array[Float]],
  window : Float,
  dt : Float,
) -> Array[Float] {
  if dt <= 0.0F {
    return []
  }
  // Determine total duration from the last spike across all neurons.
  let mut t_max = 0.0F
  for i in 0.. 0 {
      let last = spikes[n - 1]
      if last > t_max {
        t_max = last
      }
    }
  }
  let n_steps = Float::to_int(t_max / dt + 0.999F)
  let out : Array[Float] = Array::make(n_steps, 0.0F)
  let n_pop = per_neuron_spikes.length()
  for s in 0.. Float {
  let n = per_neuron_counts.length()
  if n <= 1 {
    return 0.0F
  }
  let n_f = Float::from_int(n)
  let mut sum : Float = 0.0F
  for i in 0.. Array[Float] {
  let n_pop = per_neuron_spikes.length()
  let out : Array[Float] = Array::make(n_pop, 0.0F)
  for i in 0.. Float {
  let n_pop = per_neuron_spikes.length()
  if n_pop == 0 {
    return 0.0F
  }
  let rates = per_neuron_firing_rates(per_neuron_spikes, t_start, t_end)
  let n_f = Float::from_int(n_pop)
  let mut sum : Float = 0.0F
  for i in 0..