///|
/// Statistical summary of an event catalog.
pub struct CatalogStatistics {
  count : Int
  minimum_magnitude : Double
  maximum_magnitude : Double
  mean_magnitude : Double
  median_magnitude : Double
  standard_deviation : Double
  completeness_magnitude : Double
  b_value : Double
  a_value : Double
  mean_depth_km : Double
  mean_station_count : Double
  total_duration : Double
}

///|
/// A magnitude-frequency bin.
pub struct MagnitudeBin {
  lower : Double
  upper : Double
  count : Int
  cumulative_count : Int
  log10_cumulative : Double
}

///|
/// A temporal recurrence interval.
pub struct RecurrenceInterval {
  first_event_time : Double
  second_event_time : Double
  seconds : Double
  first_magnitude : Double
  second_magnitude : Double
}

///|
/// Return catalog magnitudes in insertion order.
pub fn catalog_magnitudes(catalog : EventCatalog) -> Array[Double] {
  catalog.entries.map(entry => entry.magnitude)
}

///|
/// Return catalog origin times in insertion order.
pub fn catalog_origin_times(catalog : EventCatalog) -> Array[Double] {
  catalog.entries.map(entry => entry.origin_time)
}

///|
/// Estimate a completeness magnitude from a histogram mode.
pub fn estimate_completeness_magnitude(
  catalog : EventCatalog,
  bin_width : Double,
) -> Double {
  let width = bin_width.abs().max(0.01)
  if catalog.entries.is_empty() {
    return 0.0
  }
  let minimum = minimum_sample(catalog_magnitudes(catalog))
  let bins = Array::make(catalog.entries.length() + 1, 0)
  for entry in catalog.entries {
    let index = ((entry.magnitude - minimum) / width).to_int().max(0)
    if index < bins.length() {
      bins[index] += 1
    }
  }
  let mut best = 0
  for index = 1; index < bins.length(); index = index + 1 {
    if bins[index] > bins[best] {
      best = index
    }
  }
  minimum + best.to_double() * width
}

///|
/// Estimate Gutenberg-Richter b-value above a completeness magnitude.
pub fn catalog_b_value(catalog : EventCatalog, completeness : Double) -> Double {
  let values = catalog.entries
    .filter(entry => entry.magnitude >= completeness)
    .map(entry => entry.magnitude)
  if values.is_empty() {
    return 0.0
  }
  let mean = mean_samples(values)
  let correction = 0.05
  let denominator = mean - completeness + correction
  if denominator <= 0.0 {
    0.0
  } else {
    @math.log10(@math.exp(1.0)) / denominator
  }
}

///|
/// Estimate the Gutenberg-Richter intercept from count and b-value.
pub fn catalog_a_value(
  catalog : EventCatalog,
  completeness : Double,
  b_value : Double,
) -> Double {
  let count = catalog.entries
    .filter(entry => entry.magnitude >= completeness)
    .length()
  if count == 0 {
    0.0
  } else {
    @math.log10(count.to_double()) + b_value * completeness
  }
}

///|
/// Compute magnitude-frequency bins and cumulative counts.
pub fn magnitude_frequency_bins(
  catalog : EventCatalog,
  bin_width : Double,
) -> Array[MagnitudeBin] {
  let result = Array::new()
  let width = bin_width.abs().max(0.01)
  if catalog.entries.is_empty() {
    return result
  }
  let magnitudes = catalog_magnitudes(catalog)
  let minimum = (minimum_sample(magnitudes) / width).floor() * width
  let maximum = maximum_sample(magnitudes)
  let count = ((maximum - minimum) / width).to_int().max(0) + 1
  let counts = Array::make(count, 0)
  for magnitude in magnitudes {
    let index = ((magnitude - minimum) / width).to_int().max(0).min(count - 1)
    counts[index] += 1
  }
  let mut cumulative = 0
  for index = count - 1; index >= 0; index = index - 1 {
    cumulative += counts[index]
    result.push({
      lower: minimum + index.to_double() * width,
      upper: minimum + (index + 1).to_double() * width,
      count: counts[index],
      cumulative_count: cumulative,
      log10_cumulative: if cumulative == 0 {
        0.0
      } else {
        @math.log10(cumulative.to_double())
      },
    })
    if index == 0 {
      break
    }
  }
  result.rev_in_place()
  result
}

///|
/// Compute time differences between consecutive catalog entries.
pub fn recurrence_intervals(
  catalog : EventCatalog,
) -> Array[RecurrenceInterval] {
  let result = Array::new()
  if catalog.entries.length() < 2 {
    return result
  }
  let entries = catalog.entries.copy()
  entries.sort_by((left, right) => {
    if left.origin_time < right.origin_time {
      -1
    } else if left.origin_time > right.origin_time {
      1
    } else {
      0
    }
  })
  for index = 1; index < entries.length(); index = index + 1 {
    let left = entries[index - 1]
    let right = entries[index]
    result.push({
      first_event_time: left.origin_time,
      second_event_time: right.origin_time,
      seconds: right.origin_time - left.origin_time,
      first_magnitude: left.magnitude,
      second_magnitude: right.magnitude,
    })
  }
  result
}

