///|
fn feature_statistics(
  x : Array[Array[Double]],
) -> (Array[Double], Array[Double]) {
  if x.length() == 0 {
    return ([], [])
  }
  let width = x[0].length()
  let means = Array::make(width, 0.0)
  for row in x {
    for j in 0.. Array[Double] {
  let result = Array::new(capacity=row.length())
  for j in 0.. Double {
  let standardized = standardized_row(
    row,
    model.feature_means,
    model.feature_scales,
  )
  model.intercept + dot(standardized, model.coefficients)
}

///|
/// Fits a binary logistic model with standardized features and L2 regularization.
pub fn fit_logistic_regression(
  covariates : Array[Array[Double]],
  treatment : Array[Bool],
  learning_rate? : Double = 0.15,
  max_iterations? : Int = 600,
  l2? : Double = 1.0e-6,
) -> ModelFit {
  let n = covariates.length()
  let width = if n == 0 { 0 } else { covariates[0].length() }
  let (means, scales) = feature_statistics(covariates)
  let coefficients = Array::make(width, 0.0)
  let mut intercept = 0.0
  let mut converged = false
  let mut iterations = 0
  for iteration in 0.. max_step {
        max_step = step.abs()
      }
    }
    iterations = iteration + 1
    if max_step < 1.0e-8 {
      converged = true
      break
    }
  }
  {
    coefficients,
    intercept,
    iterations,
    converged,
    loss: logistic_log_loss(
      covariates, treatment, coefficients, intercept, means, scales, l2,
    ),
    feature_means: means,
    feature_scales: scales,
  }
}

///|
fn logistic_log_loss(
  covariates : Array[Array[Double]],
  treatment : Array[Bool],
  coefficients : Array[Double],
  intercept : Double,
  means : Array[Double],
  scales : Array[Double],
  l2 : Double,
) -> Double {
  if covariates.length() == 0 {
    return 0.0
  }
  let mut loss = 0.0
  for i in 0.. Array[Double] {
  let result = Array::new(capacity=covariates.length())
  for row in covariates {
    result.push(safe_probability(sigmoid(model_linear_predict(row, model))))
  }
  result
}

///|
pub fn model_predict(model : ModelFit, row : Array[Double]) -> Double {
  model_linear_predict(row, model)
}

///|
/// Fits an ordinary or ridge-regularized linear outcome model.
pub fn fit_linear_outcome_model(
  covariates : Array[Array[Double]],
  outcomes : Array[Double],
  ridge? : Double = 1.0e-8,
) -> ModelFit {
  let n = covariates.length()
  let width = if n == 0 { 0 } else { covariates[0].length() }
  let (means, scales) = feature_statistics(covariates)
  let design = matrix_zeros(n, width + 1)
  for i in 0.. 0 {
    normal[0][0] -= ridge
  }
  let rhs = matrix_vector_multiply(transposed, outcomes)
  let solution = solve_linear_system(normal, rhs)
  let coefficients = Array::new(capacity=width)
  let mut intercept = 0.0
  if solution.length() == width + 1 {
    intercept = solution[0]
    for j in 0.. Array[Double] {
  let result = Array::new(capacity=covariates.length())
  for row in covariates {
    result.push(model_linear_predict(row, model))
  }
  result
}

///|
pub fn residuals(
  model : ModelFit,
  covariates : Array[Array[Double]],
  outcomes : Array[Double],
) -> Array[Double] {
  let predictions = predict_outcomes(model, covariates)
  let result = Array::new(capacity=outcomes.length())
  for i in 0.. Double {
  if predictions.length() == 0 || predictions.length() != treatment.length() {
    return 0.0
  }
  let mut correct = 0
  for i in 0..= threshold
    if predicted == treatment[i] {
      correct += 1
    }
  }
  correct.to_double() / predictions.length().to_double()
}