///|
fn valid_bounds(bounds : Bounds, nonnegative~ : Bool) -> Bool {
  bounds.lower <= bounds.upper &&
  bounds.lower >= (if nonnegative { 0 } else { -1000000000 }) &&
  bounds.upper <= 1000000000
}

///|
fn scale(form : Form, factor : R) -> Form {
  {
    constant: form.constant * factor,
    coefficients: form.coefficients.map(fn(x) { x * factor }),
  }
}

///|
fn plus(a : Form, b : Form) -> Form {
  {
    constant: a.constant + b.constant,
    coefficients: Array::makei(a.coefficients.length(), fn(i) {
      a.coefficients[i] + b.coefficients[i]
    }),
  }
}

///|
fn input_form(
  id : String,
  bounds : Bounds,
  parameters : Array[Parameter],
) -> Form {
  let coefficients = Array::make(parameters.length(), r(0))
  let mut constant = r(bounds.lower)
  for i = 0; i < parameters.length(); i = i + 1 {
    if parameters[i].id == id {
      coefficients[i] = r(1)
      constant = r(0)
    }
  }
  { constant, coefficients, }
}

///|
fn inequality(id : String, form : Form, upper : R) -> Row {
  { id, coefficients: form.coefficients, upper: upper - form.constant, }
}

///|
fn node_index(bottles : Array[Bottle], name : String) -> Int32? {
  for i = 0; i < bottles.length(); i = i + 1 {
    if bottles[i].id == name {
      return Some(i)
    }
  }
  None
}

///|
fn compile_study(study : AliquotStudy) -> Result[Compiled, String] {
  if study.schema != "aqueous-aliquot/v1" {
    return Err("schema must be aqueous-aliquot/v1")
  }
  if study.source.id.length() == 0 ||
    study.source.volume_ul <= 0 ||
    study.source.volume_ul > 1000000000 ||
    study.steps.length() > 16 ||
    study.assays.length() > 32 {
    return Err("source needs an id and 1..1e9 uL; at most 16 steps/32 assays")
  }
  let source_bounds = [
    study.source.calcium_nmol_l,
    study.source.dic_nmol_l,
    study.source.ta_neq_l,
    study.source.sodium_nmol_l,
    study.source.chloride_nmol_l,
  ]
  let parameters : Array[Parameter] = []
  for i = 0; i < 5; i = i + 1 {
    let bounds = source_bounds[i]
    if !valid_bounds(bounds, nonnegative=i != 2) {
      return Err("invalid source bounds for " + analytes[i])
    }
    if bounds.lower != bounds.upper {
      parameters.push({ id: "source:" + analytes[i], bounds, })
    }
  }
  let ids = [study.source.id]
  for step in study.steps {
    if step.id.length() == 0 ||
      ids.contains(step.id) ||
      !ids.contains(step.parent) {
      return Err("step ids must be unique and parents must precede children")
    }
    if !["take", "dilute", "evaporate", "co2", "calcite", "hcl", "naoh"].contains(
        step.kind,
      ) {
      return Err("unsupported process: " + step.kind)
    }
    if !valid_bounds(
        step.amount_nmol,
        nonnegative=step.kind != "co2" && step.kind != "calcite",
      ) ||
      step.volume_ul < 0 ||
      step.volume_ul > 1000000000 {
      return Err("invalid amount or volume at " + step.id)
    }
    if ["take", "dilute", "evaporate"].contains(step.kind) {
      if step.volume_ul == 0 ||
        step.amount_nmol.lower != 0 ||
        step.amount_nmol.upper != 0 {
        return Err("volume-only steps require positive volume and zero amount")
      }
    } else if ["co2", "calcite"].contains(step.kind) && step.volume_ul != 0 {
      return Err("co2/calcite steps use zero volume in this finite model")
    }
    if step.amount_nmol.lower != step.amount_nmol.upper {
      parameters.push({ id: "step:" + step.id, bounds: step.amount_nmol, })
    }
    ids.push(step.id)
  }
  if parameters.length() > 10 {
    return Err("at most 10 uncertain source/process quantities are supported")
  }
  let bottles : Array[Bottle] = [
    {
      id: study.source.id,
      volume_ul: study.source.volume_ul,
      amounts: Array::makei(5, fn(i) {
        scale(
          input_form("source:" + analytes[i], source_bounds[i], parameters),
          r(study.source.volume_ul) / r(1000000),
        )
      }),
    },
  ]
  for step in study.steps {
    let parent = match node_index(bottles, step.parent) {
      Some(i) => bottles[i]
      None => return Err("parent not found")
    }
    let mut volume = parent.volume_ul
    let amounts = parent.amounts.copy()
    if step.kind == "take" {
      if step.volume_ul > parent.volume_ul {
        return Err("aliquot exceeds parent volume")
      }
      volume = step.volume_ul
      for i = 0; i < 5; i = i + 1 {
        amounts[i] = scale(amounts[i], r(volume) / r(parent.volume_ul))
      }
    } else if ["dilute", "hcl", "naoh"].contains(step.kind) {
      if volume > 1000000000 - step.volume_ul {
        return Err("prepared volume exceeds 1e9 uL")
      }
      volume = volume + step.volume_ul
    } else if step.kind == "evaporate" {
      volume = volume - step.volume_ul
      if volume <= 0 {
        return Err("evaporation leaves no liquid")
      }
    }
    let dose = input_form("step:" + step.id, step.amount_nmol, parameters)
    let signature = match step.kind {
      "co2" => [0, 1, 0, 0, 0]
      "calcite" => [1, 1, 2, 0, 0]
      "hcl" => [0, 0, -1, 0, 1]
      "naoh" => [0, 0, 1, 1, 0]
      _ => [0, 0, 0, 0, 0]
    }
    for i = 0; i < 5; i = i + 1 {
      amounts[i] = plus(amounts[i], scale(dose, r(signature[i])))
    }
    bottles.push({ id: step.id, volume_ul: volume, amounts, })
  }
  // Children represent actual liquid allocation, not duplicated full bottles.
  for parent in bottles {
    let mut children = 0
    let mut takes = 0
    let mut allocated = 0
    for step in study.steps {
      if step.parent == parent.id {
        children = children + 1
        if step.kind == "take" {
          takes = takes + 1
          if allocated > parent.volume_ul - step.volume_ul {
            return Err("sibling aliquots exceed their one shared parent volume")
          }
          allocated = allocated + step.volume_ul
        }
      }
    }
    if children > 1 && takes != children {
      return Err("multiple children must be explicit take steps")
    }
  }
  let rows : Array[Row] = []
  for i = 0; i < parameters.length(); i = i + 1 {
    let p = parameters[i]
    let f = input_form(p.id, p.bounds, parameters)
    rows.push(inequality("bound:" + p.id + ":upper", f, r(p.bounds.upper)))
    rows.push(
      inequality(
        "bound:" + p.id + ":lower",
        scale(f, r(-1)),
        -r(p.bounds.lower),
      ),
    )
  }
  for bottle in bottles {
    for i = 0; i < 5; i = i + 1 {
      if i != 2 {
        rows.push(
          inequality(
            "nonnegative:" + bottle.id + ":" + analytes[i],
            scale(bottle.amounts[i], r(-1)),
            r(0),
          ),
        )
      }
    }
  }
  let origin = bottles[0].amounts
  let charge = plus(
    plus(scale(origin[0], r(2)), origin[3]),
    scale(plus(origin[4], origin[2]), r(-1)),
  )
  rows.push(inequality("source:carbonate-charge:upper", charge, r(0)))
  rows.push(
    inequality("source:carbonate-charge:lower", scale(charge, r(-1)), r(0)),
  )
  let assay_ids : Array[String] = []
  for assay in study.assays {
    if assay.id.length() == 0 || assay_ids.contains(assay.id) {
      return Err("assay ids must be nonempty and unique")
    }
    assay_ids.push(assay.id)
    let bottle = match node_index(bottles, assay.node) {
      Some(i) => bottles[i]
      None => return Err("assay refers to an unknown preparation state")
    }
    let mut ai = -1
    for i = 0; i < 5; i = i + 1 {
      if analytes[i] == assay.analyte {
        ai = i
      }
    }
    if ai < 0 {
      return Err("assay must be Ca, DIC, TA, Na or Cl; pH is not additive")
    }
    if !valid_bounds(assay.bounds, nonnegative=ai != 2) {
      return Err("invalid assay bounds")
    }
    let f = scale(bottle.amounts[ai], r(1000000) / r(bottle.volume_ul))
    rows.push(
      inequality("assay:" + assay.id + ":upper", f, r(assay.bounds.upper)),
    )
    rows.push(
      inequality(
        "assay:" + assay.id + ":lower",
        scale(f, r(-1)),
        -r(assay.bounds.lower),
      ),
    )
  }
  Ok({ parameters, bottles, rows, })
}

