///|
/// Evidence collected while comparing a candidate model against observations.
pub struct ModelScore {
  name : String
  rmse : Double
  mae : Double
  complexity : Int
  consistency : Double
} derive(Debug)

///|
pub fn ModelScore::new(
  name : String,
  rmse : Double,
  mae : Double,
  complexity : Int,
  consistency : Double,
) -> ModelScore {
  {
    name,
    rmse: if rmse < 0.0 {
      0.0
    } else {
      rmse
    },
    mae: if mae < 0.0 {
      0.0
    } else {
      mae
    },
    complexity: if complexity < 0 {
      0
    } else {
      complexity
    },
    consistency: if consistency < 0.0 {
      0.0
    } else {
      consistency
    },
  }
}

///|
pub fn ModelScore::name(self : ModelScore) -> String {
  self.name
}

///|
pub fn ModelScore::rmse(self : ModelScore) -> Double {
  self.rmse
}

///|
pub fn ModelScore::mae(self : ModelScore) -> Double {
  self.mae
}

///|
pub fn ModelScore::complexity(self : ModelScore) -> Int {
  self.complexity
}

///|
pub fn ModelScore::consistency(self : ModelScore) -> Double {
  self.consistency
}

///|
pub fn ModelScore::objective(
  self : ModelScore,
  complexity_penalty : Double,
) -> Double {
  self.rmse +
  complexity_penalty * self.complexity.to_double() +
  (1.0 - self.consistency)
}

///|
pub struct ModelSelector {
  scores : Array[ModelScore]
  penalty : Double
}

///|
pub fn ModelSelector::new(penalty : Double) -> ModelSelector {
  { scores: [], penalty: if penalty < 0.0 { 0.0 } else { penalty } }
}

///|
pub fn ModelSelector::add(self : ModelSelector, score : ModelScore) -> Unit {
  self.scores.push(score)
}

///|
pub fn ModelSelector::scores(self : ModelSelector) -> Array[ModelScore] {
  self.scores.copy()
}

///|
pub fn ModelSelector::length(self : ModelSelector) -> Int {
  self.scores.length()
}

///|
pub fn ModelSelector::penalty(self : ModelSelector) -> Double {
  self.penalty
}

///|
pub fn ModelSelector::best(self : ModelSelector) -> ModelScore? {
  if self.scores.length() == 0 {
    return None
  }
  let mut index = 0
  let mut best_objective = self.scores[0].objective(self.penalty)
  for i in 1.. Array[ModelScore] {
  let result = self.scores.copy()
  for i in 1.. 0 &&
          result[index - 1].objective(self.penalty) > value_objective {
      result[index] = result[index - 1]
      index = index - 1
    }
    result[index] = value
  }
  result
}

///|
pub fn ModelSelector::clear(self : ModelSelector) -> Unit {
  self.scores.clear()
}

///|
/// Calculate common model-selection scores from residuals.
pub fn score_model(
  name : String,
  residuals : Array[Array[Double]],
  parameter_count : Int,
  consistency : Double,
) -> ModelScore {
  let zero = ModelScore::new(name, 0.0, 0.0, parameter_count, consistency)
  if residuals.length() == 0 {
    return zero
  }
  let mut squared = 0.0
  let mut absolute = 0.0
  let mut samples = 0
  for residual in residuals {
    for value in residual {
      if !value.is_nan() && !value.is_inf() {
        squared = squared + value * value
        absolute = absolute + value.abs()
        samples = samples + 1
      }
    }
  }
  if samples == 0 {
    return zero
  }
  ModelScore::new(
    name,
    (squared / samples.to_double()).sqrt(),
    absolute / samples.to_double(),
    parameter_count,
    consistency,
  )
}

///|
pub fn compare_model_scores(
  left : ModelScore,
  right : ModelScore,
  penalty : Double,
) -> Int {
  let left_objective = left.objective(penalty)
  let right_objective = right.objective(penalty)
  if left_objective < right_objective {
    -1
  } else if left_objective > right_objective {
    1
  } else {
    0
  }
}

///|
/// A reusable collection of linear models for a fixed state dimension.
pub struct ModelRegistry {
  dimension : Int
  models : Array[LinearModel]
  names : Array[String]
}

///|
pub fn ModelRegistry::new(dimension : Int) -> ModelRegistry {
  {
    dimension: if dimension < 0 {
      0
    } else {
      dimension
    },
    models: [],
    names: [],
  }
}

///|
pub fn ModelRegistry::add(
  self : ModelRegistry,
  name : String,
  model : LinearModel,
) -> Bool {
  if model.state_dimension() != self.dimension {
    return false
  }
  self.names.push(name)
  self.models.push(model)
  true
}

///|
pub fn ModelRegistry::length(self : ModelRegistry) -> Int {
  self.models.length()
}

///|
pub fn ModelRegistry::names(self : ModelRegistry) -> Array[String] {
  self.names.copy()
}

///|
pub fn ModelRegistry::model(self : ModelRegistry, index : Int) -> LinearModel? {
  self.models.get(index)
}

///|
pub fn ModelRegistry::name(self : ModelRegistry, index : Int) -> String {
  match self.names.get(index) {
    None => ""
    Some(value) => value
  }
}

///|
pub fn ModelRegistry::clear(self : ModelRegistry) -> Unit {
  self.models.clear()
  self.names.clear()
}

///|
pub fn model_complexity(model : LinearModel) -> Int {
  model.state_dimension() * model.state_dimension() +
  model.measurement_dimension() * model.state_dimension()
}

///|
pub fn model_stability_score(model : LinearModel) -> Double {
  let eigenvalues = model.transition().jacobi_eigenvalues(40)
  if eigenvalues.length() == 0 {
    return 0.0
  }
  let mut maximum = 0.0
  for value in eigenvalues {
    if value.abs() > maximum {
      maximum = value.abs()
    }
  }
  1.0 / (1.0 + maximum)
}

///|
pub fn model_is_well_formed(model : LinearModel) -> Bool {
  model.transition().rows() == model.state_dimension() &&
  model.transition().is_square() &&
  model.process_noise().rows() == model.state_dimension() &&
  model.observation().cols() == model.state_dimension() &&
  model.measurement_noise().rows() == model.measurement_dimension() &&
  covariance_is_psd(model.process_noise(), 0.001) &&
  covariance_is_psd(model.measurement_noise(), 0.001)
}