///|
fn solver_sign_change(a : Double, b : Double) -> Bool {
  (a <= 0.0 && b >= 0.0) || (a >= 0.0 && b <= 0.0)
}

///|
pub fn solve_bisection(
  low~ : Double,
  high~ : Double,
  options~ : SolverOptions,
  f : (Double) -> Double raise VleError,
) -> SolverReport raise VleError {
  if high <= low {
    raise VleError::InvalidRange(low~, high~)
  }
  let f_low = f(low)
  let f_high = f(high)
  if f_low == 0.0 {
    return SolverReport::new(
      root=low,
      residual=0.0,
      iterations=0,
      converged=true,
    )
  }
  if f_high == 0.0 {
    return SolverReport::new(
      root=high,
      residual=0.0,
      iterations=0,
      converged=true,
    )
  }
  if !solver_sign_change(f_low, f_high) {
    raise VleError::SolverDidNotBracket(low~, high~)
  }
  for i = 0, lo = low, hi = high, flo = f_low; i < options.max_iterations; {
    let mid = (lo + hi) / 2.0
    let fm = f(mid)
    if abs_double(fm) <= options.tolerance ||
      abs_double(hi - lo) <= options.tolerance {
      break SolverReport::new(
        root=mid,
        residual=abs_double(fm),
        iterations=i + 1,
        converged=true,
      )
    }
    if solver_sign_change(flo, fm) {
      continue i + 1, lo, mid, flo
    } else {
      continue i + 1, mid, hi, fm
    }
  } nobreak {
    raise VleError::SolverDidNotConverge(iterations=options.max_iterations)
  }
}

///|
pub fn solve_newton(
  initial~ : Double,
  low~ : Double,
  high~ : Double,
  options~ : SolverOptions,
  f : (Double) -> Double raise VleError,
  derivative : (Double) -> Double raise VleError,
) -> SolverReport raise VleError {
  if high <= low || initial < low || initial > high {
    raise VleError::InvalidRange(low~, high~)
  }
  let f_low = f(low)
  let f_high = f(high)
  if f_low == 0.0 {
    return SolverReport::new(
      root=low,
      residual=0.0,
      iterations=0,
      converged=true,
    )
  }
  if f_high == 0.0 {
    return SolverReport::new(
      root=high,
      residual=0.0,
      iterations=0,
      converged=true,
    )
  }
  if !solver_sign_change(f_low, f_high) {
    raise VleError::SolverDidNotBracket(low~, high~)
  }
  let initial_value = f(initial)
  for i = 0, x = initial, fx = initial_value, lo = low, hi = high, flo = f_low; i <
     options.max_iterations; {
    if abs_double(fx) <= options.tolerance {
      break SolverReport::new(
        root=x,
        residual=abs_double(fx),
        iterations=i,
        converged=true,
      )
    }
    let slope = derivative(x)
    let trial = if slope == 0.0 { (lo + hi) / 2.0 } else { x - fx / slope }
    let next = if trial <= lo || trial >= hi { (lo + hi) / 2.0 } else { trial }
    let f_next = f(next)
    if solver_sign_change(flo, f_next) {
      continue i + 1, next, f_next, lo, next, flo
    } else {
      continue i + 1, next, f_next, next, hi, f_next
    }
  } nobreak {
    raise VleError::SolverDidNotConverge(iterations=options.max_iterations)
  }
}

///|
pub fn vector_distance(
  a : Array[Double],
  b : Array[Double],
) -> Double raise VleError {
  assert_same_length(a.length(), b.length())
  @math.pow(
    for i = 0, acc = 0.0; i < a.length(); {
      continue i + 1, acc + (a[i] - b[i]) * (a[i] - b[i])
    } nobreak {
      acc
    },
    0.5,
  )
}

///|
pub fn weighted_rms(
  expected : Array[Double],
  actual : Array[Double],
) -> Double raise VleError {
  assert_same_length(expected.length(), actual.length())
  if expected.length() == 0 {
    raise VleError::EmptyMixture
  }
  @math.pow(
    for i = 0, acc = 0.0; i < expected.length(); {
      let scale = clamp_positive(abs_double(expected[i]))
      let error = (actual[i] - expected[i]) / scale
      continue i + 1, acc + error * error
    } nobreak {
      acc / expected.length().to_double()
    },
    0.5,
  )
}