///|
fn evaluate_form(form : Form, values : Array[R]) -> R {
  let mut value = form.constant
  for i = 0; i < values.length(); i = i + 1 {
    value = value + form.coefficients[i] * values[i]
  }
  value
}

///|
fn history(compiled : Compiled, values : Array[R]) -> HandlingHistory {
  {
    model: "Ca-carbonate-Na-Cl-water; exact declared volumes; no untracked ions or processes",
    parameters: Array::makei(values.length(), fn(i) {
      { parameter: compiled.parameters[i].id, value: exact(values[i]), }
    }),
    samples: compiled.bottles.map(fn(bottle) {
      let c = bottle.amounts.map(fn(f) {
        exact(evaluate_form(f, values) * r(1000000) / r(bottle.volume_ul))
      })
      {
        id: bottle.id,
        volume_ul: bottle.volume_ul,
        calcium_nmol_l: c[0],
        dic_nmol_l: c[1],
        ta_neq_l: c[2],
        sodium_nmol_l: c[3],
        chloride_nmol_l: c[4],
      }
    }),
  }
}

///|
/// Replay one complete shared preparation history, including all assay and
/// nonnegative-inventory checks. Missing/extra choices are rejected.
pub fn replay_handling(
  study : AliquotStudy,
  choices : Array[ParameterValue],
) -> Result[HandlingHistory, String] {
  let compiled = match compile_study(study) {
    Ok(c) => c
    Err(e) => return Err(e)
  }
  if choices.length() != compiled.parameters.length() {
    return Err("provide exactly one choice for every uncertain parameter")
  }
  let values = Array::make(compiled.parameters.length(), r(0))
  for i = 0; i < compiled.parameters.length(); i = i + 1 {
    let p = compiled.parameters[i]
    let matching = choices.filter(fn(x) { x.parameter == p.id })
    if matching.length() != 1 {
      return Err("missing or duplicate parameter: " + p.id)
    }
    values[i] = match parse_exact(matching[0].value) {
      Ok(value) => value
      Err(e) => return Err(e)
    }
  }
  for row in compiled.rows {
    let f = Form::{ constant: r(0), coefficients: row.coefficients, }
    if evaluate_form(f, values) > row.upper {
      return Err("history violates " + row.id)
    }
  }
  Ok(history(compiled, values))
}