///|
/// Operating-envelope utilities for translating process limits into
/// executable checks and scenario summaries.
pub(all) struct OperatingLimit {
  name : String
  lower : Double
  upper : Double
  unit : String
  critical : Bool
} derive(Debug, Eq, ToJson)

///|
pub(all) struct OperatingPoint {
  name : String
  values : Array[Double]
  feasible : Bool
  violations : Array[String]
  distance_to_limit : Double
} derive(Debug, Eq, ToJson)

///|
pub(all) struct ScenarioPoint {
  name : String
  temperature : Double
  pressure : Double
  flow : Double
  solvent_ratio : Double
  recovery : Double
  feasible : Bool
} derive(Debug, Eq, ToJson)

///|
pub(all) struct EnvelopeSummary {
  evaluated : Int
  feasible : Int
  infeasible : Int
  minimum_recovery : Double
  maximum_recovery : Double
  average_recovery : Double
  limiting_scenario : String
} derive(Debug, Eq, ToJson)

///|
pub fn OperatingLimit::new(
  name~ : String,
  lower~ : Double,
  upper~ : Double,
  unit? : String = "",
  critical? : Bool = true,
) -> OperatingLimit raise @domain.ExtractionError {
  if lower > upper {
    raise @domain.ExtractionError::InvalidArgument(
      message="operating lower limit exceeds upper limit",
    )
  }
  { name, lower, upper, unit, critical }
}

///|
pub fn OperatingLimit::contains(self : OperatingLimit, value : Double) -> Bool {
  value >= self.lower && value <= self.upper
}

///|
pub fn OperatingLimit::normalized_distance(
  self : OperatingLimit,
  value : Double,
) -> Double {
  if value < self.lower {
    self.lower - value
  } else if value > self.upper {
    value - self.upper
  } else {
    let span = self.upper - self.lower
    if span == 0.0 {
      0.0
    } else {
      let lower_distance = value - self.lower
      let upper_distance = self.upper - value
      lower_distance.min(upper_distance) / span
    }
  }
}

///|
pub fn OperatingLimit::to_text(self : OperatingLimit) -> String {
  let severity = if self.critical { "critical" } else { "advisory" }
  self.name +
  " [" +
  self.lower.to_string() +
  ", " +
  self.upper.to_string() +
  "] " +
  self.unit +
  " (" +
  severity +
  ")"
}