///|
/// Return the mean recurrence rate in events per day.
pub fn recurrence_rate_per_day(catalog : EventCatalog) -> Double {
  if catalog.entries.length() < 2 {
    return 0.0
  }
  let intervals = recurrence_intervals(catalog)
  let total = sum_samples(intervals.map(item => item.seconds))
  if total <= 0.0 {
    0.0
  } else {
    (catalog.entries.length() - 1).to_double() * 86400.0 / total
  }
}

///|
/// Calculate the catalog bounding box as lat/lon/depth extrema.
pub fn catalog_bounds(catalog : EventCatalog) -> Array[Double] {
  if catalog.entries.is_empty() {
    return [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
  }
  let first = catalog.entries[0]
  let mut min_lat = first.latitude
  let mut max_lat = first.latitude
  let mut min_lon = first.longitude
  let mut max_lon = first.longitude
  let mut min_depth = first.depth_km
  let mut max_depth = first.depth_km
  for entry in catalog.entries {
    min_lat = min_lat.min(entry.latitude)
    max_lat = max_lat.max(entry.latitude)
    min_lon = min_lon.min(entry.longitude)
    max_lon = max_lon.max(entry.longitude)
    min_depth = min_depth.min(entry.depth_km)
    max_depth = max_depth.max(entry.depth_km)
  }
  [min_lat, max_lat, min_lon, max_lon, min_depth, max_depth]
}

///|
/// Compute a full catalog statistics record.
pub fn describe_catalog(catalog : EventCatalog) -> CatalogStatistics {
  let magnitudes = catalog_magnitudes(catalog)
  let depths = catalog.entries.map(entry => entry.depth_km)
  let stations = catalog.entries.map(entry => entry.station_count.to_double())
  let completeness = estimate_completeness_magnitude(catalog, 0.1)
  let b_value = catalog_b_value(catalog, completeness)
  let times = catalog_origin_times(catalog)
  {
    count: catalog.entries.length(),
    minimum_magnitude: minimum_sample(magnitudes),
    maximum_magnitude: maximum_sample(magnitudes),
    mean_magnitude: mean_samples(magnitudes),
    median_magnitude: median_samples(magnitudes),
    standard_deviation: standard_deviation(magnitudes),
    completeness_magnitude: completeness,
    b_value,
    a_value: catalog_a_value(catalog, completeness, b_value),
    mean_depth_km: mean_samples(depths),
    mean_station_count: mean_samples(stations),
    total_duration: if times.length() < 2 {
      0.0
    } else {
      maximum_sample(times) - minimum_sample(times)
    },
  }
}

///|
/// Filter a catalog by a rectangular spatial region.
pub fn catalog_in_region(
  catalog : EventCatalog,
  latitude_min : Double,
  latitude_max : Double,
  longitude_min : Double,
  longitude_max : Double,
) -> Array[CatalogEntry] {
  catalog.entries.filter(entry => {
    entry.latitude >= latitude_min.min(latitude_max) &&
    entry.latitude <= latitude_min.max(latitude_max) &&
    entry.longitude >= longitude_min.min(longitude_max) &&
    entry.longitude <= longitude_min.max(longitude_max)
  })
}

///|
/// Return events within a radius of a geographic center.
pub fn catalog_near(
  catalog : EventCatalog,
  center : StationPosition,
  radius_km : Double,
) -> Array[CatalogEntry] {
  catalog.entries.filter(entry => {
    let position = StationPosition::new(
      "event",
      entry.latitude,
      entry.longitude,
      -entry.depth_km,
    )
    surface_distance_km(center, position) <= radius_km.abs()
  })
}

///|
/// Compute the fraction of events meeting a magnitude threshold.
pub fn magnitude_exceedance_fraction(
  catalog : EventCatalog,
  threshold : Double,
) -> Double {
  if catalog.entries.is_empty() {
    0.0
  } else {
    catalog.entries
    .filter(entry => entry.magnitude >= threshold)
    .length()
    .to_double() /
    catalog.entries.length().to_double()
  }
}

///|
/// Format catalog statistics as a single stable report line.
pub fn catalog_statistics_text(statistics : CatalogStatistics) -> String {
  "count=" +
  statistics.count.to_string() +
  " min_mag=" +
  statistics.minimum_magnitude.to_string() +
  " max_mag=" +
  statistics.maximum_magnitude.to_string() +
  " mean_mag=" +
  statistics.mean_magnitude.to_string() +
  " median_mag=" +
  statistics.median_magnitude.to_string() +
  " completeness=" +
  statistics.completeness_magnitude.to_string() +
  " b_value=" +
  statistics.b_value.to_string() +
  " a_value=" +
  statistics.a_value.to_string() +
  " mean_depth=" +
  statistics.mean_depth_km.to_string() +
  " catalog_duration=" +
  statistics.total_duration.to_string()
}