///|
pub struct BootstrapInterval {
  estimate : Double
  lower : Double
  upper : Double
  confidence : Double
  replicates : Int
}

///|
pub struct DeterministicRng {
  mut state : Int
}

///|
pub fn DeterministicRng::new(seed : Int) -> DeterministicRng {
  let normalized = if seed == 0 { 1 } else { seed }
  { state: normalized }
}

///|
pub fn DeterministicRng::next_int(self : DeterministicRng, upper : Int) -> Int {
  if upper <= 0 {
    abort("upper bound must be positive")
  }
  self.state = self.state * 1103515245 + 12345
  let candidate = if self.state < 0 { -self.state } else { self.state }
  candidate % upper
}

///|
pub fn DeterministicRng::next_unit(self : DeterministicRng) -> Double {
  self.next_int(1000000).to_double() / 1000000.0
}

///|
pub fn resample_with_replacement(
  data : Array[Double],
  sample_size : Int,
  seed : Int,
) -> Array[Double] {
  if sample_size < 0 {
    abort("sample_size must not be negative")
  }
  if data.length() == 0 {
    return []
  }
  let rng = DeterministicRng::new(seed)
  let result = []
  for index = 0; index < sample_size; index = index + 1 {
    result.push(data[rng.next_int(data.length())])
  }
  result
}

///|
pub fn bootstrap_replicates(
  data : Array[Double],
  replicates : Int,
  sample_size? : Int = -1,
  seed? : Int = 12345,
) -> Array[Array[Double]] {
  if replicates < 0 {
    abort("replicates must not be negative")
  }
  let actual_size = if sample_size < 0 { data.length() } else { sample_size }
  if actual_size < 0 {
    abort("sample_size must not be negative")
  }
  let result = []
  let rng = DeterministicRng::new(seed)
  for replicate = 0; replicate < replicates; replicate = replicate + 1 {
    let sample = []
    if data.length() > 0 {
      for index = 0; index < actual_size; index = index + 1 {
        sample.push(data[rng.next_int(data.length())])
      }
    }
    result.push(sample)
  }
  result
}

///|
fn confidence_z(confidence : Double) -> Double {
  validate_confidence(confidence)
  if confidence >= 0.999 {
    3.291
  } else if confidence >= 0.99 {
    2.576
  } else if confidence >= 0.95 {
    1.96
  } else if confidence >= 0.90 {
    1.645
  } else if confidence >= 0.80 {
    1.282
  } else {
    1.0
  }
}

///|
pub fn bootstrap_mean_interval(
  data : Array[Double],
  replicates : Int,
  confidence? : Double = 0.95,
  seed? : Int = 12345,
) -> BootstrapInterval {
  if data.length() == 0 || replicates <= 0 {
    return { estimate: 0.0, lower: 0.0, upper: 0.0, confidence, replicates: 0 }
  }
  let samples = bootstrap_replicates(data, replicates, seed~)
  let estimates = []
  for sample in samples {
    estimates.push(mean(sample))
  }
  let estimate = mean(data)
  let spread = sample_stddev(estimates)
  let z = confidence_z(confidence)
  {
    estimate,
    lower: estimate - z * spread,
    upper: estimate + z * spread,
    confidence,
    replicates,
  }
}

///|
pub fn bootstrap_median_interval(
  data : Array[Double],
  replicates : Int,
  confidence? : Double = 0.95,
  seed? : Int = 12345,
) -> BootstrapInterval {
  if data.length() == 0 || replicates <= 0 {
    return { estimate: 0.0, lower: 0.0, upper: 0.0, confidence, replicates: 0 }
  }
  let samples = bootstrap_replicates(data, replicates, seed~)
  let estimates = []
  for sample in samples {
    estimates.push(median(sample))
  }
  let estimate = median(data)
  let lower = quantile(estimates, (1.0 - confidence) / 2.0)
  let upper = quantile(estimates, (1.0 + confidence) / 2.0)
  { estimate, lower, upper, confidence, replicates }
}

///|
pub fn bootstrap_trimmed_mean_interval(
  data : Array[Double],
  trim_percent : Double,
  replicates : Int,
  confidence? : Double = 0.95,
  seed? : Int = 12345,
) -> BootstrapInterval {
  if data.length() == 0 || replicates <= 0 {
    return { estimate: 0.0, lower: 0.0, upper: 0.0, confidence, replicates: 0 }
  }
  let samples = bootstrap_replicates(data, replicates, seed~)
  let estimates = []
  for sample in samples {
    estimates.push(trimmed_mean(sample, trim_percent))
  }
  let estimate = trimmed_mean(data, trim_percent)
  let lower = quantile(estimates, (1.0 - confidence) / 2.0)
  let upper = quantile(estimates, (1.0 + confidence) / 2.0)
  { estimate, lower, upper, confidence, replicates }
}

