///|
fn observation_x(
  observation : AntoineObservation,
  c : Double,
) -> Double raise VleError {
  let temperature_c = observation.temperature_k - 273.15
  let denominator = temperature_c + c
  if denominator == 0.0 {
    raise VleError::InvalidParameter("Antoine denominator cannot be zero")
  }
  0.0 - 1.0 / denominator
}

///|
fn observation_y(observation : AntoineObservation) -> Double raise VleError {
  if observation.temperature_k <= 0.0 || observation.pressure_bar <= 0.0 {
    raise VleError::InvalidParameter("Antoine observations must be positive")
  }
  @math.log10(observation.pressure_bar / 0.001333223684)
}

///|
fn fit_linear_candidate(
  observations : Array[AntoineObservation],
  c : Double,
) -> (Double, Double) raise VleError {
  let n = observations.length().to_double()
  let (sx, sy, sxx, sxy) = for i = 0, sx = 0.0, sy = 0.0, sxx = 0.0, sxy = 0.0; i <
                              observations.length(); {
    let x = observation_x(observations[i], c)
    let y = observation_y(observations[i])
    continue i + 1, sx + x, sy + y, sxx + x * x, sxy + x * y
  } nobreak {
    (sx, sy, sxx, sxy)
  }
  let denominator = n * sxx - sx * sx
  if abs_double(denominator) <= 1.0E-18 {
    raise VleError::InvalidParameter("Antoine observations are rank deficient")
  }
  let b = (n * sxy - sx * sy) / denominator
  let a = (sy - b * sx) / n
  (a, b)
}

///|
pub fn antoine_fit_pressure_bar(
  a : Double,
  b : Double,
  c : Double,
  temperature_k : Double,
) -> Double raise VleError {
  assert_temperature(temperature_k)
  let denominator = temperature_k - 273.15 + c
  if denominator == 0.0 {
    raise VleError::InvalidParameter("Antoine denominator cannot be zero")
  }
  @math.pow(10.0, a + b * (0.0 - 1.0 / denominator)) * 0.001333223684
}

///|
pub fn fit_antoine(
  observations : Array[AntoineObservation],
  initial_c~ : Double,
  c_span~ : Double,
  c_steps~ : Int,
) -> AntoineFitResult raise VleError {
  if observations.length() < 3 {
    raise VleError::InvalidParameter("Antoine fit needs at least three points")
  }
  if c_span <= 0.0 || c_steps <= 0 {
    raise VleError::InvalidParameter(
      "Antoine search span and steps must be positive",
    )
  }
  for i = 0; i < observations.length(); i = i + 1 {
    ignore(observation_y(observations[i]))
  }
  let step_size = 2.0 * c_span / c_steps.to_double()
  let (best_a, best_b, best_c, best_sse) = for i = 0, a = 0.0, b = 0.0, c = initial_c, sse = 1.0E300; i <=
                                              c_steps; {
    let candidate_c = initial_c - c_span + i.to_double() * step_size
    let (candidate_a, candidate_b) = fit_linear_candidate(
      observations, candidate_c,
    )
    let candidate_sse = for j = 0, value = 0.0; j < observations.length(); {
      let predicted = antoine_fit_pressure_bar(
        candidate_a,
        candidate_b,
        candidate_c,
        observations[j].temperature_k,
      )
      let error = predicted - observations[j].pressure_bar
      continue j + 1, value + error * error
    } nobreak {
      value
    }
    if candidate_sse < sse {
      continue i + 1, candidate_a, candidate_b, candidate_c, candidate_sse
    } else {
      continue i + 1, a, b, c, sse
    }
  } nobreak {
    (a, b, c, sse)
  }
  let rmse = @math.pow(best_sse / observations.length().to_double(), 0.5)
  let max_error = for i = 0, value = 0.0; i < observations.length(); {
    let error = abs_double(
      antoine_fit_pressure_bar(
        best_a,
        best_b,
        best_c,
        observations[i].temperature_k,
      ) -
      observations[i].pressure_bar,
    )
    continue i + 1, if error > value { error } else { value }
  } nobreak {
    value
  }
  AntoineFitResult::{
    a: best_a,
    b: best_b,
    c: best_c,
    rmse_bar: rmse,
    max_error_bar: max_error,
    points: observations.length(),
    converged: true,
  }
}

///|
pub fn antoine_fit_residuals(
  result : AntoineFitResult,
  observations : Array[AntoineObservation],
) -> Array[Double] raise VleError {
  [
    for observation in observations => {
      antoine_fit_pressure_bar(
        result.a,
        result.b,
        result.c,
        observation.temperature_k,
      ) -
      observation.pressure_bar
    }
  ]
}

///|
pub fn antoine_fit_mean_absolute_error(
  result : AntoineFitResult,
  observations : Array[AntoineObservation],
) -> Double raise VleError {
  let residuals = antoine_fit_residuals(result, observations)
  if residuals.length() == 0 {
    raise VleError::EmptyMixture
  }
  for i = 0, value = 0.0; i < residuals.length(); {
    continue i + 1, value + abs_double(residuals[i])
  } nobreak {
    value / residuals.length().to_double()
  }
}

///|
pub fn antoine_fit_predict_many(
  result : AntoineFitResult,
  temperatures_k : Array[Double],
) -> Array[Double] raise VleError {
  [
    for temperature in temperatures_k => {
      antoine_fit_pressure_bar(result.a, result.b, result.c, temperature)
    }
  ]
}

///|
pub fn antoine_fit_is_reasonable(
  result : AntoineFitResult,
  max_rmse_bar : Double,
) -> Bool raise VleError {
  if max_rmse_bar <= 0.0 {
    raise VleError::InvalidParameter("fit threshold must be positive")
  }
  result.converged && result.rmse_bar <= max_rmse_bar
}

///|
pub fn relative_error(
  reference : Double,
  predicted : Double,
) -> Double raise VleError {
  if reference == 0.0 {
    raise VleError::InvalidParameter("relative error reference cannot be zero")
  }
  abs_double(predicted - reference) / abs_double(reference)
}

///|
pub fn mean_absolute_error(
  references : Array[Double],
  predicted : Array[Double],
) -> Double raise VleError {
  assert_same_length(references.length(), predicted.length())
  if references.length() == 0 {
    raise VleError::EmptyMixture
  }
  for i = 0, value = 0.0; i < references.length(); {
    continue i + 1, value + abs_double(predicted[i] - references[i])
  } nobreak {
    value / references.length().to_double()
  }
}

///|
pub fn root_mean_square_error(
  references : Array[Double],
  predicted : Array[Double],
) -> Double raise VleError {
  assert_same_length(references.length(), predicted.length())
  if references.length() == 0 {
    raise VleError::EmptyMixture
  }
  @math.pow(
    for i = 0, value = 0.0; i < references.length(); {
      let error = predicted[i] - references[i]
      continue i + 1, value + error * error
    } nobreak {
      value / references.length().to_double()
    },
    0.5,
  )
}

///|
pub fn maximum_absolute_error(
  references : Array[Double],
  predicted : Array[Double],
) -> Double raise VleError {
  assert_same_length(references.length(), predicted.length())
  if references.length() == 0 {
    raise VleError::EmptyMixture
  }
  for i = 0, value = 0.0; i < references.length(); {
    let error = abs_double(predicted[i] - references[i])
    continue i + 1, if error > value { error } else { value }
  } nobreak {
    value
  }
}