///|
pub(all) struct RandomStream {
  mut state : Int
} derive(Debug, ToJson)

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

///|
pub fn RandomStream::next_int(stream : RandomStream) -> Int {
  stream.state = (stream.state * 997 + 37) % 1000003
  stream.state
}

///|
pub fn RandomStream::next_unit(stream : RandomStream) -> Double {
  stream.next_int().to_double() / 1000003.0
}

///|
pub fn RandomStream::next_signed(stream : RandomStream) -> Double {
  2.0 * stream.next_unit() - 1.0
}

///|
pub fn sample_uniform(
  count : Int,
  seed : Int,
  low : Double,
  high : Double,
) -> Array[Double] {
  let stream = RandomStream::new(seed)
  let n = if count < 0 { 0 } else { count }
  Array::makei(n, fn(_) { lerp(low, high, stream.next_unit()) })
}

///|
pub fn sample_uniform_grid(
  grid : Grid1D,
  count : Int,
  seed : Int,
) -> Array[Double] {
  sample_uniform(count, seed, 0.0, grid.length)
}

///|
pub fn sample_normal(
  count : Int,
  seed : Int,
  mean_value : Double,
  standard_deviation : Double,
) -> Array[Double] {
  let stream = RandomStream::new(seed)
  let n = if count < 0 { 0 } else { count }
  let output = zeros(n)
  for i in 0.. Array[Double] {
  sample_normal(count, seed, drift, spread)
}

///|
pub fn sample_two_beam(
  count : Int,
  seed : Int,
  drift : Double,
  spread : Double,
) -> Array[Double] {
  let stream = RandomStream::new(seed)
  let n = if count < 0 { 0 } else { count }
  let output = zeros(n)
  for i in 0.. Array[Double] {
  let n = if count < 0 { 0 } else { count }
  Array::makei(n, fn(i) {
    let fraction = (i.to_double() + 0.5) / n.max(1).to_double()
    (fraction * length + phase).wrap_period(length)
  })
}

///|
fn Double::wrap_period(value : Double, period : Double) -> Double {
  if period <= 0.0 {
    value
  } else {
    let turns = @math.floor(value / period)
    let result = value - turns * period
    if result < 0.0 {
      result + period
    } else if result >= period {
      result - period
    } else {
      result
    }
  }
}

///|
pub fn sample_profile(
  grid : Grid1D,
  profile : Profile1D,
  count : Int,
  seed : Int,
) -> Array[Double] {
  let positions = sample_uniform_grid(grid, count, seed)
  positions.map(fn(position) { profile_value(profile, position) })
}

///|
pub fn sample_rejection(
  grid : Grid1D,
  profile : Profile1D,
  count : Int,
  seed : Int,
  upper : Double,
) -> Array[Double] {
  let stream = RandomStream::new(seed)
  let output : Array[Double] = []
  let limit = if count < 0 { 0 } else { count }
  let mut attempts = 0
  while output.length() < limit && attempts < limit * 100 + 1 {
    let x = stream.next_unit() * grid.length
    let y = stream.next_unit() * upper
    if y <= profile_value(profile, x) {
      output.push(x)
    }
    attempts = attempts + 1
  }
  output
}

///|
pub fn sample_quantiles(
  values : ArrayView[Double],
  count : Int,
) -> Array[Double] {
  if values.length() == 0 || count <= 0 {
    []
  } else {
    let sorted = values.to_owned()
    sorted.sort()
    Array::makei(count, fn(i) {
      let fraction = if count == 1 {
        0.0
      } else {
        i.to_double() / (count - 1).to_double()
      }
      let index = (fraction * (sorted.length() - 1).to_double()).to_int()
      sorted[clamp_int(index, 0, sorted.length() - 1)]
    })
  }
}

///|
pub fn sample_histogram(
  values : ArrayView[Double],
  low : Double,
  high : Double,
  bins : Int,
) -> Array[Int] {
  let output = Array::make(bins.max(0), 0)
  if high > low && bins > 0 {
    for value in values {
      let position = ((value - low) / (high - low) * bins.to_double()).to_int()
      let index = clamp_int(position, 0, bins - 1)
      output[index] = output[index] + 1
    }
  }
  output
}

///|
pub fn histogram_edges(
  low : Double,
  high : Double,
  bins : Int,
) -> Array[Double] {
  if bins <= 0 {
    []
  } else {
    Array::makei(bins + 1, fn(i) {
      lerp(low, high, i.to_double() / bins.to_double())
    })
  }
}

///|
pub fn sample_mean(values : ArrayView[Double]) -> Double {
  mean(values)
}

///|
pub fn sample_variance(values : ArrayView[Double]) -> Double {
  velocity_variance(values.map(fn(value) { Particle::new(x=0.0, v=value) }))
}

///|
pub fn sample_span(values : ArrayView[Double]) -> Double {
  max_value(values, 0.0) - min_value(values, 0.0)
}

///|
pub fn deterministic_seed(name : String, salt : Int) -> Int {
  let mut hash = 2166136
  for i in 0.. Array[Particle] {
  let n = if count < 0 { 0 } else { count }
  Array::makei(n, fn(i) {
    Particle::new(x=grid.position(i % grid.cells), v=velocity)
  })
}

///|
pub fn sample_particles(
  grid : Grid1D,
  count : Int,
  seed : Int,
  drift : Double,
  spread : Double,
) -> Array[Particle] {
  let positions = sample_uniform_grid(grid, count, seed)
  let velocities = sample_normal(count, seed + 17, drift, spread)
  Array::makei(positions.length(), fn(i) {
    Particle::new(x=positions[i], v=velocities[i])
  })
}

///|
pub fn sample_weighted_particles(
  grid : Grid1D,
  count : Int,
  seed : Int,
  profile : Profile1D,
) -> Array[Particle] {
  let positions = sample_uniform_grid(grid, count, seed)
  Array::makei(positions.length(), fn(i) {
    let weight = profile_value(profile, positions[i]).max(0.0)
    Particle::new(x=positions[i], v=0.0, weight~)
  })
}

///|
pub fn sample_clone(particles : ArrayView[Particle]) -> Array[Particle] {
  particles.map(fn(p) {
    Particle::new(x=p.x, v=p.v, weight=p.weight, charge=p.charge, mass=p.mass)
  })
}