///|
/// Interpolation strategy used by waveform resampling and gap repair.
pub enum InterpolationMethod {
  LinearInterpolation
  NearestInterpolation
  HoldInterpolation
}

///|
/// Policy for replacing missing or invalid samples.
pub enum GapFillPolicy {
  FillWithZero
  FillWithMean
  FillWithLinear
  FillWithPrevious
}

///|
/// An inclusive run of samples identified as a gap or invalid segment.
pub struct GapInterval {
  start_index : Int
  end_index : Int
  length : Int
  start_time : Double
  end_time : Double
}

///|
/// Diagnostics returned by a resampling operation.
pub struct ResampleReport {
  old_rate : Double
  new_rate : Double
  input_samples : Int
  output_samples : Int
  duration : Double
  interpolation : InterpolationMethod
  endpoint_error : Double
}

///|
/// Three-component sample at one time index.
pub struct ComponentSample {
  north : Double
  east : Double
  vertical : Double
}

///|
/// Three-component waveform bundle with aligned metadata.
pub struct ComponentTriplet {
  station : String
  sampling_rate : Double
  start_time : Double
  north : Array[Double]
  east : Array[Double]
  vertical : Array[Double]
}

///|
/// Construct an aligned three-component waveform bundle.
pub fn ComponentTriplet::new(
  station : String,
  sampling_rate : Double,
  start_time : Double,
  north : Array[Double],
  east : Array[Double],
  vertical : Array[Double],
) -> ComponentTriplet {
  { station, sampling_rate, start_time, north, east, vertical }
}

///|
/// Horizontal components after rotation into radial and transverse axes.
pub struct RotatedComponents {
  radial : Array[Double]
  transverse : Array[Double]
  back_azimuth : Double
}

///|
/// Result of a deterministic trace transformation.
pub struct TransformReport {
  input_samples : Int
  output_samples : Int
  input_mean : Double
  output_mean : Double
  input_rms : Double
  output_rms : Double
  peak_before : Double
  peak_after : Double
}

///|
/// Return all interpolation strategies supported by the portable backend.
pub fn supported_interpolation_methods() -> Array[InterpolationMethod] {
  [LinearInterpolation, NearestInterpolation, HoldInterpolation]
}

///|
/// Return the linear interpolation strategy.
pub fn interpolation_linear() -> InterpolationMethod {
  LinearInterpolation
}

///|
/// Return the nearest-neighbor interpolation strategy.
pub fn interpolation_nearest() -> InterpolationMethod {
  NearestInterpolation
}

///|
/// Return all gap policies supported by the QC pipeline.
pub fn supported_gap_policies() -> Array[GapFillPolicy] {
  [FillWithZero, FillWithMean, FillWithLinear, FillWithPrevious]
}

///|
/// Return the linear gap-fill policy.
pub fn gap_fill_linear_policy() -> GapFillPolicy {
  FillWithLinear
}

///|
/// Return a reversed copy without mutating the input array.
pub fn reverse_samples(data : Array[Double]) -> Array[Double] {
  let result = data.copy()
  result.rev_in_place()
  result
}

///|
/// Return a copy with a constant offset removed.
pub fn remove_offset(data : Array[Double], offset : Double) -> Array[Double] {
  let result = data.copy()
  for index = 0; index < result.length(); index = index + 1 {
    result[index] -= offset
  }
  result
}

///|
/// Remove a least-squares line from a copied signal.
pub fn detrend_copy(data : Array[Double]) -> Array[Double] {
  let result = data.copy()
  if result.length() <= 1 {
    return result
  }
  let n = result.length()
  let mean_y = mean_samples(result)
  let mean_x = (n - 1).to_double() / 2.0
  let mut covariance = 0.0
  let mut variance = 0.0
  for index = 0; index < n; index = index + 1 {
    let x = index.to_double() - mean_x
    covariance += x * (result[index] - mean_y)
    variance += x * x
  }
  if variance <= 0.0 {
    return result
  }
  let slope = covariance / variance
  let intercept = mean_y - slope * mean_x
  for index = 0; index < n; index = index + 1 {
    result[index] -= slope * index.to_double() + intercept
  }
  result
}

