///|
/// A reproducible pseudo-random generator for uncertainty studies.
pub(all) struct DeterministicRng {
  state : Int64
} derive(Debug, Eq, ToJson)

///|
pub fn DeterministicRng::new(seed : Int64) -> DeterministicRng {
  { state: if seed == 0 { 88172645463393265L } else { seed.abs() } }
}

///|
pub fn DeterministicRng::next(
  self : DeterministicRng,
) -> (DeterministicRng, Double) {
  let next = self.state * 2862933555777941757L + 3037000493L
  let positive = next.abs()
  let value = (positive % 9007199254740991L).to_double() / 9007199254740991.0
  ({ state: next }, value)
}

///|
pub fn DeterministicRng::uniform(
  self : DeterministicRng,
  lower~ : Double,
  upper~ : Double,
) -> (DeterministicRng, Double) raise @domain.ExtractionError {
  @domain.require_finite("uniform lower", lower)
  @domain.require_finite("uniform upper", upper)
  if upper < lower {
    raise @domain.ExtractionError::InvalidArgument(
      message="uniform upper bound is below lower bound",
    )
  }
  let (next, unit) = self.next()
  (next, lower + (upper - lower) * unit)
}

///|
pub fn DeterministicRng::integer(
  self : DeterministicRng,
  lower~ : Int,
  upper~ : Int,
) -> (DeterministicRng, Int) raise @domain.ExtractionError {
  if upper < lower {
    raise @domain.ExtractionError::InvalidArgument(
      message="integer upper bound is below lower bound",
    )
  }
  let (next, unit) = self.next()
  let width = upper - lower + 1
  (next, lower + (unit * width.to_double()).floor().to_int().min(width - 1))
}

///|
pub fn DeterministicRng::normal_approx(
  self : DeterministicRng,
  mean~ : Double,
  standard_deviation~ : Double,
) -> (DeterministicRng, Double) raise @domain.ExtractionError {
  @domain.require_finite("normal mean", mean)
  @domain.require_non_negative("normal standard deviation", standard_deviation)
  let mut generator = self
  let mut sum = 0.0
  for _ in 0..<12 {
    let (next, value) = generator.next()
    generator = next
    sum += value
  }
  (generator, mean + (sum - 6.0) * standard_deviation)
}

///|
pub fn DeterministicRng::triangular(
  self : DeterministicRng,
  lower~ : Double,
  mode~ : Double,
  upper~ : Double,
) -> (DeterministicRng, Double) raise @domain.ExtractionError {
  @domain.require_finite("triangular lower", lower)
  @domain.require_finite("triangular mode", mode)
  @domain.require_finite("triangular upper", upper)
  if mode < lower || mode > upper || upper <= lower {
    raise @domain.ExtractionError::InvalidArgument(
      message="triangular parameters are not ordered",
    )
  }
  let (next, unit) = self.next()
  let split = (mode - lower) / (upper - lower)
  if unit <= split {
    (next, lower + (unit * (upper - lower) * (mode - lower)).max(0.0).sqrt())
  } else {
    (
      next,
      upper - ((1.0 - unit) * (upper - lower) * (upper - mode)).max(0.0).sqrt(),
    )
  }
}

///|
pub(all) enum Distribution {
  Uniform(lower~ : Double, upper~ : Double)
  Triangular(lower~ : Double, mode~ : Double, upper~ : Double)
  Normal(mean~ : Double, standard_deviation~ : Double)
} derive(Debug, Eq, ToJson)

///|
fn sample_distribution(
  generator : DeterministicRng,
  distribution : Distribution,
) -> (DeterministicRng, Double) raise @domain.ExtractionError {
  match distribution {
    Uniform(lower~, upper~) => generator.uniform(lower~, upper~)
    Triangular(lower~, mode~, upper~) =>
      generator.triangular(lower~, mode~, upper~)
    Normal(mean~, standard_deviation~) =>
      generator.normal_approx(mean~, standard_deviation~)
  }
}

///|
pub(all) struct UncertainInput {
  name : String
  distribution : Distribution
} derive(Debug, Eq, ToJson)

