///|
/// Run deterministic workloads used by the local benchmark harness.
fn main {
  let matrix_checksum = benchmark_matrix()
  let interpolation_checksum = benchmark_interpolation()
  let statistics_checksum = benchmark_statistics()
  let calibration_checksum = benchmark_calibration()
  println("matrix_checksum=" + matrix_checksum.to_string())
  println("interpolation_checksum=" + interpolation_checksum.to_string())
  println("statistics_checksum=" + statistics_checksum.to_string())
  println("calibration_checksum=" + calibration_checksum.to_string())
}

///|
fn benchmark_matrix() -> Double {
  let mut checksum = 0.0
  for iteration = 0; iteration < 32; iteration = iteration + 1 {
    let matrix = make_matrix(10, iteration.to_double())
    let transpose = matrix.transpose()
    let product = matrix.mul(transpose).unwrap()
    let solved = product
      .solve_matrix(@matrix.Matrix::identity(10).unwrap())
      .unwrap()
    checksum = checksum + solved.trace().unwrap() + product.frobenius_norm()
  }
  checksum
}

///|
fn make_matrix(size : Int, offset : Double) -> @matrix.Matrix {
  let values = Array::make(size * size, 0.0)
  for row = 0; row < size; row = row + 1 {
    for col = 0; col < size; col = col + 1 {
      let diagonal = if row == col { 12.0 } else { 0.05 }
      values[row * size + col] = diagonal +
        (row + col).to_double() * 0.01 +
        offset * 0.001
    }
  }
  @matrix.Matrix::from_array(size, size, values).unwrap()
}

///|
fn benchmark_interpolation() -> Double {
  let x = [0.0, 1.0, 2.0, 3.0, 4.0, 5.0]
  let y = [0.0, 1.2, 3.8, 8.9, 16.1, 25.2]
  let spline = @interpolation.CubicSpline::fit(x, y).unwrap()
  let barycentric = @interpolation.BarycentricInterpolator::fit(x, y).unwrap()
  let grid = @interpolation.RegularGrid2D::fit(
    [0.0, 1.0, 2.0, 3.0],
    [0.0, 1.0, 2.0],
    [0.0, 1.0, 4.0, 9.0, 1.0, 2.0, 5.0, 10.0, 2.0, 3.0, 6.0, 11.0],
  ).unwrap()
  let mut checksum = 0.0
  for iteration = 0; iteration < 4000; iteration = iteration + 1 {
    let query = (iteration % 500).to_double() / 100.0
    checksum = checksum + spline.eval(query).unwrap()
    checksum = checksum + barycentric.eval(query)
    checksum = checksum + grid.eval(query / 2.0, query / 3.0).unwrap()
  }
  checksum
}

///|
fn benchmark_statistics() -> Double {
  let mut checksum = 0.0
  let observations = Array::make(32, [0.0, 0.0, 0.0])
  for row = 0; row < observations.length(); row = row + 1 {
    let value = row.to_double()
    observations[row] = [value, value * 0.5 + 1.0, value * value * 0.01]
  }
  for iteration = 0; iteration < 1000; iteration = iteration + 1 {
    checksum = checksum +
      @statistics.summarize([1.0, 2.0, 3.0, 4.0, 5.0]).unwrap().iqr
    checksum = checksum +
      @statistics.covariance([1.0, 2.0, 3.0], [2.0, 4.0, 6.0]).unwrap()
    checksum = checksum +
      @statistics.covariance_matrix(observations).unwrap().trace().unwrap()
  }
  checksum
}

///|
fn benchmark_calibration() -> Double {
  let x = [0.0, 1.0, 2.0, 3.0, 4.0, 5.0]
  let y = [1.0, 2.1, 5.2, 10.3, 17.4, 26.5]
  let curve = @calibration.CalibrationCurve::fit(x, y).unwrap()
  let polynomial = @calibration.PolynomialCalibration::fit(x, y, 2).unwrap()
  let mut checksum = 0.0
  for iteration = 0; iteration < 500; iteration = iteration + 1 {
    let query = (iteration % 500).to_double() / 100.0
    checksum = checksum + curve.predict(query).unwrap()
    checksum = checksum + polynomial.predict(query)
  }
  checksum + curve.leave_one_out_rmse().unwrap()
}