///|
/// A response-surface observation for two process variables.
pub(all) struct SurfaceObservation {
  first : Double
  second : Double
  response : Double
  weight : Double
} derive(Debug, Eq, ToJson)

///|
pub fn SurfaceObservation::new(
  first~ : Double,
  second~ : Double,
  response~ : Double,
  weight? : Double = 1.0,
) -> SurfaceObservation raise @domain.ExtractionError {
  @domain.require_finite("surface first variable", first)
  @domain.require_finite("surface second variable", second)
  @domain.require_finite("surface response", response)
  @domain.require_positive("surface weight", weight)
  { first, second, response, weight }
}

///|
pub(all) struct ResponseSurfaceModel {
  coefficients : Array[Double]
  observations : Int
  residual_sum_squares : Double
  r_squared : Double
  first_center : Double
  second_center : Double
  first_scale : Double
  second_scale : Double
} derive(Debug, Eq, ToJson)

///|
fn surface_basis(first : Double, second : Double) -> Array[Double] {
  [1.0, first, second, first * first, first * second, second * second]
}

///|
fn solve_surface_system(
  matrix : Array[Array[Double]],
  rhs : Array[Double],
) -> Array[Double] raise @domain.ExtractionError {
  let n = rhs.length()
  for pivot in 0.. magnitude {
        best = row
        magnitude = matrix[row][pivot].abs()
      }
    }
    if magnitude <= 1.0e-14 {
      raise @domain.ExtractionError::InvalidArgument(
        message="surface matrix is singular",
      )
    }
    if best != pivot {
      let row = matrix[pivot]
      matrix[pivot] = matrix[best]
      matrix[best] = row
      let value = rhs[pivot]
      rhs[pivot] = rhs[best]
      rhs[best] = value
    }
    for row in (pivot + 1).. ResponseSurfaceModel raise @domain.ExtractionError {
  if observations.length() < 6 {
    raise @domain.ExtractionError::InsufficientData(
      required=6,
      actual=observations.length(),
    )
  }
  let first_center = observations.fold(init=0.0, (sum, item) => sum + item.first) /
    observations.length().to_double()
  let second_center = observations.fold(init=0.0, (sum, item) => {
      sum + item.second
    }) /
    observations.length().to_double()
  let first_scale = observations
    .fold(init=0.0, (maximum, item) => {
      maximum.max((item.first - first_center).abs())
    })
    .max(1.0)
  let second_scale = observations
    .fold(init=0.0, (maximum, item) => {
      maximum.max((item.second - second_center).abs())
    })
    .max(1.0)
  let normal = [ for _ in 0..<6 => [ for _ in 0..<6 => 0.0 ] ]
  let rhs = Array::make(6, 0.0)
  let mut weight_total = 0.0
  let mut response_mean = 0.0
  for item in observations {
    let first = (item.first - first_center) / first_scale
    let second = (item.second - second_center) / second_scale
    let basis = surface_basis(first, second)
    for row in 0..<6 {
      rhs[row] += item.weight * item.response * basis[row]
      for column in 0..<6 {
        normal[row][column] += item.weight * basis[row] * basis[column]
      }
    }
    weight_total += item.weight
    response_mean += item.weight * item.response
  }
  let coefficients = solve_surface_system(normal, rhs)
  let mean = response_mean / weight_total
  let mut residual_sum = 0.0
  let mut total_sum = 0.0
  for item in observations {
    let first = (item.first - first_center) / first_scale
    let second = (item.second - second_center) / second_scale
    let prediction = coefficients
      .zip(surface_basis(first, second))
      .fold(init=0.0, (sum, pair) => sum + pair.0 * pair.1)
    residual_sum += item.weight *
      (item.response - prediction) *
      (item.response - prediction)
    total_sum += item.weight * (item.response - mean) * (item.response - mean)
  }
  {
    coefficients,
    observations: observations.length(),
    residual_sum_squares: residual_sum,
    r_squared: if total_sum <= 1.0e-20 {
      1.0
    } else {
      1.0 - residual_sum / total_sum
    },
    first_center,
    second_center,
    first_scale,
    second_scale,
  }
}