///|
pub fn UncertainInput::new(
  name~ : String,
  distribution~ : Distribution,
) -> UncertainInput raise @domain.ExtractionError {
  if name.trim().length() == 0 {
    raise @domain.ExtractionError::InvalidArgument(
      message="uncertain input name cannot be blank",
    )
  }
  match distribution {
    Uniform(lower~, upper~) => {
      @domain.require_finite("uniform lower", lower)
      @domain.require_finite("uniform upper", upper)
      if upper < lower {
        raise @domain.ExtractionError::InvalidArgument(
          message="uniform bounds are not ordered",
        )
      }
    }
    Triangular(lower~, mode~, upper~) => {
      @domain.require_finite("triangular lower", lower)
      @domain.require_finite("triangular mode", mode)
      @domain.require_finite("triangular upper", upper)
      if mode < lower || mode > upper || upper <= lower {
        raise @domain.ExtractionError::InvalidArgument(
          message="triangular bounds are not ordered",
        )
      }
    }
    Normal(mean~, standard_deviation~) => {
      @domain.require_finite("normal mean", mean)
      @domain.require_non_negative(
        "normal standard deviation", standard_deviation,
      )
    }
  }
  { name, distribution }
}

///|
pub(all) struct SampleSummary {
  count : Int
  minimum : Double
  q05 : Double
  q25 : Double
  median : Double
  q75 : Double
  q95 : Double
  maximum : Double
  mean : Double
  standard_deviation : Double
} derive(Debug, Eq, ToJson)

///|
pub fn summarize_samples(
  values : ArrayView[Double],
) -> SampleSummary raise @domain.ExtractionError {
  if values.length() == 0 {
    raise @domain.ExtractionError::InsufficientData(required=1, actual=0)
  }
  for value in values {
    @domain.require_finite("sample output", value)
  }
  let sorted = values.to_owned()
  sorted.sort()
  let mut sum = 0.0
  for value in sorted {
    sum += value
  }
  let mean = sum / sorted.length().to_double()
  let variance = sorted.fold(init=0.0, (total, value) => {
      total + (value - mean) * (value - mean)
    }) /
    sorted.length().to_double()
  {
    count: sorted.length(),
    minimum: sorted[0],
    q05: uncertainty_quantile(sorted, 0.05),
    q25: uncertainty_quantile(sorted, 0.25),
    median: uncertainty_quantile(sorted, 0.5),
    q75: uncertainty_quantile(sorted, 0.75),
    q95: uncertainty_quantile(sorted, 0.95),
    maximum: sorted[sorted.length() - 1],
    mean,
    standard_deviation: variance.sqrt(),
  }
}

///|
fn uncertainty_quantile(
  values : ArrayView[Double],
  probability : Double,
) -> Double {
  let position = probability * (values.length() - 1).to_double()
  let lower = position.floor().to_int()
  let upper = position.ceil().to_int()
  if lower == upper {
    values[lower]
  } else {
    values[lower] +
    (position - lower.to_double()) * (values[upper] - values[lower])
  }
}

///|
pub(all) struct MonteCarloResult {
  seed : Int64
  outputs : Array[Double]
  summary : SampleSummary
  input_names : Array[String]
  input_samples : Array[Array[Double]]
} derive(Debug, Eq, ToJson)

