///|
pub fn is_finite(value : Double) -> Bool {
  value == value && value.abs() < 1.7976931348623157e308
}

///|
pub fn approx_equal(
  a : Double,
  b : Double,
  tolerance? : Double = 0.0000001,
) -> Bool {
  let scale = 1.0 + a.abs().max(b.abs())
  (a - b).abs() <= tolerance * scale
}

///|
pub fn clamp_double(value : Double, lower : Double, upper : Double) -> Double {
  if lower > upper {
    upper
  } else {
    value.max(lower).min(upper)
  }
}

///|
pub fn lerp(a : Double, b : Double, t : Double) -> Double {
  a + (b - a) * t
}

///|
pub fn Vector2::new(x~ : Double, y~ : Double) -> Vector2 {
  { x, y }
}

///|
pub fn Vector2::add(a : Vector2, b : Vector2) -> Vector2 {
  { x: a.x + b.x, y: a.y + b.y }
}

///|
pub fn Vector2::sub(a : Vector2, b : Vector2) -> Vector2 {
  { x: a.x - b.x, y: a.y - b.y }
}

///|
pub fn Vector2::scale(v : Vector2, factor : Double) -> Vector2 {
  { x: v.x * factor, y: v.y * factor }
}

///|
pub fn Vector2::dot(a : Vector2, b : Vector2) -> Double {
  a.x * b.x + a.y * b.y
}

///|
pub fn Vector2::cross(a : Vector2, b : Vector2) -> Double {
  a.x * b.y - a.y * b.x
}

///|
pub fn Vector2::length(v : Vector2) -> Double {
  (v.x * v.x + v.y * v.y).sqrt()
}

///|
pub fn Vector2::normalized(v : Vector2) -> Vector2 {
  let length = v.length()
  if length <= 0.0000000001 {
    { x: 0.0, y: 0.0 }
  } else {
    v.scale(1.0 / length)
  }
}

///|
pub fn Vector2::distance(a : Vector2, b : Vector2) -> Double {
  a.sub(b).length()
}

///|
pub fn Vector2::angle(v : Vector2) -> Double {
  if v.x > 0.0 {
    @math.atan(v.y / v.x)
  } else if v.x < 0.0 && v.y >= 0.0 {
    @math.atan(v.y / v.x) + 3.141592653589793
  } else if v.x < 0.0 {
    @math.atan(v.y / v.x) - 3.141592653589793
  } else if v.y > 0.0 {
    1.5707963267948966
  } else if v.y < 0.0 {
    -1.5707963267948966
  } else {
    0.0
  }
}

///|
pub fn Matrix2::new(
  m00~ : Double,
  m01~ : Double,
  m10~ : Double,
  m11~ : Double,
) -> Matrix2 {
  { m00, m01, m10, m11 }
}

///|
pub fn Matrix2::identity() -> Matrix2 {
  { m00: 1.0, m01: 0.0, m10: 0.0, m11: 1.0 }
}

///|
pub fn Matrix2::determinant(matrix : Matrix2) -> Double {
  matrix.m00 * matrix.m11 - matrix.m01 * matrix.m10
}

///|
pub fn Matrix2::trace(matrix : Matrix2) -> Double {
  matrix.m00 + matrix.m11
}

///|
pub fn Matrix2::add(a : Matrix2, b : Matrix2) -> Matrix2 {
  {
    m00: a.m00 + b.m00,
    m01: a.m01 + b.m01,
    m10: a.m10 + b.m10,
    m11: a.m11 + b.m11,
  }
}

///|
pub fn Matrix2::scale(matrix : Matrix2, factor : Double) -> Matrix2 {
  {
    m00: matrix.m00 * factor,
    m01: matrix.m01 * factor,
    m10: matrix.m10 * factor,
    m11: matrix.m11 * factor,
  }
}

///|
pub fn Matrix2::multiply(a : Matrix2, b : Matrix2) -> Matrix2 {
  {
    m00: a.m00 * b.m00 + a.m01 * b.m10,
    m01: a.m00 * b.m01 + a.m01 * b.m11,
    m10: a.m10 * b.m00 + a.m11 * b.m10,
    m11: a.m10 * b.m01 + a.m11 * b.m11,
  }
}

///|
pub fn Matrix2::apply(matrix : Matrix2, vector : Vector2) -> Vector2 {
  {
    x: matrix.m00 * vector.x + matrix.m01 * vector.y,
    y: matrix.m10 * vector.x + matrix.m11 * vector.y,
  }
}

///|
pub fn Matrix2::inverse(matrix : Matrix2) -> Matrix2 {
  let determinant = matrix.determinant()
  if determinant.abs() <= 0.0000000001 {
    Matrix2::identity()
  } else {
    {
      m00: matrix.m11 / determinant,
      m01: -matrix.m01 / determinant,
      m10: -matrix.m10 / determinant,
      m11: matrix.m00 / determinant,
    }
  }
}

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

///|
pub fn mean_value(values : Array[Double]) -> Double {
  if values.is_empty() {
    0.0
  } else {
    sum_values(values) / Double::from_int(values.length())
  }
}