///|
pub fn OperatingPoint::new(
  name~ : String,
  values~ : ArrayView[Double],
  limits~ : ArrayView[OperatingLimit],
) -> OperatingPoint raise @domain.ExtractionError {
  if values.length() != limits.length() {
    raise @domain.ExtractionError::InvalidArgument(
      message="values and limits have different widths",
    )
  }
  let violations = []
  let mut distance = 1.0e30
  for index in 0.. Bool {
  for index in 0.. Array[Double] {
  [
    for index in 0.. {
      let limit = limits[index]
      if self.values[index] < limit.lower {
        self.values[index] - limit.lower
      } else if self.values[index] > limit.upper {
        limit.upper - self.values[index]
      } else {
        (self.values[index] - limit.lower).min(limit.upper - self.values[index])
      }
    }
  ]
}

///|
pub fn evaluate_grid(
  names : ArrayView[String],
  axes : ArrayView[Array[Double]],
  limits : ArrayView[OperatingLimit],
) -> Array[OperatingPoint] raise @domain.ExtractionError {
  if names.length() != axes.length() || names.length() != limits.length() {
    raise @domain.ExtractionError::InvalidArgument(
      message="grid axes do not match limits",
    )
  }
  if names.length() == 0 {
    return []
  }
  let points = []
  let indices = Array::make(names.length(), 0)
  let mut finished = false
  while !finished {
    let values = [
      for index in 0.. axes[index][indices[index]]
    ]
    points.push(OperatingPoint::new(name=names.join("/"), values~, limits~))
    let mut carry = indices.length() - 1
    while carry >= 0 {
      indices[carry] += 1
      if indices[carry] < axes[carry].length() {
        break
      }
      indices[carry] = 0
      carry -= 1
    }
    if carry < 0 {
      finished = true
    }
  }
  points
}

///|
pub fn feasible_points(
  points : ArrayView[OperatingPoint],
) -> Array[OperatingPoint] {
  let selected = []
  for point in points {
    if point.feasible {
      selected.push(point)
    }
  }
  selected
}

///|
pub fn infeasible_points(
  points : ArrayView[OperatingPoint],
) -> Array[OperatingPoint] {
  let selected = []
  for point in points {
    if !point.feasible {
      selected.push(point)
    }
  }
  selected
}

///|
pub fn worst_margin(points : ArrayView[OperatingPoint]) -> Double {
  if points.length() == 0 {
    return 0.0
  }
  let mut margin = 1.0e30
  for point in points {
    if point.distance_to_limit < margin {
      margin = point.distance_to_limit
    }
  }
  margin
}

///|
pub fn envelope_recovery(
  temperature : Double,
  pressure : Double,
  flow : Double,
  solvent_ratio : Double,
  equilibrium_factor : Double,
  stage_count : Int,
) -> Double raise @domain.ExtractionError {
  if temperature <= 0.0 ||
    pressure <= 0.0 ||
    flow <= 0.0 ||
    solvent_ratio < 0.0 ||
    equilibrium_factor <= 0.0 ||
    stage_count < 1 {
    raise @domain.ExtractionError::InvalidArgument(
      message="envelope scenario is invalid",
    )
  }
  let contact = equilibrium_factor * solvent_ratio / (1.0 + solvent_ratio)
  let stage_factor = 1.0 - 1.0 / ((1.0 + contact) * stage_count.to_double())
  let pressure_factor = 1.0 -
    0.02 * ((pressure - 1.0).abs() / pressure.max(1.0))
  let temperature_factor = 1.0 - 0.001 * (temperature - 298.15).abs().min(100.0)
  (stage_factor * pressure_factor * temperature_factor).max(0.0).min(1.0)
}

///|
pub fn make_scenario(
  name~ : String,
  temperature~ : Double,
  pressure~ : Double,
  flow~ : Double,
  solvent_ratio~ : Double,
  equilibrium_factor~ : Double,
  stage_count~ : Int,
) -> ScenarioPoint raise @domain.ExtractionError {
  let recovery = envelope_recovery(
    temperature, pressure, flow, solvent_ratio, equilibrium_factor, stage_count,
  )
  {
    name,
    temperature,
    pressure,
    flow,
    solvent_ratio,
    recovery,
    feasible: recovery >= 0.5,
  }
}

///|
pub fn summarize_scenarios(
  scenarios : ArrayView[ScenarioPoint],
) -> EnvelopeSummary raise @domain.ExtractionError {
  if scenarios.length() == 0 {
    raise @domain.ExtractionError::InsufficientData(required=1, actual=0)
  }
  let mut feasible = 0
  let mut minimum = scenarios[0].recovery
  let mut maximum = scenarios[0].recovery
  let mut total = 0.0
  let mut limiting_name = scenarios[0].name
  for scenario in scenarios {
    if scenario.feasible {
      feasible += 1
    }
    if scenario.recovery < minimum {
      minimum = scenario.recovery
      limiting_name = scenario.name
    }
    if scenario.recovery > maximum {
      maximum = scenario.recovery
    }
    total += scenario.recovery
  }
  {
    evaluated: scenarios.length(),
    feasible,
    infeasible: scenarios.length() - feasible,
    minimum_recovery: minimum,
    maximum_recovery: maximum,
    average_recovery: total / scenarios.length().to_double(),
    limiting_scenario: limiting_name,
  }
}

///|
pub fn scenario_recovery_values(
  scenarios : ArrayView[ScenarioPoint],
) -> Array[Double] {
  [
    for scenario in scenarios => scenario.recovery
  ]
}

///|
pub fn scenario_names(scenarios : ArrayView[ScenarioPoint]) -> Array[String] {
  [
    for scenario in scenarios => scenario.name
  ]
}

///|
pub fn scenario_markdown(scenarios : ArrayView[ScenarioPoint]) -> String {
  let lines = [
    "| Scenario | T | P | Flow | S/F | Recovery | Feasible |", "|---|---:|---:|---:|---:|---:|---|",
  ]
  for scenario in scenarios {
    lines.push(
      "| " +
      scenario.name +
      " | " +
      scenario.temperature.to_string() +
      " | " +
      scenario.pressure.to_string() +
      " | " +
      scenario.flow.to_string() +
      " | " +
      scenario.solvent_ratio.to_string() +
      " | " +
      scenario.recovery.to_string() +
      " | " +
      scenario.feasible.to_string() +
      " |",
    )
  }
  lines.join("\n")
}

///|
pub fn recovery_quantile(
  scenarios : ArrayView[ScenarioPoint],
  probability : Double,
) -> Double raise @domain.ExtractionError {
  if scenarios.length() == 0 || probability < 0.0 || probability > 1.0 {
    raise @domain.ExtractionError::InvalidArgument(
      message="recovery quantile input is invalid",
    )
  }
  let values = scenario_recovery_values(scenarios)
  values.sort()
  let index = @math.floor(probability * (values.length() - 1).to_double()).to_int()
  values[index]
}

///|
pub fn recovery_percentile_band(
  scenarios : ArrayView[ScenarioPoint],
  lower_probability? : Double = 0.05,
  upper_probability? : Double = 0.95,
) -> (Double, Double) raise @domain.ExtractionError {
  if lower_probability < 0.0 ||
    upper_probability > 1.0 ||
    lower_probability > upper_probability {
    raise @domain.ExtractionError::InvalidArgument(
      message="recovery percentile band is invalid",
    )
  }
  (
    recovery_quantile(scenarios, lower_probability),
    recovery_quantile(scenarios, upper_probability),
  )
}

///|
pub fn robust_scenario(
  scenarios : ArrayView[ScenarioPoint],
  minimum_recovery : Double,
) -> Bool {
  for scenario in scenarios {
    if scenario.recovery < minimum_recovery {
      return false
    }
  }
  true
}

///|
pub fn scenario_sensitivity(
  baseline : ScenarioPoint,
  perturbations : ArrayView[ScenarioPoint],
) -> Array[Double] {
  [
    for scenario in perturbations => {
      (scenario.recovery - baseline.recovery) /
      (scenario.flow - baseline.flow).abs().max(1.0e-12)
    }
  ]
}

///|
pub fn utility_limit_report(limits : ArrayView[OperatingLimit]) -> String {
  let lines = [
    "| Name | Lower | Upper | Unit | Critical |", "|---|---:|---:|---|---|",
  ]
  for limit in limits {
    lines.push(
      "| " +
      limit.name +
      " | " +
      limit.lower.to_string() +
      " | " +
      limit.upper.to_string() +
      " | " +
      limit.unit +
      " | " +
      limit.critical.to_string() +
      " |",
    )
  }
  lines.join("\n")
}

///|
pub fn limits_minimum_span(limits : ArrayView[OperatingLimit]) -> Double {
  if limits.length() == 0 {
    return 0.0
  }
  let mut span = 1.0e30
  for limit in limits {
    let current = limit.upper - limit.lower
    if current < span {
      span = current
    }
  }
  span
}

///|
pub fn point_is_near_limit(point : OperatingPoint, threshold : Double) -> Bool {
  point.distance_to_limit <= threshold
}

///|
pub fn near_limit_points(
  points : ArrayView[OperatingPoint],
  threshold : Double,
) -> Array[OperatingPoint] {
  let selected = []
  for point in points {
    if point_is_near_limit(point, threshold) {
      selected.push(point)
    }
  }
  selected
}

///|
pub fn operating_point_csv(point : OperatingPoint) -> String {
  point.name +
  "," +
  point.values.map(value => value.to_string()).join(",") +
  "," +
  point.feasible.to_string() +
  "," +
  point.distance_to_limit.to_string()
}

///|
pub fn operating_points_csv(points : ArrayView[OperatingPoint]) -> String {
  let lines = ["name,values,feasible,distance_to_limit"]
  for point in points {
    lines.push(operating_point_csv(point))
  }
  lines.join("\n")
}