///|
/// Return a moving arithmetic mean with edge windows clipped to the signal.
pub fn moving_average_samples(
  data : Array[Double],
  radius : Int,
) -> Array[Double] {
  let result = Array::make(data.length(), 0.0)
  let safe_radius = if radius < 0 { 0 } else { radius }
  for index = 0; index < data.length(); index = index + 1 {
    let first = if index < safe_radius { 0 } else { index - safe_radius }
    let last = if index + safe_radius + 1 > data.length() {
      data.length()
    } else {
      index + safe_radius + 1
    }
    let mut total = 0.0
    for cursor = first; cursor < last; cursor = cursor + 1 {
      total += data[cursor]
    }
    result[index] = if last == first {
      0.0
    } else {
      total / (last - first).to_double()
    }
  }
  result
}

///|
/// Return a moving root-mean-square amplitude.
pub fn moving_rms(data : Array[Double], radius : Int) -> Array[Double] {
  let squares = data.map(value => value * value)
  moving_average_samples(squares, radius).map(value => value.sqrt())
}

///|
/// Return a moving median, useful for impulsive noise suppression.
pub fn moving_median(data : Array[Double], radius : Int) -> Array[Double] {
  let result = Array::make(data.length(), 0.0)
  let safe_radius = if radius < 0 { 0 } else { radius }
  for index = 0; index < data.length(); index = index + 1 {
    let first = if index < safe_radius { 0 } else { index - safe_radius }
    let last = if index + safe_radius + 1 > data.length() {
      data.length()
    } else {
      index + safe_radius + 1
    }
    result[index] = median_samples(slice_samples(data, first, last))
  }
  result
}

///|
/// Return the median absolute deviation in each moving window.
pub fn moving_mad(data : Array[Double], radius : Int) -> Array[Double] {
  let result = Array::make(data.length(), 0.0)
  let safe_radius = if radius < 0 { 0 } else { radius }
  for index = 0; index < data.length(); index = index + 1 {
    let first = if index < safe_radius { 0 } else { index - safe_radius }
    let last = if index + safe_radius + 1 > data.length() {
      data.length()
    } else {
      index + safe_radius + 1
    }
    let values = slice_samples(data, first, last)
    let center = median_samples(values)
    let distances = values.map(value => (value - center).abs())
    result[index] = median_samples(distances)
  }
  result
}

///|
/// Apply a deadband while preserving signal polarity.
pub fn apply_deadband(
  data : Array[Double],
  threshold : Double,
) -> Array[Double] {
  let result = data.copy()
  let limit = threshold.abs()
  for index = 0; index < result.length(); index = index + 1 {
    if result[index].abs() < limit {
      result[index] = 0.0
    }
  }
  result
}

///|
/// Winsorize extreme values to quantile limits.
pub fn winsorize_samples(
  data : Array[Double],
  lower_probability : Double,
  upper_probability : Double,
) -> Array[Double] {
  let result = data.copy()
  let low = quantile_samples(data, lower_probability)
  let high = quantile_samples(data, upper_probability.max(lower_probability))
  for index = 0; index < result.length(); index = index + 1 {
    result[index] = result[index].max(low).min(high)
  }
  result
}

///|
/// Return indexes whose absolute amplitude exceeds a threshold.
pub fn threshold_indexes(
  data : Array[Double],
  threshold : Double,
) -> Array[Int] {
  let result = Array::new()
  for index = 0; index < data.length(); index = index + 1 {
    if data[index].abs() >= threshold.abs() {
      result.push(index)
    }
  }
  result
}

///|
/// Find contiguous runs of exact zero values.
pub fn zero_runs(data : Array[Double]) -> Array[(Int, Int)] {
  let result = Array::new()
  let mut start = -1
  for index = 0; index < data.length(); index = index + 1 {
    if data[index] == 0.0 && start < 0 {
      start = index
    }
    if (data[index] != 0.0 || index == data.length() - 1) && start >= 0 {
      let end = if data[index] != 0.0 { index } else { index + 1 }
      result.push((start, end))
      start = -1
    }
  }
  result
}

///|
/// Locate non-finite placeholders represented by a configured sentinel.
pub fn find_sentinel_gaps(
  data : Array[Double],
  sentinel : Double,
) -> Array[(Int, Int)] {
  let result = Array::new()
  let mut start = -1
  for index = 0; index < data.length(); index = index + 1 {
    if data[index] == sentinel && start < 0 {
      start = index
    }
    if (data[index] != sentinel || index == data.length() - 1) && start >= 0 {
      let end = if data[index] != sentinel { index } else { index + 1 }
      result.push((start, end))
      start = -1
    }
  }
  result
}