///|
/// Run a deterministic Monte Carlo study. The callback receives one value per input.
pub fn monte_carlo(
  inputs : ArrayView[UncertainInput],
  sample_count~ : Int,
  seed~ : Int64,
  evaluate : (ArrayView[Double]) -> Double raise @domain.ExtractionError,
) -> MonteCarloResult raise @domain.ExtractionError {
  if inputs.length() == 0 {
    raise @domain.ExtractionError::InsufficientData(required=1, actual=0)
  }
  if sample_count < 1 {
    raise @domain.ExtractionError::InsufficientData(
      required=1,
      actual=sample_count,
    )
  }
  let names = [ for input in inputs => input.name ]
  let all_inputs : Array[Array[Double]] = [ for _ in 0.. [] ]
  let outputs : Array[Double] = []
  let mut generator = DeterministicRng::new(seed)
  for _ in 0.. Array[Double] raise @domain.ExtractionError {
  if values.length() == 0 {
    raise @domain.ExtractionError::InsufficientData(required=1, actual=0)
  }
  if replicates < 1 {
    raise @domain.ExtractionError::InsufficientData(
      required=1,
      actual=replicates,
    )
  }
  for value in values {
    @domain.require_finite("bootstrap value", value)
  }
  let mut generator = DeterministicRng::new(seed)
  let means : Array[Double] = []
  for _ in 0.. ConfidenceInterval raise @domain.ExtractionError {
  if confidence_level <= 0.0 || confidence_level >= 1.0 {
    raise @domain.ExtractionError::InvalidArgument(
      message="confidence level must be between zero and one",
    )
  }
  let means = bootstrap_mean(values, replicates~, seed~)
  let sorted = means.copy()
  sorted.sort()
  let alpha = (1.0 - confidence_level) / 2.0
  let estimate = values.fold(init=0.0, (sum, value) => sum + value) /
    values.length().to_double()
  {
    estimate,
    lower: uncertainty_quantile(sorted, alpha),
    upper: uncertainty_quantile(sorted, 1.0 - alpha),
    confidence_level,
    sample_count: values.length(),
  }
}

///|
/// Estimate rank-free standardized sensitivity from paired input/output samples.
pub fn standardized_sensitivity(
  inputs : ArrayView[Double],
  outputs : ArrayView[Double],
) -> Double raise @domain.ExtractionError {
  if inputs.length() != outputs.length() || inputs.length() < 2 {
    raise @domain.ExtractionError::InsufficientData(
      required=2,
      actual=inputs.length(),
    )
  }
  let input_mean = inputs.fold(init=0.0, (sum, value) => sum + value) /
    inputs.length().to_double()
  let output_mean = outputs.fold(init=0.0, (sum, value) => sum + value) /
    outputs.length().to_double()
  let mut covariance = 0.0
  let mut input_variance = 0.0
  let mut output_variance = 0.0
  for index in 0.. Double raise @domain.ExtractionError {
  if index < 0 || base < 2 {
    raise @domain.ExtractionError::InvalidArgument(
      message="radical inverse requires non-negative index and base >= 2",
    )
  }
  if index == 0 {
    0.0
  } else {
    let mut remaining = index
    let mut denominator = 1.0
    let mut result = 0.0
    while remaining > 0 {
      denominator *= base.to_double()
      let digit = remaining % base
      result += digit.to_double() / denominator
      remaining /= base
    }
    result
  }
}

///|
pub fn quasi_random_grid(
  ranges : ArrayView[Range],
  count~ : Int,
) -> Array[Array[Double]] raise @domain.ExtractionError {
  if ranges.length() == 0 {
    raise @domain.ExtractionError::InsufficientData(required=1, actual=0)
  }
  if count < 1 {
    raise @domain.ExtractionError::InsufficientData(required=1, actual=count)
  }
  [
    for index in 0.. {
      [
        for dimension, range in ranges => {
          let unit = radical_inverse(index + 1, 2 + dimension)
          range.lower + unit * (range.upper - range.lower)
        }
      ]
    }
  ]
}

///|
pub(all) struct RecoveryUncertainty {
  nominal : Double
  summary : SampleSummary
  lower_recovery : Double
  upper_recovery : Double
  failure_probability : Double
} derive(Debug, Eq, ToJson)

///|
/// Propagate uncertain feed, solvent and distribution coefficient through one stage.
pub fn single_stage_recovery_uncertainty(
  feed_flow : UncertainInput,
  feed_concentration : UncertainInput,
  solvent_flow : UncertainInput,
  distribution : UncertainInput,
  sample_count~ : Int,
  seed~ : Int64,
  target_recovery~ : Double,
) -> RecoveryUncertainty raise @domain.ExtractionError {
  let inputs = [feed_flow, feed_concentration, solvent_flow, distribution]
  let study = monte_carlo(inputs, sample_count~, seed~, values => {
    let result = @liquid.single_stage_constant_k(
      feed_flow=values[0],
      feed_concentration=values[1],
      solvent_flow=values[2],
      distribution=values[3],
    )
    result.recovery
  })
  let mut failures = 0
  for value in study.outputs {
    if value < target_recovery {
      failures += 1
    }
  }
  let nominal_inputs = [
    match feed_flow.distribution {
      Uniform(lower~, upper~) => (lower + upper) / 2.0
      _ => 0.0
    },
    match feed_concentration.distribution {
      Uniform(lower~, upper~) => (lower + upper) / 2.0
      _ => 0.0
    },
    match solvent_flow.distribution {
      Uniform(lower~, upper~) => (lower + upper) / 2.0
      _ => 0.0
    },
    match distribution.distribution {
      Uniform(lower~, upper~) => (lower + upper) / 2.0
      _ => 0.0
    },
  ]
  let nominal_result = @liquid.single_stage_constant_k(
    feed_flow=nominal_inputs[0],
    feed_concentration=nominal_inputs[1],
    solvent_flow=nominal_inputs[2],
    distribution=nominal_inputs[3],
  )
  {
    nominal: nominal_result.recovery,
    summary: study.summary,
    lower_recovery: study.summary.q05,
    upper_recovery: study.summary.q95,
    failure_probability: failures.to_double() /
    study.outputs.length().to_double(),
  }
}