///|
pub(all) struct PhaseSpaceGrid {
  x_bins : Int
  v_bins : Int
  x_min : Double
  x_max : Double
  v_min : Double
  v_max : Double
  counts : Array[Int]
} derive(Debug, ToJson)

///|
pub fn PhaseSpaceGrid::new(
  x_bins : Int,
  v_bins : Int,
  x_min : Double,
  x_max : Double,
  v_min : Double,
  v_max : Double,
) -> PhaseSpaceGrid {
  let nx = x_bins.max(0)
  let nv = v_bins.max(0)
  {
    x_bins: nx,
    v_bins: nv,
    x_min,
    x_max,
    v_min,
    v_max,
    counts: Array::make(nx * nv, 0),
  }
}

///|
pub fn PhaseSpaceGrid::x_width(grid : PhaseSpaceGrid) -> Double {
  safe_ratio(grid.x_max - grid.x_min, grid.x_bins.to_double(), 0.0)
}

///|
pub fn PhaseSpaceGrid::v_width(grid : PhaseSpaceGrid) -> Double {
  safe_ratio(grid.v_max - grid.v_min, grid.v_bins.to_double(), 0.0)
}

///|
pub fn PhaseSpaceGrid::x_index(grid : PhaseSpaceGrid, x : Double) -> Int {
  if grid.x_bins == 0 || grid.x_max <= grid.x_min {
    0
  } else {
    clamp_int(
      ((x - grid.x_min) / (grid.x_max - grid.x_min) * grid.x_bins.to_double()).to_int(),
      0,
      grid.x_bins - 1,
    )
  }
}

///|
pub fn PhaseSpaceGrid::v_index(grid : PhaseSpaceGrid, v : Double) -> Int {
  if grid.v_bins == 0 || grid.v_max <= grid.v_min {
    0
  } else {
    clamp_int(
      ((v - grid.v_min) / (grid.v_max - grid.v_min) * grid.v_bins.to_double()).to_int(),
      0,
      grid.v_bins - 1,
    )
  }
}

///|
pub fn PhaseSpaceGrid::index(
  grid : PhaseSpaceGrid,
  x_index : Int,
  v_index : Int,
) -> Int {
  clamp_int(x_index, 0, grid.x_bins - 1) * grid.v_bins +
  clamp_int(v_index, 0, grid.v_bins - 1)
}

///|
pub fn PhaseSpaceGrid::total_count(grid : PhaseSpaceGrid) -> Int {
  grid.counts.fold(init=0, fn(acc, value) { acc + value })
}

///|
pub fn histogram_particles(
  grid : PhaseSpaceGrid,
  particles : ArrayView[Particle],
) -> PhaseSpaceGrid {
  for particle in particles {
    let index = grid.index(grid.x_index(particle.x), grid.v_index(particle.v))
    if index >= 0 && index < grid.counts.length() {
      grid.counts[index] = grid.counts[index] + 1
    }
  }
  grid
}

///|
pub fn weighted_histogram_particles(
  grid : PhaseSpaceGrid,
  particles : ArrayView[Particle],
) -> Array[Double] {
  let output = zeros(grid.counts.length())
  for particle in particles {
    let index = grid.index(grid.x_index(particle.x), grid.v_index(particle.v))
    if index >= 0 && index < output.length() {
      output[index] = output[index] + particle.weight
    }
  }
  output
}

///|
pub fn phase_space_density(
  grid : PhaseSpaceGrid,
  particles : ArrayView[Particle],
) -> Array[Double] {
  let weighted = weighted_histogram_particles(grid, particles)
  let cell_volume = grid.x_width() * grid.v_width()
  if cell_volume == 0.0 {
    weighted
  } else {
    weighted.map(fn(value) { value / cell_volume })
  }
}

///|
pub fn phase_space_marginal_x(
  grid : PhaseSpaceGrid,
  histogram : ArrayView[Double],
) -> Array[Double] {
  let output = zeros(grid.x_bins)
  for ix in 0.. Array[Double] {
  let output = zeros(grid.v_bins)
  for iv in 0.. Double {
  let marginal = phase_space_marginal_v(grid, histogram)
  let mut numerator = 0.0
  let mut denominator = 0.0
  for iv in 0.. Double {
  let marginal = phase_space_marginal_v(grid, histogram)
  let mut energy = 0.0
  for iv in 0.. String {
  let output = StringBuilder()
  output.write_string("x_bin,v_bin,x,v,value\n")
  for ix in 0.. (Int, Int) {
  let mut best = 0.0
  let mut best_x = 0
  let mut best_v = 0
  for ix in 0.. best {
        best = value
        best_x = ix
        best_v = iv
      }
    }
  }
  (best_x, best_v)
}

///|
pub fn phase_space_entropy(histogram : ArrayView[Double]) -> Double {
  let total = sum_values(histogram)
  if total <= 0.0 {
    0.0
  } else {
    histogram.fold(init=0.0, fn(acc, value) {
      if value <= 0.0 {
        acc
      } else {
        let p = value / total
        acc - p * approximate_log(p)
      }
    })
  }
}