///|
/// Frozen training reference. Neither transform nor public access mutates it.
pub struct ScatterModel {
  reference : Spectrum
  reference_mean : Double
  reference_energy : Double
}

///|
pub(all) struct ScatterResult {
  spectrum : Spectrum
  intercept : Double
  slope : Double
  residual_rms : Double
}

///|
fn require_same_axis(
  left : Spectrum,
  right : Spectrum,
  tolerance : Double,
) -> Unit raise SpectralError {
  finite(tolerance, "axis tolerance")
  if tolerance < 0.0 || tolerance > 1.0e-3 {
    raise Invalid("axis tolerance outside [0,1e-3]")
  }
  if left.length() != right.length() {
    raise Invalid("spectral lengths differ")
  }
  for i in 0.. tolerance * scale {
      raise Invalid("axes differ; align spectra before scatter correction")
    }
  }
}

///|
/// Construct a model from an explicitly supplied reference spectrum.
pub fn scatter_reference(
  reference : Spectrum,
) -> ScatterModel raise SpectralError {
  let center = mean(reference.values)
  let mut energy = 0.0
  for v in reference.values {
    let delta = v - center
    energy = energy + delta * delta
  }
  if energy <= 0.0 ||
    sample_deviation(reference.values) <=
    1.0e-14 * root_mean_square(reference.values) {
    raise Numerical("scatter reference must be nonconstant")
  }
  {
    reference: spectrum(reference.axis, reference.values),
    reference_mean: center,
    reference_energy: checked_result(energy, "reference energy"),
  }
}

///|
/// Fit a mean reference using training data only, on an exactly common axis.
/// At most 256 spectra and 1,048,576 total samples are accepted per fit.
pub fn fit_scatter(
  training : Array[Spectrum],
) -> ScatterModel raise SpectralError {
  if training.length() < 1 || training.length() > 256 {
    raise Budget("training batch must contain 1..256 spectra")
  }
  let n = training[0].length()
  if n * training.length() > 1048576 {
    raise Budget("training batch exceeds sample budget")
  }
  for item in training {
    require_same_axis(training[0], item, 0.0)
  }
  let average = Array::make(n, 0.0)
  let correction = Array::make(n, 0.0)
  for item in training {
    for i in 0.. Spectrum raise SpectralError {
  spectrum(self.reference.axis, self.reference.values)
}

///|
/// Regress sample = intercept + slope*reference, then return
/// (sample-intercept)/slope. Negative or almost-zero physical slopes fail.
pub fn ScatterModel::transform(
  self : ScatterModel,
  input : Spectrum,
  axis_tolerance~ : Double = 0.0,
) -> ScatterResult raise SpectralError {
  require_same_axis(self.reference, input, axis_tolerance)
  let center = mean(input.values)
  let mut covariance = 0.0
  let mut correction = 0.0
  for i in 0.. 1.0e12 {
    raise Numerical("scatter slope must lie in [1e-12,1e12]")
  }
  let intercept = checked_result(
    center - slope * self.reference_mean,
    "scatter intercept",
  )
  let output = Array::make(input.length(), 0.0)
  let residual = Array::make(input.length(), 0.0)
  for i in 0..