///|
/// Convert half-open gap indexes into physical metadata.
pub fn describe_gap(
  start : Int,
  end : Int,
  sample_rate : Double,
  start_time : Double,
) -> GapInterval {
  let safe_start = start.max(0)
  let safe_end = end.max(safe_start)
  let rate = sample_rate.max(1.0e-12)
  {
    start_index: safe_start,
    end_index: safe_end,
    length: safe_end - safe_start,
    start_time: start_time + safe_start.to_double() / rate,
    end_time: start_time + safe_end.to_double() / rate,
  }
}

///|
/// Find runs of a sentinel and attach sample timing.
pub fn find_gaps(
  data : Array[Double],
  sentinel : Double,
  sample_rate : Double,
  start_time : Double,
) -> Array[GapInterval] {
  find_sentinel_gaps(data, sentinel).map(item => {
    describe_gap(item.0, item.1, sample_rate, start_time)
  })
}

///|
/// Interpolate one sample position from a source array.
pub fn interpolate_sample(
  data : Array[Double],
  position : Double,
  strategy : InterpolationMethod,
) -> Double {
  if data.is_empty() {
    return 0.0
  }
  if position <= 0.0 {
    return data[0]
  }
  if position >= (data.length() - 1).to_double() {
    return data[data.length() - 1]
  }
  let left = position.to_int()
  let fraction = position - left.to_double()
  match strategy {
    LinearInterpolation => data[left] + fraction * (data[left + 1] - data[left])
    NearestInterpolation =>
      if fraction < 0.5 {
        data[left]
      } else {
        data[left + 1]
      }
    HoldInterpolation => data[left]
  }
}

///|
/// Resample using a selected interpolation strategy.
pub fn resample_with_strategy(
  data : Array[Double],
  old_rate : Double,
  new_rate : Double,
  strategy : InterpolationMethod,
) -> Array[Double] {
  let result = Array::new()
  if data.is_empty() || old_rate <= 0.0 || new_rate <= 0.0 {
    return result
  }
  if data.length() == 1 {
    result.push(data[0])
    return result
  }
  let duration = (data.length() - 1).to_double() / old_rate
  let count = (duration * new_rate).to_int() + 1
  for index = 0; index < count; index = index + 1 {
    let position = index.to_double() * old_rate / new_rate
    result.push(interpolate_sample(data, position, strategy))
  }
  result
}

///|
/// Resample a trace and report endpoint preservation error.
pub fn resample_trace_report(
  trace : Trace,
  new_rate : Double,
  strategy : InterpolationMethod,
) -> (Trace, ResampleReport) {
  let values = resample_with_strategy(
    trace.data,
    trace.sampling_rate,
    new_rate,
    strategy,
  )
  let safe_rate = if new_rate > 0.0 { new_rate } else { trace.sampling_rate }
  let output = Trace::new(
    trace.station,
    trace.channel,
    safe_rate,
    trace.start_time,
    values,
  )
  let endpoint_error = if trace.data.is_empty() || values.is_empty() {
    0.0
  } else {
    (values[values.length() - 1] - trace.data[trace.data.length() - 1]).abs()
  }
  let report = {
    old_rate: trace.sampling_rate,
    new_rate: safe_rate,
    input_samples: trace.data.length(),
    output_samples: values.length(),
    duration: trace.duration(),
    interpolation: strategy,
    endpoint_error,
  }
  (output, report)
}

///|
/// Fill a half-open interval in place with linear endpoint interpolation.
pub fn fill_gap_linear(data : Array[Double], start : Int, end : Int) -> Unit {
  if data.is_empty() {
    return
  }
  let first = start.max(0).min(data.length())
  let last = end.max(first).min(data.length())
  if first >= last {
    return
  }
  let left_index = first - 1
  let right_index = last
  let left = if left_index >= 0 { data[left_index] } else { 0.0 }
  let right = if right_index < data.length() { data[right_index] } else { left }
  let span = (last - first + 1).to_double()
  for index = first; index < last; index = index + 1 {
    data[index] = left + (right - left) * (index - first + 1).to_double() / span
  }
}

