///|
/// Summary of one wavelet scale.
pub struct WaveletScaleSummary {
  scale : Double
  pseudo_frequency : Double
  energy : Double
  maximum_magnitude : Double
  mean_absolute : Double
  active_fraction : Double
}

///|
/// Multi-resolution event features extracted from a scalogram.
pub struct WaveletEventFeatures {
  dominant_scale : Double
  dominant_frequency : Double
  peak_time_index : Int
  peak_magnitude : Double
  entropy : Double
  onset_index : Int
  duration_samples : Int
  coherence : Double
}

///|
/// Convert a scale into a pseudo-frequency for a sample rate.
pub fn wavelet_pseudo_frequency(
  scale : Double,
  sample_rate : Double,
  family : WaveletFamily,
) -> Double {
  let center = match family {
    RickerWavelet => 0.225
    MorletWavelet => 0.796
    GaussianDerivative => 0.159
  }
  if scale <= 0.0 || sample_rate <= 0.0 {
    0.0
  } else {
    center * sample_rate / scale
  }
}

///|
/// Summarize every scale with energy and occupancy measures.
pub fn summarize_wavelet_scales(
  scalogram : Scalogram,
  sample_rate : Double,
  active_threshold : Double,
) -> Array[WaveletScaleSummary] {
  let result = Array::new(capacity=scalogram.scales.length())
  for scale_index = 0
      scale_index < scalogram.scales.length()
      scale_index = scale_index + 1 {
    let mut energy = 0.0
    let mut maximum = 0.0
    let mut absolute = 0.0
    let mut active = 0
    for sample = 0; sample < scalogram.sample_count; sample = sample + 1 {
      let value = scalogram_at(scalogram, scale_index, sample)
      let magnitude = value.abs()
      energy += value * value
      absolute += magnitude
      maximum = maximum.max(magnitude)
      if magnitude >= active_threshold {
        active += 1
      }
    }
    let count = scalogram.sample_count.to_double().max(1.0)
    result.push({
      scale: scalogram.scales[scale_index],
      pseudo_frequency: wavelet_pseudo_frequency(
        scalogram.scales[scale_index],
        sample_rate,
        scalogram.family,
      ),
      energy,
      maximum_magnitude: maximum,
      mean_absolute: absolute / count,
      active_fraction: active.to_double() / count,
    })
  }
  result
}

///|
/// Normalize wavelet scale energies to a probability distribution.
pub fn wavelet_energy_probabilities(scalogram : Scalogram) -> Array[Double] {
  let energy = wavelet_scale_energy(scalogram)
  let total = sum_samples(energy)
  if total <= 0.0 {
    Array::make(energy.length(), 0.0)
  } else {
    energy.map(value => value / total)
  }
}

///|
/// Compute Shannon entropy over wavelet scale energy.
pub fn wavelet_energy_entropy(scalogram : Scalogram) -> Double {
  let probabilities = wavelet_energy_probabilities(scalogram)
  let mut entropy = 0.0
  for probability in probabilities {
    if probability > 0.0 {
      entropy -= probability * @math.ln(probability)
    }
  }
  entropy
}

///|
/// Track the strongest scale at every sample.
pub fn wavelet_ridge(scalogram : Scalogram) -> Array[WaveletCoefficient] {
  let result = Array::new(capacity=scalogram.sample_count)
  for sample = 0; sample < scalogram.sample_count; sample = sample + 1 {
    let mut best_scale = 0
    let mut best_value = 0.0
    for scale_index = 0
        scale_index < scalogram.scales.length()
        scale_index = scale_index + 1 {
      let value = scalogram_at(scalogram, scale_index, sample)
      if value.abs() > best_value.abs() {
        best_scale = scale_index
        best_value = value
      }
    }
    result.push({
      sample_index: sample,
      scale: if scalogram.scales.is_empty() {
        0.0
      } else {
        scalogram.scales[best_scale]
      },
      value: best_value,
      magnitude: best_value.abs(),
    })
  }
  result
}

///|
/// Smooth a wavelet ridge scale with a median filter.
pub fn smooth_wavelet_ridge(
  ridge : Array[WaveletCoefficient],
  radius : Int,
) -> Array[WaveletCoefficient] {
  let result = Array::new(capacity=ridge.length())
  let window = radius.max(0)
  for index = 0; index < ridge.length(); index = index + 1 {
    let first = (index - window).max(0)
    let last = (index + window + 1).min(ridge.length())
    let scales = slice_samples(ridge.map(item => item.scale), first, last)
    let current = ridge[index]
    result.push({
      sample_index: current.sample_index,
      scale: median_samples(scales),
      value: current.value,
      magnitude: current.magnitude,
    })
  }
  result
}

