///|
/// Reproducible design-of-experiments and robust statistics for process studies.
pub(all) struct Range {
  lower : Double
  upper : Double
} derive(Debug, Eq, ToJson)

///|
pub fn Range::new(
  lower~ : Double,
  upper~ : Double,
) -> Range raise @domain.ExtractionError {
  @domain.require_finite("range lower", lower)
  @domain.require_finite("range upper", upper)
  if upper < lower {
    raise @domain.ExtractionError::InvalidArgument(
      message="range upper must not be below lower",
    )
  }
  { lower, upper }
}

///|
pub fn linspace(
  range : Range,
  count : Int,
) -> Array[Double] raise @domain.ExtractionError {
  if count < 1 {
    raise @domain.ExtractionError::InsufficientData(required=1, actual=count)
  }
  if count == 1 {
    [range.lower]
  } else {
    [
      for index in 0.. {
        range.lower +
        (range.upper - range.lower) *
        index.to_double() /
        (count - 1).to_double()
      }
    ]
  }
}

///|
pub fn quantile(
  values : ArrayView[Double],
  probability : Double,
) -> Double raise @domain.ExtractionError {
  if values.length() == 0 {
    raise @domain.ExtractionError::InsufficientData(required=1, actual=0)
  }
  if probability < 0.0 || probability > 1.0 {
    raise @domain.ExtractionError::InvalidArgument(
      message="quantile probability must be in [0,1]",
    )
  }
  let sorted = values.to_owned()
  sorted.sort()
  let position = probability * (sorted.length() - 1).to_double()
  let lower = position.floor().to_int()
  let upper = position.ceil().to_int()
  if lower == upper {
    sorted[lower]
  } else {
    sorted[lower] +
    (position - lower.to_double()) * (sorted[upper] - sorted[lower])
  }
}

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

///|
pub fn summarize(
  values : ArrayView[Double],
) -> Statistics raise @domain.ExtractionError {
  if values.length() == 0 {
    raise @domain.ExtractionError::InsufficientData(required=1, actual=0)
  }
  for value in values {
    @domain.require_finite("sample", 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 mut variance = 0.0
  for value in sorted {
    variance += (value - mean) * (value - mean)
  }
  {
    count: sorted.length(),
    minimum: sorted[0],
    q25: quantile(sorted, 0.25),
    median: quantile(sorted, 0.5),
    q75: quantile(sorted, 0.75),
    maximum: sorted[sorted.length() - 1],
    mean,
    standard_deviation: (variance / sorted.length().to_double()).sqrt(),
  }
}

///|
pub fn latin_hypercube(
  ranges : ArrayView[Range],
  samples : Int,
) -> Array[Array[Double]] raise @domain.ExtractionError {
  if ranges.length() == 0 {
    raise @domain.ExtractionError::InsufficientData(required=1, actual=0)
  }
  if samples < 1 {
    raise @domain.ExtractionError::InsufficientData(required=1, actual=samples)
  }
  [
    for sample in 0.. {
      [
        for range in ranges => {
          range.lower +
          (range.upper - range.lower) *
          (sample.to_double() + 0.5) /
          samples.to_double()
        }
      ]
    }
  ]
}

///|
pub fn cartesian_grid(
  ranges : ArrayView[Range],
  counts : ArrayView[Int],
) -> Array[Array[Double]] raise @domain.ExtractionError {
  if ranges.length() != counts.length() {
    raise @domain.ExtractionError::InvalidArgument(
      message="range and count dimensions differ",
    )
  }
  let result = [[]]
  for dimension, range in ranges {
    let count = counts[dimension]
    if count < 1 {
      raise @domain.ExtractionError::InvalidArgument(
        message="grid counts must be positive",
      )
    }
    let next = []
    for prefix in result {
      for value in linspace(range, count) {
        let point = prefix.copy()
        point.push(value)
        next.push(point)
      }
    }
    result.clear()
    for point in next {
      result.push(point)
    }
  }
  result
}