///|
pub struct Spillover {
  channels : Array[String]
  coefficients : Array[Double]
} derive(Debug, ToJson)

///|
pub fn Dataset::channel_info(self : Dataset) -> Array[Channel] {
  self.channels.copy()
}

///|
pub fn Dataset::keywords(self : Dataset) -> Array[(String, String)] {
  self.metadata.copy()
}

///|
pub fn Dataset::warnings(self : Dataset) -> Array[String] {
  self.diagnostics.copy()
}

///|
pub fn Dataset::analysis_keywords(self : Dataset) -> Array[(String, String)] {
  self.analysis.copy()
}

///|
pub fn Dataset::keyword(self : Dataset, key : String) -> String? {
  lookup(self.metadata, key)
}

///|
pub fn Dataset::spillover(self : Dataset) -> Spillover? raise FcsError {
  let modern = self.keyword("$SPILLOVER")
  let legacy = self.keyword("$SPILL")
  if modern is Some(a) && legacy is Some(b) && a != b {
    raise Invalid("conflicting SPILL and SPILLOVER")
  }
  let value = match modern {
    Some(v) => v
    None =>
      match legacy {
        Some(v) => v
        None => return None
      }
  }
  let words = value.split(",").map(s => s.trim().to_owned()).to_array()
  let n = natural(words[0])
  if n < 1 || n > 128 {
    raise Limit("spillover dimension must be 1..128")
  }
  if words.length() != 1 + n + n * n {
    raise Invalid("spillover matrix dimension mismatch")
  }
  let channels = []
  for i = 0; i < n; i = i + 1 {
    let name = words[i + 1]
    if channels.contains(name) {
      raise Invalid("duplicate spillover channel")
    }
    let mut found = 0
    for c in self.channels {
      if c.name == name {
        found = found + 1
      }
    }
    if found != 1 {
      raise Invalid(
        "spillover channel must identify exactly one parameter: \{name}",
      )
    }
    channels.push(name)
  }
  let coefficients = Array::makei(n * n, i => number(words[1 + n + i]))
  Some({ channels, coefficients, })
}

///|
fn Dataset::check_channel(self : Dataset, channel : Int) -> Unit raise FcsError {
  if channel < 0 || channel >= self.channels.length() {
    raise Invalid("channel index out of bounds")
  }
}

///|
/// Raw instrument value; integer padding bits above ceil(log2(PnR)) are masked.
pub fn Dataset::value(
  self : Dataset,
  event : Int,
  channel : Int,
  scaled? : Bool = false,
) -> Double raise FcsError {
  self.check_channel(channel)
  if event < 0 || event >= self.events {
    raise Invalid("event index out of bounds")
  }
  let c = self.channels[channel]
  let bits = uword(
    self.raw,
    event * self.width + self.offsets[channel],
    c.bits / 8,
    self.order,
  )
  let x = match self.kind {
    Float32 => Float::reinterpret_from_uint(bits.to_uint()).to_double()
    Float64 => bits.reinterpret_as_double()
    Integer => {
      let mut ceiling = 1UL
      while ceiling.to_double() < c.range {
        ceiling = ceiling << 1
      }
      (bits & (ceiling - 1UL)).to_double()
    }
  }
  if !scaled {
    return x
  }
  // Time is an explicitly named channel; amplifier gain does not apply to time.
  if c.name.to_upper() == "TIME" {
    let step = number(self.keyword("$TIMESTEP").unwrap_or("1"))
    if step <= 0.0 {
      raise Invalid("TIMESTEP must be positive")
    }
    return x * step
  }
  if c.gain <= 0.0 {
    raise Invalid("non-positive gain prevents scaled access")
  }
  (if c.decades > 0.0 {
    c.zero * @math.pow(10.0, c.decades * x / c.range)
  } else {
    x
  }) /
  c.gain
}