///|
/// Estimate cross-scale coherence at each time sample.
pub fn wavelet_cross_scale_coherence(
  scalogram : Scalogram,
  low_scale_index : Int,
  high_scale_index : Int,
) -> Array[Double] {
  let result = Array::make(scalogram.sample_count, 0.0)
  for sample = 0; sample < scalogram.sample_count; sample = sample + 1 {
    let low = scalogram_at(scalogram, low_scale_index, sample)
    let high = scalogram_at(scalogram, high_scale_index, sample)
    result[sample] = if low.abs() + high.abs() <= 1.0e-20 {
      0.0
    } else {
      2.0 * low.abs().min(high.abs()) / (low.abs() + high.abs())
    }
  }
  result
}

///|
/// Find the first sustained ridge activation above a magnitude threshold.
pub fn wavelet_onset_index(
  ridge : Array[WaveletCoefficient],
  threshold : Double,
  required_samples : Int,
) -> Int {
  let required = required_samples.max(1)
  if ridge.length() < required {
    return -1
  }
  for index = 0; index <= ridge.length() - required; index = index + 1 {
    let mut active = true
    for offset = 0; offset < required; offset = offset + 1 {
      if ridge[index + offset].magnitude < threshold {
        active = false
        break
      }
    }
    if active {
      return index
    }
  }
  -1
}

///|
/// Build event features from a scalogram and its ridge.
pub fn wavelet_event_features(
  scalogram : Scalogram,
  sample_rate : Double,
  magnitude_threshold : Double,
) -> WaveletEventFeatures {
  let ridge = wavelet_ridge(scalogram)
  let strongest = strongest_wavelet_coefficient(scalogram)
  let onset = wavelet_onset_index(ridge, magnitude_threshold, 3)
  let mut last = onset
  if onset >= 0 {
    while last < ridge.length() && ridge[last].magnitude >= magnitude_threshold {
      last += 1
    }
  }
  let coherence = if scalogram.scales.length() < 2 {
    0.0
  } else {
    mean_samples(
      wavelet_cross_scale_coherence(scalogram, 0, scalogram.scales.length() - 1),
    )
  }
  {
    dominant_scale: dominant_wavelet_scale(scalogram),
    dominant_frequency: wavelet_pseudo_frequency(
      dominant_wavelet_scale(scalogram),
      sample_rate,
      scalogram.family,
    ),
    peak_time_index: strongest.sample_index,
    peak_magnitude: strongest.magnitude,
    entropy: wavelet_energy_entropy(scalogram),
    onset_index: onset,
    duration_samples: if onset < 0 {
      0
    } else {
      last - onset
    },
    coherence,
  }
}

///|
/// Return a soft threshold based on median absolute coefficients.
pub fn wavelet_universal_threshold(scalogram : Scalogram) -> Double {
  let values = scalogram.coefficients.map(value => value.abs())
  let noise = median_samples(values) / 0.6745
  noise * (2.0 * @math.ln(scalogram.sample_count.max(2).to_double())).sqrt()
}

///|
/// Export scale summaries as CSV.
pub fn wavelet_scale_summary_to_csv(
  summaries : Array[WaveletScaleSummary],
) -> String {
  let mut output = "scale,pseudo_frequency,energy,maximum_magnitude,mean_absolute,active_fraction\n"
  for summary in summaries {
    output = output +
      summary.scale.to_string() +
      "," +
      summary.pseudo_frequency.to_string() +
      "," +
      summary.energy.to_string() +
      "," +
      summary.maximum_magnitude.to_string() +
      "," +
      summary.mean_absolute.to_string() +
      "," +
      summary.active_fraction.to_string() +
      "\n"
  }
  output
}

///|
/// Format multi-resolution event features.
pub fn wavelet_event_features_text(features : WaveletEventFeatures) -> String {
  "scale=" +
  features.dominant_scale.to_string() +
  " frequency=" +
  features.dominant_frequency.to_string() +
  " peak_index=" +
  features.peak_time_index.to_string() +
  " peak=" +
  features.peak_magnitude.to_string() +
  " entropy=" +
  features.entropy.to_string() +
  " onset=" +
  features.onset_index.to_string() +
  " duration=" +
  features.duration_samples.to_string() +
  " coherence=" +
  features.coherence.to_string()
}