///|
pub fn minimum_value(values : Array[Double]) -> Double {
  if values.is_empty() {
    0.0
  } else {
    let mut result = values[0]
    for value in values[1:] {
      result = result.min(value)
    }
    result
  }
}

///|
pub fn maximum_value(values : Array[Double]) -> Double {
  if values.is_empty() {
    0.0
  } else {
    let mut result = values[0]
    for value in values[1:] {
      result = result.max(value)
    }
    result
  }
}

///|
pub fn variance_value(values : Array[Double]) -> Double {
  if values.length() < 2 {
    0.0
  } else {
    let average = mean_value(values)
    let mut squared = 0.0
    for value in values {
      let difference = value - average
      squared += difference * difference
    }
    squared / Double::from_int(values.length() - 1)
  }
}

///|
fn insertion_sorted(values : Array[Double]) -> Array[Double] {
  let result : Array[Double] = []
  for value in values {
    let mut index = 0
    while index < result.length() && result[index] <= value {
      index += 1
    }
    result.insert(index, value)
  }
  result
}

///|
fn array_slice(
  values : Array[Double],
  first : Int,
  last : Int,
) -> Array[Double] {
  let result : Array[Double] = []
  for i in first.max(0).. Double {
  let sorted = insertion_sorted(values)
  if sorted.is_empty() {
    0.0
  } else if sorted.length() % 2 == 1 {
    sorted[sorted.length() / 2]
  } else {
    let right = sorted.length() / 2
    (sorted[right - 1] + sorted[right]) / 2.0
  }
}

///|
pub fn percentile_value(values : Array[Double], q : Double) -> Double {
  let sorted = insertion_sorted(values)
  if sorted.is_empty() {
    0.0
  } else {
    let position = clamp_double(q, 0.0, 1.0) *
      Double::from_int(sorted.length() - 1)
    let lower = position.floor().to_int()
    let upper = position.ceil().to_int()
    if lower == upper {
      sorted[lower]
    } else {
      let fraction = position - Double::from_int(lower)
      lerp(sorted[lower], sorted[upper], fraction)
    }
  }
}

///|
pub fn describe_samples(values : Array[Double]) -> StatsSummary {
  let ordered = insertion_sorted(values)
  let count = values.length()
  let average = mean_value(values)
  let variance = variance_value(values)
  {
    count,
    sum: sum_values(values),
    mean: average,
    variance,
    standard_deviation: variance.sqrt(),
    minimum: minimum_value(values),
    maximum: maximum_value(values),
    median: median_value(values),
    ordered,
  }
}

///|
pub fn robust_mean(
  values : Array[Double],
  trim_fraction? : Double = 0.1,
) -> Double {
  let ordered = insertion_sorted(values)
  if ordered.is_empty() {
    0.0
  } else {
    let fraction = clamp_double(trim_fraction, 0.0, 0.49)
    let trim = (Double::from_int(ordered.length()) * fraction).floor().to_int()
    let first = trim
    let last = ordered.length() - trim
    if first >= last {
      median_value(ordered)
    } else {
      mean_value(array_slice(ordered, first, last))
    }
  }
}

///|
pub fn mean_absolute_deviation(values : Array[Double]) -> Double {
  let center = median_value(values)
  if values.is_empty() {
    0.0
  } else {
    let deviations = values.map(fn(value) { (value - center).abs() })
    median_value(deviations)
  }
}

///|
pub fn weighted_mean(values : Array[Double], weights : Array[Double]) -> Double {
  if values.is_empty() || values.length() != weights.length() {
    0.0
  } else {
    let mut numerator = 0.0
    let mut denominator = 0.0
    for i in 0.. Double {
  if a.length() < 2 || a.length() != b.length() {
    0.0
  } else {
    let mean_a = mean_value(a)
    let mean_b = mean_value(b)
    let mut total = 0.0
    for i in 0.. Double {
  let denominator = variance_value(a).sqrt() * variance_value(b).sqrt()
  if denominator <= 0.0000000001 {
    0.0
  } else {
    covariance_values(a, b) / denominator
  }
}

///|
pub fn round_to(value : Double, decimals : Int) -> Double {
  let mut factor = 1.0
  if decimals > 0 {
    for _ in 0.. Array[Double] {
  let total = sum_values(values)
  if total.abs() <= 0.0000000001 {
    values.map(fn(_) { 0.0 })
  } else {
    values.map(fn(value) { value / total })
  }
}

///|
pub fn moving_average(values : Array[Double], window : Int) -> Array[Double] {
  let result : Array[Double] = []
  if window <= 0 {
    result
  } else {
    for i in 0.. Array[Double] {
  let result : Array[Double] = []
  if values.is_empty() {
    result
  } else {
    let coefficient = clamp_double(alpha, 0.0, 1.0)
    let mut previous = values[0]
    result.push(previous)
    for value in values[1:] {
      previous = coefficient * value + (1.0 - coefficient) * previous
      result.push(previous)
    }
    result
  }
}