///|
/// Diagnostics for a sampled calibration curve.
pub struct CalibrationCurveDiagnostics {
  valid : Bool
  point_count : Int
  minimum_frequency : Double
  maximum_frequency : Double
  minimum_sensitivity : Double
  maximum_sensitivity : Double
  dynamic_range_db : Double
  monotonic : Bool
}

///|
/// A calibration application result with trace statistics.
pub struct CalibrationApplicationReport {
  input_peak : Double
  output_peak : Double
  input_rms : Double
  output_rms : Double
  representative_frequency : Double
  sensitivity : Double
  scale_factor : Double
  stable : Bool
}

///|
/// Inspect a calibration curve for dynamic range and ordering.
pub fn calibration_curve_diagnostics(
  curve : CalibrationCurve,
) -> CalibrationCurveDiagnostics {
  let frequencies = curve.frequencies
  let sensitivities = curve.sensitivities
  let mut monotonic = true
  for index = 1; index < frequencies.length(); index = index + 1 {
    if frequencies[index] <= frequencies[index - 1] {
      monotonic = false
    }
  }
  let minimum = minimum_sample(sensitivities)
  let maximum = maximum_sample(sensitivities)
  {
    valid: frequencies.length() == sensitivities.length() &&
    !frequencies.is_empty() &&
    monotonic &&
    minimum > 0.0,
    point_count: frequencies.length(),
    minimum_frequency: minimum_sample(frequencies),
    maximum_frequency: maximum_sample(frequencies),
    minimum_sensitivity: minimum,
    maximum_sensitivity: maximum,
    dynamic_range_db: if minimum <= 0.0 {
      0.0
    } else {
      20.0 * @math.log10(maximum / minimum)
    },
    monotonic,
  }
}

///|
/// Compute log-frequency slopes between curve points.
pub fn calibration_log_slopes(curve : CalibrationCurve) -> Array[Double] {
  let result = Array::new()
  for index = 1
      index < curve.frequencies.length() && index < curve.sensitivities.length()
      index = index + 1 {
    let frequency_ratio = @math.ln(
      curve.frequencies[index] / curve.frequencies[index - 1],
    )
    let sensitivity_ratio = @math.ln(
      curve.sensitivities[index] / curve.sensitivities[index - 1],
    )
    result.push(
      if frequency_ratio.abs() <= 1.0e-12 {
        0.0
      } else {
        sensitivity_ratio / frequency_ratio
      },
    )
  }
  result
}

///|
/// Smooth a calibration curve with a centered moving average in log space.
pub fn smooth_calibration_curve(
  curve : CalibrationCurve,
  radius : Int,
) -> CalibrationCurve {
  let result = curve.sensitivities.copy()
  let window = radius.max(0)
  for index = 0; index < result.length(); index = index + 1 {
    let first = (index - window).max(0)
    let last = (index + window + 1).min(result.length())
    let values = slice_samples(result, first, last)
    result[index] = mean_samples(values).max(1.0e-20)
  }
  {
    frequencies: curve.frequencies.copy(),
    sensitivities: result,
    unit: curve.unit,
  }
}

///|
/// Apply a curve and calculate a calibration report.
pub fn apply_calibration_report(
  trace : Trace,
  curve : CalibrationCurve,
  representative_frequency : Double,
) -> (Trace, CalibrationApplicationReport) {
  let sensitivity = curve.sensitivity(representative_frequency)
  let output = calibrate_with_curve(trace, curve, representative_frequency)
  let report = {
    input_peak: absolute_peak(trace.data),
    output_peak: absolute_peak(output.data),
    input_rms: rms_samples(trace.data),
    output_rms: rms_samples(output.data),
    representative_frequency,
    sensitivity,
    scale_factor: if sensitivity <= 0.0 {
      0.0
    } else {
      1.0 / sensitivity
    },
    stable: sensitivity > 0.0 && trace_is_well_formed(output),
  }
  (output, report)
}

///|
/// Combine two calibration curves at shared frequencies.
pub fn combine_calibration_curves(
  left : CalibrationCurve,
  right : CalibrationCurve,
  right_weight : Double,
) -> CalibrationCurve {
  let count = left.frequencies.length().min(right.frequencies.length())
  let frequencies = Array::new(capacity=count)
  let sensitivities = Array::new(capacity=count)
  let weight = right_weight.max(0.0).min(1.0)
  for index = 0; index < count; index = index + 1 {
    frequencies.push((left.frequencies[index] + right.frequencies[index]) / 2.0)
    sensitivities.push(
      left.sensitivities[index] * (1.0 - weight) +
      right.sensitivities[index] * weight,
    )
  }
  match CalibrationCurve::new(frequencies, sensitivities, left.unit) {
    Ok(curve) => curve
    Err(_) => left
  }
}

///|
/// Export calibration points and derived slope diagnostics.
pub fn calibration_curve_to_csv(curve : CalibrationCurve) -> String {
  let slopes = calibration_log_slopes(curve)
  let mut output = "frequency,sensitivity,log_slope\n"
  for index = 0; index < curve.frequencies.length(); index = index + 1 {
    let slope = if index == 0 || index - 1 >= slopes.length() {
      0.0
    } else {
      slopes[index - 1]
    }
    output = output +
      curve.frequencies[index].to_string() +
      "," +
      curve.sensitivities[index].to_string() +
      "," +
      slope.to_string() +
      "\n"
  }
  output
}

///|
/// Format calibration diagnostics for the release report.
pub fn calibration_curve_diagnostics_text(
  diagnostics : CalibrationCurveDiagnostics,
) -> String {
  "valid=" +
  diagnostics.valid.to_string() +
  " points=" +
  diagnostics.point_count.to_string() +
  " frequency_min=" +
  diagnostics.minimum_frequency.to_string() +
  " frequency_max=" +
  diagnostics.maximum_frequency.to_string() +
  " sensitivity_min=" +
  diagnostics.minimum_sensitivity.to_string() +
  " sensitivity_max=" +
  diagnostics.maximum_sensitivity.to_string() +
  " dynamic_range_db=" +
  diagnostics.dynamic_range_db.to_string() +
  " monotonic=" +
  diagnostics.monotonic.to_string()
}