///|
/// Fill a half-open interval in place according to a production QC policy.
pub fn fill_gap(
  data : Array[Double],
  start : Int,
  end : Int,
  policy : GapFillPolicy,
) -> Unit {
  let first = start.max(0).min(data.length())
  let last = end.max(first).min(data.length())
  match policy {
    FillWithZero =>
      for index = first; index < last; index = index + 1 {
        data[index] = 0.0
      }
    FillWithMean => {
      let value = mean_samples(data)
      for index = first; index < last; index = index + 1 {
        data[index] = value
      }
    }
    FillWithPrevious => {
      let value = if first > 0 { data[first - 1] } else { 0.0 }
      for index = first; index < last; index = index + 1 {
        data[index] = value
      }
    }
    FillWithLinear => fill_gap_linear(data, first, last)
  }
}

///|
/// Fill all sentinel gaps and return a repaired copy.
pub fn repair_sentinel_gaps(
  data : Array[Double],
  sentinel : Double,
  policy : GapFillPolicy,
) -> Array[Double] {
  let result = data.copy()
  for gap in find_sentinel_gaps(result, sentinel) {
    fill_gap(result, gap.0, gap.1, policy)
  }
  result
}

///|
/// Apply a symmetric integer delay to a component with zero padding.
pub fn delay_component(data : Array[Double], delay : Int) -> Array[Double] {
  let result = Array::make(data.length(), 0.0)
  for index = 0; index < data.length(); index = index + 1 {
    let source = index - delay
    if source >= 0 && source < data.length() {
      result[index] = data[source]
    }
  }
  result
}

///|
/// Rotate north/east components into radial/transverse components.
pub fn rotate_horizontal(
  north : Array[Double],
  east : Array[Double],
  back_azimuth_degrees : Double,
) -> RotatedComponents {
  let count = north.length().min(east.length())
  let radial = Array::make(count, 0.0)
  let transverse = Array::make(count, 0.0)
  let angle = back_azimuth_degrees * @math.PI / 180.0
  let cosine = @math.cos(angle)
  let sine = @math.sin(angle)
  for index = 0; index < count; index = index + 1 {
    radial[index] = north[index] * cosine + east[index] * sine
    transverse[index] = -north[index] * sine + east[index] * cosine
  }
  { radial, transverse, back_azimuth: back_azimuth_degrees }
}

///|
/// Return component-wise vector amplitude for a three-component record.
pub fn component_vector_amplitude(sample : ComponentSample) -> Double {
  (sample.north * sample.north +
  sample.east * sample.east +
  sample.vertical * sample.vertical).sqrt()
}

///|
/// Convert aligned component arrays into vector amplitudes.
pub fn triplet_vector_amplitude(triplet : ComponentTriplet) -> Array[Double] {
  let count = triplet.north
    .length()
    .min(triplet.east.length())
    .min(triplet.vertical.length())
  let result = Array::make(count, 0.0)
  for index = 0; index < count; index = index + 1 {
    result[index] = component_vector_amplitude({
      north: triplet.north[index],
      east: triplet.east[index],
      vertical: triplet.vertical[index],
    })
  }
  result
}

///|
/// Normalize each component by one shared vector peak.
pub fn normalize_triplet(triplet : ComponentTriplet) -> ComponentTriplet {
  let amplitudes = triplet_vector_amplitude(triplet)
  let scale = maximum_sample(amplitudes)
  let safe_scale = if scale <= 0.0 { 1.0 } else { scale }
  {
    station: triplet.station,
    sampling_rate: triplet.sampling_rate,
    start_time: triplet.start_time,
    north: triplet.north.map(value => value / safe_scale),
    east: triplet.east.map(value => value / safe_scale),
    vertical: triplet.vertical.map(value => value / safe_scale),
  }
}

///|
/// Produce a report comparing two signal arrays.
pub fn compare_transform(
  input : Array[Double],
  output : Array[Double],
) -> TransformReport {
  {
    input_samples: input.length(),
    output_samples: output.length(),
    input_mean: mean_samples(input),
    output_mean: mean_samples(output),
    input_rms: rms_samples(input),
    output_rms: rms_samples(output),
    peak_before: maximum_sample(input).abs().max(minimum_sample(input).abs()),
    peak_after: maximum_sample(output).abs().max(minimum_sample(output).abs()),
  }
}

