///|
fn finite_sample_pair(a : FieldSample, b : FieldSample) -> Bool {
  is_finite(a.value) && is_finite(b.value)
}

///|
pub fn sample_pair_count(samples : ArrayView[FieldSample]) -> Int {
  let mut count = 0
  for i in 0.. Double {
  if index < 0 || index >= samples.length() {
    0.0
  } else {
    let mut best = 1.0e308
    for i in 0.. Double {
  if samples.length() < 2 {
    0.0
  } else {
    let mut total = 0.0
    for i in 0.. Array[VariogramBin] {
  let safe_bins = Int::max(1, bins)
  let safe_distance = Double::max(0.0, maximum_distance_m)
  let width = if safe_distance == 0.0 {
    1.0
  } else {
    safe_distance / Double::from_int(safe_bins)
  }
  let pair_counts : Array[Int] = []
  let distance_totals : Array[Double] = []
  let semivariance_totals : Array[Double] = []
  for _ in 0.. Double {
  let distance = Double::max(0.0, distance_m)
  if model.range_m <= 0.0 || distance >= model.range_m {
    model.nugget + model.sill
  } else {
    let ratio = distance / model.range_m
    model.nugget + model.sill * (1.5 * ratio - 0.5 * ratio * ratio * ratio)
  }
}

///|
pub fn variogram_exponential(
  distance_m : Double,
  model : VariogramModel,
) -> Double {
  if model.range_m <= 0.0 {
    model.nugget + model.sill
  } else {
    model.nugget +
    model.sill * (1.0 - @math.exp(-3.0 * distance_m.abs() / model.range_m))
  }
}

///|
pub fn variogram_gaussian(
  distance_m : Double,
  model : VariogramModel,
) -> Double {
  if model.range_m <= 0.0 {
    model.nugget + model.sill
  } else {
    let ratio = distance_m.abs() / model.range_m
    model.nugget + model.sill * (1.0 - @math.exp(-3.0 * ratio * ratio))
  }
}

///|
pub fn variogram_model_from_bins(
  bins : ArrayView[VariogramBin],
  range_m : Double,
) -> VariogramModel {
  let mut minimum = 1.0e308
  let mut maximum = 0.0
  for bin in bins {
    if bin.pairs > 0 && is_finite(bin.semivariance) {
      minimum = Double::min(minimum, bin.semivariance)
      maximum = Double::max(maximum, bin.semivariance)
    }
  }
  let baseline = if minimum == 1.0e308 { 0.0 } else { minimum }
  {
    nugget: if minimum == 1.0e308 {
      0.0
    } else {
      minimum
    },
    sill: Double::max(0.0, maximum - baseline),
    range_m: Double::max(0.0, range_m),
  }
}

///|
pub fn kriging_like_weights(
  samples : ArrayView[FieldSample],
  point : GeoPoint,
  model : VariogramModel,
) -> Array[Double] {
  let raw : Array[Double] = []
  let mut total = 0.0
  for sample in samples {
    let distance = haversine_distance_m(sample.point, point)
    let weight = if distance == 0.0 {
      1.0e12
    } else {
      1.0 / Double::max(1.0e-12, variogram_spherical(distance, model))
    }
    raw.push(weight)
    total += weight
  }
  if total == 0.0 {
    raw.map(fn(_) { 0.0 })
  } else {
    raw.map(fn(weight) { weight / total })
  }
}

///|
pub fn kriging_like_estimate(
  samples : ArrayView[FieldSample],
  point : GeoPoint,
  model : VariogramModel,
) -> PointEstimate {
  let weights = kriging_like_weights(samples, point, model)
  let mut value = 0.0
  let mut neighbours = 0
  for i in 0.. 0.0 {
      value += samples[i].value * weights[i]
      neighbours += 1
    }
  }
  { point, value: if neighbours == 0 { 0.0 / 0.0 } else { value }, neighbours }
}

///|
pub fn leave_one_out_idw(
  samples : ArrayView[FieldSample],
  power : Double,
) -> CrossValidationReport {
  let predictions : Array[Double] = []
  let residuals : Array[Double] = []
  let safe_power = Double::max(0.1, power)
  for i in 0.. 0.0 {
          let weight = 1.0 / @math.pow(distance, safe_power)
          weighted += weight * samples[j].value
          total_weight += weight
        }
      }
    }
    let prediction = if total_weight == 0.0 {
      0.0 / 0.0
    } else {
      weighted / total_weight
    }
    predictions.push(prediction)
    residuals.push(samples[i].value - prediction)
  }
  {
    samples: samples.length(),
    predictions,
    residuals,
    rmse: rms(residuals),
    mean_absolute_error: mean_absolute_deviation(residuals),
  }
}

///|
pub fn local_linear_trend(
  samples : ArrayView[FieldSample],
  origin : GeoPoint,
) -> LeastSquaresResult {
  let x : Array[Double] = []
  let y : Array[Double] = []
  for sample in samples {
    if is_finite(sample.value) {
      let offset = local_enu_delta_m(origin, sample.point)
      x.push(offset.x)
      y.push(sample.value)
    }
  }
  polynomial_fit(x, y, 1)
}

///|
pub fn detrend_spatial_samples(
  samples : ArrayView[FieldSample],
  origin : GeoPoint,
) -> Array[FieldSample] {
  let fit = local_linear_trend(samples, origin)
  let output : Array[FieldSample] = []
  for sample in samples {
    let offset = local_enu_delta_m(origin, sample.point)
    let trend = polynomial_evaluate(fit.coefficients, offset.x)
    output.push({ ..sample, value: sample.value - trend })
  }
  output
}

///|
pub fn moran_spatial_autocorrelation(
  samples : ArrayView[FieldSample],
  radius_m : Double,
) -> Double {
  let usable : Array[FieldSample] = []
  for sample in samples {
    if is_finite(sample.value) {
      usable.push(sample)
    }
  }
  if usable.length() < 2 {
    0.0
  } else {
    let values = usable.map(fn(sample) { sample.value })
    let mean = sum_finite(values) / Double::from_int(values.length())
    let mut numerator = 0.0
    let mut denominator = 0.0
    let mut weights = 0.0
    for i in 0.. Array[Double] {
  let values : Array[Double] = []
  for sample in samples {
    if is_finite(sample.value) {
      values.push(sample.value)
    }
  }
  let output : Array[Double] = []
  for fraction in fractions {
    output.push(percentile(values, fraction))
  }
  output
}

///|
pub fn sample_value_range(samples : ArrayView[FieldSample]) -> NumericRange {
  let values : Array[Double] = []
  for sample in samples {
    if is_finite(sample.value) {
      values.push(sample.value)
    }
  }
  let summary = summary(values)
  { minimum: summary.minimum, maximum: summary.maximum }
}