///|
pub fn ResponseSurfaceModel::predict(
  self : ResponseSurfaceModel,
  first : Double,
  second : Double,
) -> Double {
  let normalized_first = (first - self.first_center) / self.first_scale
  let normalized_second = (second - self.second_center) / self.second_scale
  let basis = surface_basis(normalized_first, normalized_second)
  self.coefficients
  .zip(basis)
  .fold(init=0.0, (sum, pair) => sum + pair.0 * pair.1)
}

///|
pub fn ResponseSurfaceModel::coefficient(
  self : ResponseSurfaceModel,
  index : Int,
) -> Double raise @domain.ExtractionError {
  match self.coefficients.get(index) {
    Some(value) => value
    None =>
      raise @domain.ExtractionError::InvalidArgument(
        message="surface coefficient index out of range",
      )
  }
}

///|
pub fn ResponseSurfaceModel::gradient(
  self : ResponseSurfaceModel,
  first : Double,
  second : Double,
) -> (Double, Double) {
  let x = (first - self.first_center) / self.first_scale
  let y = (second - self.second_center) / self.second_scale
  let first_gradient = (
      self.coefficients[1] +
      2.0 * self.coefficients[3] * x +
      self.coefficients[4] * y
    ) /
    self.first_scale
  let second_gradient = (
      self.coefficients[2] +
      2.0 * self.coefficients[5] * y +
      self.coefficients[4] * x
    ) /
    self.second_scale
  (first_gradient, second_gradient)
}

///|
pub(all) struct SurfaceGridPoint {
  first : Double
  second : Double
  response : Double
  first_gradient : Double
  second_gradient : Double
} derive(Debug, Eq, ToJson)

///|
pub fn ResponseSurfaceModel::grid(
  self : ResponseSurfaceModel,
  first_lower~ : Double,
  first_upper~ : Double,
  second_lower~ : Double,
  second_upper~ : Double,
  count~ : Int,
) -> Array[SurfaceGridPoint] raise @domain.ExtractionError {
  if count < 2 {
    raise @domain.ExtractionError::InsufficientData(required=2, actual=count)
  }
  if first_upper <= first_lower || second_upper <= second_lower {
    raise @domain.ExtractionError::InvalidArgument(
      message="surface grid bounds are invalid",
    )
  }
  let points : Array[SurfaceGridPoint] = []
  for first_index in 0.. Array[SurfaceResidual] {
  [
    for item in observations => {
      let predicted = model.predict(item.first, item.second)
      {
        first: item.first,
        second: item.second,
        observed: item.response,
        predicted,
        residual: item.response - predicted,
      }
    }
  ]
}

///|
pub fn surface_rmse(
  model : ResponseSurfaceModel,
  observations : ArrayView[SurfaceObservation],
) -> Double raise @domain.ExtractionError {
  if observations.length() == 0 {
    raise @domain.ExtractionError::InsufficientData(required=1, actual=0)
  }
  let residuals = response_surface_residuals(model, observations)
  (residuals.fold(init=0.0, (sum, item) => sum + item.residual * item.residual) /
  residuals.length().to_double()).sqrt()
}

///|
pub fn select_surface_optimum(
  model : ResponseSurfaceModel,
  observations : ArrayView[SurfaceObservation],
  maximize? : Bool = true,
) -> SurfaceObservation raise @domain.ExtractionError {
  if observations.length() == 0 {
    raise @domain.ExtractionError::InsufficientData(required=1, actual=0)
  }
  let mut best = observations[0]
  for item in observations {
    let best_value = model.predict(best.first, best.second)
    let value = model.predict(item.first, item.second)
    if (if maximize { value > best_value } else { value < best_value }) {
      best = item
    }
  }
  best
}