///|
/// Clip a trace to a physical time interval using sample-grid rounding.
pub fn clip_trace_time(
  trace : Trace,
  start_time : Double,
  end_time : Double,
) -> Trace {
  if trace.sampling_rate <= 0.0 {
    return trace.clone()
  }
  let start = ((start_time - trace.start_time) * trace.sampling_rate).to_int()
  let end = ((end_time - trace.start_time) * trace.sampling_rate).to_int() + 1
  trace_slice(trace, start, end)
}

///|
/// Pad a trace to a requested duration, preserving the start time.
pub fn pad_trace_duration(
  trace : Trace,
  duration : Double,
  fill : Double,
) -> Trace {
  if trace.sampling_rate <= 0.0 || duration <= 0.0 {
    return trace.clone()
  }
  let required = (duration * trace.sampling_rate).to_int() + 1
  if required <= trace.npts() {
    return trace.clone()
  }
  let data = trace.data.copy()
  for index = trace.npts(); index < required; index = index + 1 {
    data.push(fill)
  }
  trace_with_data(trace, data)
}

///|
/// Normalize a trace with robust median/MAD statistics.
pub fn robust_normalize_trace(trace : Trace) -> Trace {
  let center = median_samples(trace.data)
  let deviations = trace.data.map(value => (value - center).abs())
  let scale = median_samples(deviations) * 1.4826
  let safe_scale = if scale <= 0.0 { 1.0 } else { scale }
  trace_with_data(trace, trace.data.map(value => (value - center) / safe_scale))
}

///|
/// Apply a Tukey-like soft limiter to reduce single-sample spikes.
pub fn soft_limit(data : Array[Double], limit : Double) -> Array[Double] {
  let result = data.copy()
  let safe_limit = limit.abs()
  if safe_limit <= 0.0 {
    result.fill(0.0)
    return result
  }
  for index = 0; index < result.length(); index = index + 1 {
    let value = result[index]
    let magnitude = value.abs()
    if magnitude > safe_limit {
      result[index] = value.signum() *
        safe_limit *
        (1.0 + @math.ln(1.0 + (magnitude - safe_limit) / safe_limit))
    }
  }
  result
}

///|
/// Return the component covariance terms for an aligned triplet window.
pub fn component_covariance(triplet : ComponentTriplet) -> Array[Double] {
  let count = triplet.north
    .length()
    .min(triplet.east.length())
    .min(triplet.vertical.length())
  if count == 0 {
    return [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
  }
  let north = mean_samples(triplet.north)
  let east = mean_samples(triplet.east)
  let vertical = mean_samples(triplet.vertical)
  let mut nn = 0.0
  let mut ee = 0.0
  let mut vv = 0.0
  let mut ne = 0.0
  let mut nv = 0.0
  let mut ev = 0.0
  for index = 0; index < count; index = index + 1 {
    let n = triplet.north[index] - north
    let e = triplet.east[index] - east
    let v = triplet.vertical[index] - vertical
    nn += n * n
    ee += e * e
    vv += v * v
    ne += n * e
    nv += n * v
    ev += e * v
  }
  let scale = 1.0 / count.to_double()
  [nn * scale, ee * scale, vv * scale, ne * scale, nv * scale, ev * scale]
}

///|
/// Quantize a sample array to a fixed decimal resolution for stable exports.
pub fn quantize_samples(data : Array[Double], decimals : Int) -> Array[Double] {
  let digits = if decimals < 0 { 0 } else { decimals }
  let scale = @math.pow(10.0, digits.to_double())
  if scale <= 0.0 {
    return data.copy()
  }
  data.map(value => (value * scale).round() / scale)
}

///|
/// Create a stable transform report line for QC logs.
pub fn transform_report_text(report : TransformReport) -> String {
  "input=" +
  report.input_samples.to_string() +
  " output=" +
  report.output_samples.to_string() +
  " input_mean=" +
  report.input_mean.to_string() +
  " output_mean=" +
  report.output_mean.to_string() +
  " input_rms=" +
  report.input_rms.to_string() +
  " output_rms=" +
  report.output_rms.to_string() +
  " peak_before=" +
  report.peak_before.to_string() +
  " peak_after=" +
  report.peak_after.to_string()
}