///|
/// Stage-by-stage profile and sensitivity utilities for liquid extraction.
pub(all) struct CascadeProfilePoint {
  stage : Int
  raffinate_in : Double
  raffinate_out : Double
  extract_in : Double
  extract_out : Double
  local_recovery : Double
  cumulative_recovery : Double
  driving_force : Double
  balance_error : Double
} derive(Debug, Eq, ToJson)

///|
pub(all) struct CascadeProfile {
  stages : Array[CascadeProfilePoint]
  total_recovery : Double
  total_balance_error : Double
  minimum_driving_force : Double
  bottleneck_stage : Int
  monotonic : Bool
} derive(Debug, Eq, ToJson)

///|
pub(all) struct CascadeSensitivityPoint {
  parameter : String
  baseline : Double
  perturbed : Double
  absolute_change : Double
  relative_change : Double
  elasticity : Double
} derive(Debug, Eq, ToJson)

///|
pub fn build_cascade_profile(
  raffinate_inlets : ArrayView[Double],
  raffinate_outlets : ArrayView[Double],
  extract_inlets : ArrayView[Double],
  extract_outlets : ArrayView[Double],
  tolerance : Double,
) -> CascadeProfile raise @domain.ExtractionError {
  let count = raffinate_inlets.length()
  if count == 0 ||
    raffinate_outlets.length() != count ||
    extract_inlets.length() != count ||
    extract_outlets.length() != count {
    raise @domain.ExtractionError::InvalidArgument(
      message="cascade profile arrays have different lengths",
    )
  }
  let points = []
  let mut cumulative = 0.0
  let mut minimum_force = 1.0e30
  let mut bottleneck = 0
  let mut total_error = 0.0
  let mut monotonic = true
  let mut previous_out = raffinate_inlets[0]
  for index in 0.. 0 && outgoing > previous_out + tolerance {
      monotonic = false
    }
    previous_out = outgoing
    points.push({
      stage: index + 1,
      raffinate_in: incoming,
      raffinate_out: outgoing,
      extract_in,
      extract_out,
      local_recovery,
      cumulative_recovery: cumulative,
      driving_force: force,
      balance_error: balance,
    })
  }
  let total_feed = raffinate_inlets[0]
  {
    stages: points,
    total_recovery: if total_feed == 0.0 {
      0.0
    } else {
      (total_feed - raffinate_outlets[count - 1]) / total_feed
    },
    total_balance_error: total_error,
    minimum_driving_force: if minimum_force == 1.0e30 {
      0.0
    } else {
      minimum_force
    },
    bottleneck_stage: bottleneck,
    monotonic,
  }
}

///|
pub fn profile_from_stage_data(
  feed : Double,
  stage_recoveries : ArrayView[Double],
  solvent_flow : Double,
) -> CascadeProfile raise @domain.ExtractionError {
  if feed < 0.0 || solvent_flow < 0.0 || stage_recoveries.length() == 0 {
    raise @domain.ExtractionError::InvalidArgument(
      message="stage profile inputs are invalid",
    )
  }
  let raffinate_in = []
  let raffinate_out = []
  let extract_in = []
  let extract_out = []
  let mut current = feed
  let mut solvent = 0.0
  for recovery in stage_recoveries {
    if recovery < 0.0 || recovery > 1.0 {
      raise @domain.ExtractionError::InvalidFraction(
        name="stage recovery",
        value=recovery,
      )
    }
    raffinate_in.push(current)
    raffinate_out.push(current * (1.0 - recovery))
    extract_in.push(solvent)
    let transfer = current * recovery
    solvent += transfer
    extract_out.push(solvent)
    current = current * (1.0 - recovery)
  }
  build_cascade_profile(
    raffinate_in, raffinate_out, extract_in, extract_out, 1.0e-12,
  )
}

///|
pub fn profile_recoveries(profile : CascadeProfile) -> Array[Double] {
  [
    for point in profile.stages => point.local_recovery
  ]
}

///|
pub fn profile_driving_forces(profile : CascadeProfile) -> Array[Double] {
  [
    for point in profile.stages => point.driving_force
  ]
}

///|
pub fn profile_balance_errors(profile : CascadeProfile) -> Array[Double] {
  [
    for point in profile.stages => point.balance_error
  ]
}

///|
pub fn profile_stage(
  profile : CascadeProfile,
  stage : Int,
) -> CascadeProfilePoint raise @domain.ExtractionError {
  if stage < 1 || stage > profile.stages.length() {
    raise @domain.ExtractionError::OutOfRange(
      name="stage",
      value=stage.to_double(),
      lower=1.0,
      upper=profile.stages.length().to_double(),
    )
  }
  profile.stages[stage - 1]
}

///|
pub fn profile_average_recovery(profile : CascadeProfile) -> Double {
  if profile.stages.length() == 0 {
    return 0.0
  }
  let mut total = 0.0
  for point in profile.stages {
    total += point.local_recovery
  }
  total / profile.stages.length().to_double()
}