///|
/// Welford sample variance (n-1); non-finite values are counted separately.
pub fn Dataset::statistics(
  self : Dataset,
  channel : Int,
  scaled? : Bool = false,
) -> ChannelStats raise FcsError {
  self.check_channel(channel)
  let mut count = 0
  let mut nonfinite = 0
  let mut lo = 0.0
  let mut hi = 0.0
  let mut mean = 0.0
  let mut m2 = 0.0
  for e = 0; e < self.events; e = e + 1 {
    let x = self.value(e, channel, scaled~)
    if !finite(x) {
      nonfinite = nonfinite + 1
      continue
    }
    count = count + 1
    if count == 1 {
      lo = x
      hi = x
    }
    if x < lo {
      lo = x
    }
    if x > hi {
      hi = x
    }
    let delta = x - mean
    mean = mean + delta / count.to_double()
    m2 = m2 + delta * (x - mean)
  }
  {
    count,
    nonfinite,
    minimum: if count > 0 {
      Some(lo)
    } else {
      None
    },
    maximum: if count > 0 {
      Some(hi)
    } else {
      None
    },
    mean: if count > 0 && finite(mean) {
      Some(mean)
    } else {
      None
    },
    variance: if count > 1 && finite(m2) {
      Some(m2 / (count - 1).to_double())
    } else {
      None
    },
  }
}

///|
/// Intersection of half-open intervals [lower, upper). NaN and infinities excluded.
pub fn Dataset::rectangle(
  self : Dataset,
  bounds : Array[Bounds],
  scaled? : Bool = false,
) -> Array[Int] raise FcsError {
  if bounds.is_empty() {
    raise Invalid("at least one gate interval is required")
  }
  for b in bounds {
    self.check_channel(b.channel)
    if !finite(b.lower) || !finite(b.upper) || b.lower >= b.upper {
      raise Invalid("invalid gate bounds")
    }
  }
  let selected = []
  for e = 0; e < self.events; e = e + 1 {
    let mut keep = true
    for b in bounds {
      let v = self.value(e, b.channel, scaled~)
      if !finite(v) || v < b.lower || v >= b.upper {
        keep = false
        break
      }
    }
    if keep {
      selected.push(e)
    }
  }
  selected
}

///|
/// Even-odd polygon gate; exact floating-point edge points included, no epsilon.
pub fn Dataset::polygon(
  self : Dataset,
  x_channel : Int,
  y_channel : Int,
  vertices : Array[(Double, Double)],
  scaled? : Bool = false,
) -> Array[Int] raise FcsError {
  self.check_channel(x_channel)
  self.check_channel(y_channel)
  if x_channel == y_channel || vertices.length() < 3 || vertices.length() > 4096 {
    raise Invalid("invalid polygon dimensions")
  }
  for (x, y) in vertices {
    if !finite(x) || !finite(y) || x.abs() > 1.0e100 || y.abs() > 1.0e100 {
      raise Invalid("invalid polygon vertex")
    }
  }
  let selected = []
  for e = 0; e < self.events; e = e + 1 {
    let x = self.value(e, x_channel, scaled~)
    let y = self.value(e, y_channel, scaled~)
    if !finite(x) || !finite(y) || x.abs() > 1.0e100 || y.abs() > 1.0e100 {
      continue
    }
    let mut inside = false
    for i = 0; i < vertices.length(); i = i + 1 {
      let (ax, ay) = vertices[i]
      let (bx, by) = vertices[(i + 1) % vertices.length()]
      if (x - ax) * (by - ay) == (y - ay) * (bx - ax) &&
        x >= ax.min(bx) &&
        x <= ax.max(bx) &&
        y >= ay.min(by) &&
        y <= ay.max(by) {
        inside = true
        break
      }
      if (ay > y) != (by > y) && x < (bx - ax) * (y - ay) / (by - ay) + ax {
        inside = !inside
      }
    }
    if inside {
      selected.push(e)
    }
  }
  selected
}