///|
pub fn bootstrap_mad_interval(
  data : Array[Double],
  replicates : Int,
  confidence? : Double = 0.95,
  seed? : Int = 12345,
) -> BootstrapInterval {
  if data.length() == 0 || replicates <= 0 {
    return { estimate: 0.0, lower: 0.0, upper: 0.0, confidence, replicates: 0 }
  }
  let samples = bootstrap_replicates(data, replicates, seed~)
  let estimates = []
  for sample in samples {
    estimates.push(mad(sample))
  }
  let estimate = mad(data)
  let lower = quantile(estimates, (1.0 - confidence) / 2.0)
  let upper = quantile(estimates, (1.0 + confidence) / 2.0)
  { estimate, lower, upper, confidence, replicates }
}

///|
pub fn permutation_difference(
  first : Array[Double],
  second : Array[Double],
  replicates : Int,
  seed? : Int = 12345,
) -> Array[Double] {
  if replicates < 0 {
    abort("replicates must not be negative")
  }
  if first.length() == 0 || second.length() == 0 {
    return []
  }
  let pooled = []
  for value in first {
    pooled.push(value)
  }
  for value in second {
    pooled.push(value)
  }
  let rng = DeterministicRng::new(seed)
  let result = []
  for replicate = 0; replicate < replicates; replicate = replicate + 1 {
    let selected_first = []
    let selected_second = []
    for index = 0; index < pooled.length(); index = index + 1 {
      if rng.next_unit() <
        first.length().to_double() / pooled.length().to_double() {
        selected_first.push(pooled[index])
      } else {
        selected_second.push(pooled[index])
      }
    }
    if selected_first.length() == 0 || selected_second.length() == 0 {
      result.push(0.0)
    } else {
      result.push(mean(selected_first) - mean(selected_second))
    }
  }
  result
}

///|
pub fn bootstrap_difference_interval(
  first : Array[Double],
  second : Array[Double],
  replicates : Int,
  confidence? : Double = 0.95,
  seed? : Int = 12345,
) -> BootstrapInterval {
  if first.length() == 0 || second.length() == 0 || replicates <= 0 {
    return { estimate: 0.0, lower: 0.0, upper: 0.0, confidence, replicates: 0 }
  }
  let first_samples = bootstrap_replicates(first, replicates, seed~)
  let second_samples = bootstrap_replicates(second, replicates, seed=seed + 1)
  let estimates = []
  for index = 0; index < replicates; index = index + 1 {
    estimates.push(mean(first_samples[index]) - mean(second_samples[index]))
  }
  let estimate = mean(first) - mean(second)
  let lower = quantile(estimates, (1.0 - confidence) / 2.0)
  let upper = quantile(estimates, (1.0 + confidence) / 2.0)
  { estimate, lower, upper, confidence, replicates }
}

///|
pub fn jackknife_estimates(data : Array[Double]) -> Array[Double] {
  let result = []
  if data.length() <= 1 {
    return result
  }
  for omitted = 0; omitted < data.length(); omitted = omitted + 1 {
    let sample = []
    for index = 0; index < data.length(); index = index + 1 {
      if index != omitted {
        sample.push(data[index])
      }
    }
    result.push(mean(sample))
  }
  result
}

///|
pub fn jackknife_bias_corrected_mean(data : Array[Double]) -> Double {
  if data.length() <= 1 {
    return mean(data)
  }
  let estimates = jackknife_estimates(data)
  let average = mean(estimates)
  let data_mean = mean(data)
  data_mean + data_mean - average
}

///|
pub fn monte_carlo_mean(
  data : Array[Double],
  draws : Int,
  seed : Int,
) -> Double {
  if data.length() == 0 || draws <= 0 {
    return 0.0
  }
  let rng = DeterministicRng::new(seed)
  let mut total = 0.0
  for draw = 0; draw < draws; draw = draw + 1 {
    total += data[rng.next_int(data.length())]
  }
  total / draws.to_double()
}

///|
pub fn bootstrap_standard_error(
  data : Array[Double],
  replicates : Int,
  seed? : Int = 12345,
) -> Double {
  if data.length() == 0 || replicates <= 1 {
    return 0.0
  }
  let samples = bootstrap_replicates(data, replicates, seed~)
  let estimates = []
  for sample in samples {
    estimates.push(mean(sample))
  }
  sample_stddev(estimates)
}