///|
pub fn profile_recovery_stddev(profile : CascadeProfile) -> Double {
  if profile.stages.length() == 0 {
    return 0.0
  }
  let mean = profile_average_recovery(profile)
  let mut total = 0.0
  for point in profile.stages {
    let difference = point.local_recovery - mean
    total += difference * difference
  }
  (total / profile.stages.length().to_double()).sqrt()
}

///|
pub fn profile_last_stage_share(profile : CascadeProfile) -> Double {
  if profile.stages.length() == 0 || profile.total_recovery == 0.0 {
    return 0.0
  }
  profile.stages[profile.stages.length() - 1].local_recovery /
  profile.total_recovery
}

///|
pub fn profile_first_stage_share(profile : CascadeProfile) -> Double {
  if profile.stages.length() == 0 || profile.total_recovery == 0.0 {
    return 0.0
  }
  profile.stages[0].local_recovery / profile.total_recovery
}

///|
pub fn profile_is_well_balanced(
  profile : CascadeProfile,
  maximum_balance_error : Double,
  minimum_force : Double,
) -> Bool {
  profile.total_balance_error <= maximum_balance_error &&
  profile.minimum_driving_force >= minimum_force
}

///|
pub fn profile_is_acceptable(
  profile : CascadeProfile,
  minimum_recovery : Double,
  maximum_error : Double,
) -> Bool {
  profile.total_recovery >= minimum_recovery &&
  profile.total_balance_error <= maximum_error &&
  profile.monotonic
}

///|
pub fn profile_sensitivity(
  baseline : CascadeProfile,
  perturbed : CascadeProfile,
  parameter : String,
) -> CascadeSensitivityPoint {
  let change = perturbed.total_recovery - baseline.total_recovery
  let relative = if baseline.total_recovery == 0.0 {
    0.0
  } else {
    change / baseline.total_recovery
  }
  {
    parameter,
    baseline: baseline.total_recovery,
    perturbed: perturbed.total_recovery,
    absolute_change: change,
    relative_change: relative,
    elasticity: relative,
  }
}

///|
pub fn profile_sensitivity_matrix(
  baseline : CascadeProfile,
  variants : ArrayView[CascadeProfile],
  names : ArrayView[String],
) -> Array[CascadeSensitivityPoint] {
  let result = []
  let count = variants.length().min(names.length())
  for index in 0.. String {
  let lines = [
    "stage,raffinate_in,raffinate_out,extract_in,extract_out,local_recovery,cumulative_recovery,driving_force,balance_error",
  ]
  for point in profile.stages {
    lines.push(
      point.stage.to_string() +
      "," +
      point.raffinate_in.to_string() +
      "," +
      point.raffinate_out.to_string() +
      "," +
      point.extract_in.to_string() +
      "," +
      point.extract_out.to_string() +
      "," +
      point.local_recovery.to_string() +
      "," +
      point.cumulative_recovery.to_string() +
      "," +
      point.driving_force.to_string() +
      "," +
      point.balance_error.to_string(),
    )
  }
  lines.join("\n")
}

///|
pub fn profile_markdown(profile : CascadeProfile) -> String {
  let lines = [
    "# Cascade profile", "", "| Stage | Raffinate in | Raffinate out | Extract in | Extract out | Recovery | Driving force | Error |",
    "|---:|---:|---:|---:|---:|---:|---:|---:|",
  ]
  for point in profile.stages {
    lines.push(
      "| " +
      point.stage.to_string() +
      " | " +
      point.raffinate_in.to_string() +
      " | " +
      point.raffinate_out.to_string() +
      " | " +
      point.extract_in.to_string() +
      " | " +
      point.extract_out.to_string() +
      " | " +
      point.local_recovery.to_string() +
      " | " +
      point.driving_force.to_string() +
      " | " +
      point.balance_error.to_string() +
      " |",
    )
  }
  lines.push("")
  lines.push("Total recovery: " + profile.total_recovery.to_string())
  lines.push("Bottleneck stage: " + profile.bottleneck_stage.to_string())
  lines.join("\n")
}

///|
pub fn profile_resample(
  profile : CascadeProfile,
  target_stages : Int,
) -> CascadeProfile raise @domain.ExtractionError {
  if target_stages < 1 {
    raise @domain.ExtractionError::InvalidArgument(
      message="target stages must be positive",
    )
  }
  if profile.stages.length() == target_stages {
    return profile
  }
  let recoveries = profile_recoveries(profile)
  let resized = Array::make(target_stages, 0.0)
  for index in 0.. Array[String] {
  let alerts = []
  for point in profile.stages {
    if point.driving_force < minimum_force {
      alerts.push("stage " + point.stage.to_string() + ": low driving force")
    }
    if point.balance_error > maximum_error {
      alerts.push("stage " + point.stage.to_string() + ": balance error")
    }
  }
  if !profile.monotonic {
    alerts.push("raffinate profile is not monotonic")
  }
  alerts
}

///|
pub fn profile_compare(
  first : CascadeProfile,
  second : CascadeProfile,
) -> String {
  let recovery_change = second.total_recovery - first.total_recovery
  let error_change = second.total_balance_error - first.total_balance_error
  "recovery_change=" +
  recovery_change.to_string() +
  ",balance_error_change=" +
  error_change.to_string()
}