///|
/// Dot product with a safe length check.
pub fn vector_dot(left : Array[Double], right : Array[Double]) -> Double {
  if left.length() != right.length() {
    return 0.0
  }
  let mut total = 0.0
  for i in 0.. Double {
  let mut total = 0.0
  for value in values {
    total = total + value
  }
  total
}

///|
pub fn vector_mean(values : Array[Double]) -> Double {
  if values.length() == 0 {
    return 0.0
  }
  vector_sum(values) / values.length().to_double()
}

///|
pub fn vector_l1_norm(values : Array[Double]) -> Double {
  let mut total = 0.0
  for value in values {
    total = total + value.abs()
  }
  total
}

///|
pub fn vector_l2_norm(values : Array[Double]) -> Double {
  vector_dot(values, values).sqrt()
}

///|
pub fn vector_linf_norm(values : Array[Double]) -> Double {
  let mut result = 0.0
  for value in values {
    let magnitude = value.abs()
    if magnitude > result {
      result = magnitude
    }
  }
  result
}

///|
pub fn vector_distance(left : Array[Double], right : Array[Double]) -> Double {
  if left.length() != right.length() {
    return 0.0
  }
  let mut total = 0.0
  for i in 0.. Array[Double] {
  if left.length() != right.length() {
    return []
  }
  Array::makei(left.length(), i => left[i] + right[i])
}

///|
pub fn vector_sub(left : Array[Double], right : Array[Double]) -> Array[Double] {
  if left.length() != right.length() {
    return []
  }
  Array::makei(left.length(), i => left[i] - right[i])
}

///|
pub fn vector_scale(values : Array[Double], factor : Double) -> Array[Double] {
  values.map(value => value * factor)
}

///|
pub fn vector_axpy(
  alpha : Double,
  x : Array[Double],
  y : Array[Double],
) -> Array[Double] {
  if x.length() != y.length() {
    return []
  }
  Array::makei(x.length(), i => alpha * x[i] + y[i])
}

///|
pub fn vector_hadamard(
  left : Array[Double],
  right : Array[Double],
) -> Array[Double] {
  if left.length() != right.length() {
    return []
  }
  Array::makei(left.length(), i => left[i] * right[i])
}

///|
/// Normalize a vector.  Zero vectors stay zero instead of producing NaNs.
pub fn vector_normalize(values : Array[Double]) -> Array[Double] {
  let norm = vector_l2_norm(values)
  if norm <= 0.000000000001 {
    return values.map(_ => 0.0)
  }
  vector_scale(values, 1.0 / norm)
}

///|
pub fn vector_clamp(
  values : Array[Double],
  lower : Double,
  upper : Double,
) -> Array[Double] {
  let low = if lower <= upper { lower } else { upper }
  let high = if lower <= upper { upper } else { lower }
  values.map(value => {
    if value < low {
      low
    } else if value > high {
      high
    } else {
      value
    }
  })
}

///|
pub fn vector_lerp(
  left : Array[Double],
  right : Array[Double],
  amount : Double,
) -> Array[Double] {
  if left.length() != right.length() {
    return []
  }
  Array::makei(left.length(), i => left[i] + amount * (right[i] - left[i]))
}

///|
pub fn vector_weighted_mean(
  values : Array[Double],
  weights : Array[Double],
) -> Double {
  if values.length() != weights.length() || values.length() == 0 {
    return 0.0
  }
  let mut numerator = 0.0
  let mut denominator = 0.0
  for i in 0.. Double {
  if values.length() < 2 {
    return 0.0
  }
  let mean = vector_mean(values)
  let mut total = 0.0
  for value in values {
    let difference = value - mean
    total = total + difference * difference
  }
  total / (values.length() - 1).to_double()
}

///|
pub fn vector_weighted_variance(
  values : Array[Double],
  weights : Array[Double],
) -> Double {
  if values.length() != weights.length() || values.length() == 0 {
    return 0.0
  }
  let mean = vector_weighted_mean(values, weights)
  let mut numerator = 0.0
  let mut denominator = 0.0
  for i in 0.. Double {
  if actual.length() != expected.length() || actual.length() == 0 {
    return 0.0
  }
  let mut total = 0.0
  for i in 0.. Double {
  if actual.length() != expected.length() || actual.length() == 0 {
    return 0.0
  }
  let mut total = 0.0
  for i in 0.. Double {
  if values.length() == 0 {
    return 0.0
  }
  let sorted = values.copy()
  sorted.sort()
  let middle = sorted.length() / 2
  if sorted.length() % 2 == 1 {
    sorted[middle]
  } else {
    (sorted[middle - 1] + sorted[middle]) * 0.5
  }
}

///|
/// Linear-interpolated quantile in the closed interval [0, 1].
pub fn vector_quantile(values : Array[Double], probability : Double) -> Double {
  if values.length() == 0 {
    return 0.0
  }
  let sorted = values.copy()
  sorted.sort()
  let p = if probability < 0.0 {
    0.0
  } else if probability > 1.0 {
    1.0
  } else {
    probability
  }
  let position = p * (sorted.length() - 1).to_double()
  let lower = position.to_int()
  let upper = if lower + 1 >= sorted.length() { lower } else { lower + 1 }
  sorted[lower] +
  (position - lower.to_double()) * (sorted[upper] - sorted[lower])
}

///|
pub fn vector_mad(values : Array[Double]) -> Double {
  let median = vector_median(values)
  let deviations = values.map(value => (value - median).abs())
  vector_median(deviations)
}

///|
pub fn vector_is_finite(values : Array[Double]) -> Bool {
  for value in values {
    if value.is_nan() || value.is_inf() {
      return false
    }
  }
  true
}

///|
pub fn vector_replace_non_finite(
  values : Array[Double],
  fallback : Double,
) -> Array[Double] {
  values.map(value => {
    if value.is_nan() || value.is_inf() {
      fallback
    } else {
      value
    }
  })
}

///|
pub fn vector_mean_center(values : Array[Double]) -> Array[Double] {
  let mean = vector_mean(values)
  values.map(value => value - mean)
}

///|
pub fn vector_covariance(samples : Array[Array[Double]]) -> Matrix {
  if samples.length() < 2 {
    return Matrix::zeros(0, 0)
  }
  let width = samples[0].length()
  for sample in samples {
    if sample.length() != width {
      return Matrix::zeros(0, 0)
    }
  }
  let mean = Array::makei(width, j => {
    let mut total = 0.0
    for sample in samples {
      total = total + sample[j]
    }
    total / samples.length().to_double()
  })
  let result = Matrix::zeros(width, width)
  for sample in samples {
    let centered = vector_sub(sample, mean)
    let contribution = Matrix::outer(centered, centered)
    for i in 0.. ignore
      }
    }
  }
  result.scale(1.0 / (samples.length() - 1).to_double())
}

///|
pub fn vector_project(
  value : Array[Double],
  basis : Array[Double],
) -> Array[Double] {
  let denominator = vector_dot(basis, basis)
  if denominator <= 0.000000000001 {
    return basis.map(_ => 0.0)
  }
  vector_scale(basis, vector_dot(value, basis) / denominator)
}

///|
pub fn vector_reject(
  value : Array[Double],
  basis : Array[Double],
) -> Array[Double] {
  vector_sub(value, vector_project(value, basis))
}

///|
pub fn vector_wrap(value : Double, period : Double) -> Double {
  if period <= 0.0 {
    return value
  }
  let mut result = value
  let half = period * 0.5
  while result > half {
    result = result - period
  }
  while result <= -half {
    result = result + period
  }
  result
}

///|
pub fn angle_difference(target : Double, source : Double) -> Double {
  vector_wrap(target - source, 6.283185307179586)
}

///|
pub fn vector_all_close(
  left : Array[Double],
  right : Array[Double],
  tolerance : Double,
) -> Bool {
  if left.length() != right.length() {
    return false
  }
  for i in 0.. tolerance {
      return false
    }
  }
  true
}