moonopt_test.mbt

///|
/// Blackbox tests: they exercise the public surface exactly as a user would, via
/// `@moonopt` and `@model`.
test "production plan is solved to the known optimum" {
  // max 5x + 4y  s.t.  6x + 4y <= 24,  x + 2y <= 6,  x, y >= 0  ->  21 at (3, 1.5)
  let m = @model.Model::new(@model.Sense::Maximize)
  let x = m.add_var("x")
  let y = m.add_var("y")
  m.set_objective([(x, 5.0), (y, 4.0)])
  m.add_named_constraint(
    [(x, 6.0), (y, 4.0)],
    @model.Rel::LessEqual,
    24.0,
    "machine",
  )
  m.add_named_constraint(
    [(x, 1.0), (y, 2.0)],
    @model.Rel::LessEqual,
    6.0,
    "labour",
  )
  let sol = @moonopt.solve(m)
  assert_eq(sol.status, @moonopt.SolveStatus::Optimal)
  assert_true(@core.approx_eq(sol.objective, 21.0))
  assert_eq(sol.values.length(), 2)
  assert_true(@core.approx_eq(sol.values[0], 3.0))
  assert_true(@core.approx_eq(sol.values[1], 1.5))
  assert_eq(sol.message, "")
}

///|
test "minimisation is handled by sign convention" {
  // min 2x + 3y  s.t.  x + y >= 4,  x, y >= 0  ->  8 at (4, 0)
  let m = @model.Model::new(@model.Sense::Minimize)
  let x = m.add_var("x")
  let y = m.add_var("y")
  m.set_objective([(x, 2.0), (y, 3.0)])
  m.add_constraint([(x, 1.0), (y, 1.0)], @model.Rel::GreaterEqual, 4.0)
  let sol = @moonopt.solve(m)
  assert_eq(sol.status, @moonopt.SolveStatus::Optimal)
  assert_true(@core.approx_eq(sol.objective, 8.0))
  assert_true(@core.approx_eq(sol.values[0], 4.0))
  assert_true(@core.approx_eq(sol.values[1], 0.0))
}

///|
test "finite upper bounds become explicit rows" {
  // max x + y  s.t.  x + y <= 10,  x <= 2,  y <= 3  ->  5
  let m = @model.Model::new(@model.Sense::Maximize)
  let x = m.add_var("x", ub=2.0)
  let y = m.add_var("y", ub=3.0)
  m.set_objective([(x, 1.0), (y, 1.0)])
  m.add_constraint([(x, 1.0), (y, 1.0)], @model.Rel::LessEqual, 10.0)
  let sol = @moonopt.solve(m)
  assert_eq(sol.status, @moonopt.SolveStatus::Optimal)
  assert_true(@core.approx_eq(sol.objective, 5.0))
}

///|
test "infeasible model reports Infeasible, not a wrong answer" {
  let m = @model.Model::new(@model.Sense::Maximize)
  let x = m.add_var("x", ub=1.0)
  m.set_objective([(x, 1.0)])
  m.add_constraint([(x, 1.0)], @model.Rel::GreaterEqual, 2.0)
  let sol = @moonopt.solve(m)
  assert_eq(sol.status, @moonopt.SolveStatus::Infeasible)
  assert_eq(sol.values.length(), 0)
  // The wording depends on which layer proved it: the reduction says it proved
  // infeasibility at a named row, the kernel says no feasible point exists. Both
  // have to say something, and both have to be about feasibility.
  assert_true(sol.message.contains("feasible"))
}

///|
test "non-zero lower bounds are supported by the kernel" {
  // min x  s.t.  x >= 0,  x >= 1 (from the bound)  ->  1
  let m = @model.Model::new(@model.Sense::Minimize)
  let x = m.add_var("x", lb=1.0)
  m.set_objective([(x, 1.0)])
  let sol = @moonopt.solve(m)
  assert_eq(sol.status, @moonopt.SolveStatus::Optimal)
  assert_true(@core.approx_eq(sol.objective, 1.0))
  assert_true(@core.approx_eq(sol.values[0], 1.0))
}

///|
test "free variables are supported by the kernel" {
  // max x  s.t.  x <= -2,  x free  ->  -2
  let m = @model.Model::new(@model.Sense::Maximize)
  let x = m.add_var("x", lb=-@model.infinite_bound, ub=@model.infinite_bound)
  m.set_objective([(x, 1.0)])
  m.add_constraint([(x, 1.0)], @model.Rel::LessEqual, -2.0)
  let sol = @moonopt.solve(m)
  assert_eq(sol.status, @moonopt.SolveStatus::Optimal)
  assert_true(@core.approx_eq(sol.objective, -2.0))
  assert_true(@core.approx_eq(sol.values[0], -2.0))
}

///|
/// `max 5x + 4y  s.t.  6x + 4y <= 13`, `x, y` integers in `[0, 2]`.
///
/// The relaxation fills `y` to its bound of 2 and gives `x` what is left — `(5/6, 2)` for
/// 12.1667 — while the integer optimum is 10 at `(2, 0)`. That gap is what makes this a
/// model where the answer has to come from the branch and bound rather than from the
/// relaxation: a relaxation reported as an answer would say 12.1667 and be wrong.
fn integer_mip() -> @model.Model {
  let m = @model.Model::new(@model.Sense::Maximize)
  let x = m.add_var("x", ub=2.0, is_int=true)
  let y = m.add_var("y", ub=2.0, is_int=true)
  m.set_objective([(x, 5.0), (y, 4.0)])
  m.add_constraint([(x, 6.0), (y, 4.0)], @model.Rel::LessEqual, 13.0)
  m
}

///|
/// `max x + y  s.t.  x + y <= 1`, both continuous: the LP the two fields are pinned on.
fn continuous_lp() -> @model.Model {
  let m = @model.Model::new(@model.Sense::Maximize)
  let x = m.add_var("x", ub=1.0)
  let y = m.add_var("y", ub=1.0)
  m.set_objective([(x, 1.0), (y, 1.0)])
  m.add_constraint([(x, 1.0), (y, 1.0)], @model.Rel::LessEqual, 1.0)
  m
}

///|
/// `relaxation` answers the other question a caller can ask about an integer model: not
/// "what is the best whole point" but "what is the bound every whole point is measured
/// against", which is what makes a plan someone already has checkable. The relaxation of
/// `integer_mip()` is the same model without integrality, so its optimum (73/6) bounds the
/// integer answer (10) — from above, because that model maximises; the direction is the
/// model's sense, and getting it wrong is how a "gap" comes out negative. The variables keep
/// their order, so the caller's own point still lives in the same space.
test "the relaxation bounds the whole points and keeps the caller's space" {
  let model = integer_mip()
  let relaxed = @moonopt.relaxation(model)
  assert_eq(relaxed.num_vars(), model.num_vars())
  assert_eq(relaxed.num_constraints(), model.num_constraints())
  for j = 0; j < relaxed.num_vars(); j = j + 1 {
    assert_false(relaxed.get_var(j).is_int)
    assert_true(relaxed.get_var(j).name == model.get_var(j).name)
    assert_true(relaxed.get_var(j).ub == model.get_var(j).ub)
  }
  let bound = @moonopt.solve(relaxed)
  assert_eq(bound.status, @moonopt.SolveStatus::Optimal)
  assert_true(
    @core.approx_eq(bound.objective, 73.0 / 6.0, rel=1.0e-9, abs=1.0e-9),
  )
  // A linear solve is its own bound, which is what makes the two numbers a gap.
  assert_true(
    @core.approx_eq(bound.bound, bound.objective, rel=1.0e-12, abs=0.0),
  )
  // A maximisation's relaxation is an upper bound, and here it is strictly above the integer
  // answer: that difference is the number the relaxation exists to produce.
  let answer = @moonopt.solve(model)
  assert_eq(answer.status, @moonopt.SolveStatus::Optimal)
  assert_true(bound.objective > answer.objective + 1.0e-6)
  // The same model minimised instead: there the relaxation is a lower bound. Both directions
  // are asserted because "a relaxation bounds the whole points" is all the helper promises.
  let minimised = @model.Model::new(@model.Sense::Minimize)
  let a = minimised.add_var("a", ub=2.0, is_int=true)
  let b = minimised.add_var("b", ub=2.0, is_int=true)
  minimised.set_objective([(a, 5.0), (b, 4.0)])
  minimised.add_constraint([(a, 6.0), (b, 4.0)], @model.Rel::GreaterEqual, 13.0)
  let low = @moonopt.solve(@moonopt.relaxation(minimised))
  let whole = @moonopt.solve(minimised)
  assert_eq(low.status, @moonopt.SolveStatus::Optimal)
  assert_eq(whole.status, @moonopt.SolveStatus::Optimal)
  assert_true(low.objective < whole.objective - 1.0e-6)
  // The point of the plan the caller wrote down reads back on the relaxation: the same
  // variable indices, the same rows, so the comparison is meaningful in both directions.
  let plan : Array[Double] = [2.0, 0.0]
  assert_true(
    @core.approx_eq(
      @presolve.objective_value(relaxed, plan),
      10.0,
      rel=1.0e-9,
      abs=1.0e-9,
    ),
  )
  assert_true(@presolve.max_row_violation(relaxed, plan) <= 1.0e-9)
  assert_true(@presolve.max_bound_violation(relaxed, plan) <= 1.0e-9)
}

///|
/// A model with integer variables is answered by the branch and bound, not by its
/// relaxation: the relaxation's optimum here is 73/6 and the answer is 10.
test "an integer model is solved to the integer optimum, not to its relaxation" {
  let sol = @moonopt.solve(integer_mip())
  assert_eq(sol.status, @moonopt.SolveStatus::Optimal)
  assert_true(@core.approx_eq(sol.objective, 10.0, rel=1.0e-9, abs=1.0e-9))
  assert_eq(sol.values.length(), 2)
  // The point is whole, which is the one thing a relaxation answer would not be.
  assert_true(@core.approx_eq(sol.values[0], 2.0, rel=1.0e-9, abs=1.0e-9))
  assert_true(@core.approx_eq(sol.values[1], 0.0, rel=1.0e-9, abs=1.0e-9))
  // Nothing is left open, and the two fields say so.
  assert_true(@core.approx_eq(sol.bound, 10.0, rel=1.0e-9, abs=1.0e-9))
  assert_true(@core.approx_eq(sol.gap, 0.0, rel=0.0, abs=1.0e-12))
  assert_true(sol.nodes > 0)
  // A linear solve is its own bound, so the same two fields mean the same thing there.
  let linear = @moonopt.solve(continuous_lp())
  assert_eq(linear.status, @moonopt.SolveStatus::Optimal)
  assert_true(
    @core.approx_eq(linear.bound, linear.objective, rel=1.0e-12, abs=0.0),
  )
  assert_true(@core.approx_eq(linear.gap, 0.0, rel=0.0, abs=1.0e-12))
  assert_eq(linear.nodes, 0)
}

///|
/// A node budget is a budget, not an answer: the run stops as `NodeLimit` with whatever integer
/// point it has and the bound still open.
///
/// The point may exist even at one relaxation — the rounding heuristic reads one out of the root's
/// own relaxation — which is why what is asserted here is the contract rather than the emptiness:
/// the status is a limit, the bound is still the relaxation's, and a point, when there is one, is a
/// point of the model and not an optimality claim.
test "a node budget is reported as a budget rather than as an optimality claim" {
  // Cuts off: this is about what one relaxation's budget leaves behind, and a round of
  // cuts at the root would spend a second relaxation to tighten the same node.
  let model = integer_mip()
  let sol = @moonopt.solve_with(model, {
    max_iterations: 10000,
    presolve: true,
    max_nodes: 1,
    cut_rounds: 0,
  })
  assert_eq(sol.status, @moonopt.SolveStatus::NodeLimit)
  assert_true(@core.approx_eq(sol.bound, 73.0 / 6.0, rel=1.0e-9, abs=1.0e-9))
  assert_true(sol.message.length() > 0)
  assert_eq(sol.nodes, 1)
  if sol.values.length() > 0 {
    // Whatever it has is a real integer point of the model, checked without the search's own
    // routines, and its objective is that point's — not the relaxation's.
    assert_eq(sol.values.length(), model.num_vars())
    for j = 0; j < model.num_vars(); j = j + 1 {
      if model.get_var(j).is_int {
        // Whole to the last place, measured the way the search measures it: the distance to the
        // nearest whole number, which is what `fractional_count` in the kernel decides a node on.
        let truncated = Double::from_int(sol.values[j].to_int())
        let distance = @core.fmin(
          @core.fabs(sol.values[j] - truncated),
          @core.fabs(sol.values[j] - truncated - 1.0),
        )
        assert_true(distance <= 1.0e-9)
      }
    }
    let mut acc = 0.0
    for term in model.objective() {
      let (index, coefficient) = term
      acc = acc + coefficient * sol.values[index]
    }
    assert_true(@core.approx_eq(sol.objective, acc, rel=1.0e-9, abs=1.0e-9))
    // And the gap is the distance between those two numbers, not a claim that they meet.
    assert_true(sol.gap > 0.0)
  } else {
    assert_true(sol.gap > 1.0e299)
  }
}

///|
/// A model whose relaxation is feasible and whose integer points are not: the tree
/// exhausts itself, so `Infeasible` is a proof here in a way it never is for a single
/// relaxation.
test "an integer model with no integer point is proved infeasible" {
  let m = @model.Model::new(@model.Sense::Minimize)
  let x = m.add_var("x", ub=4.0, is_int=true)
  let y = m.add_var("y", ub=4.0, is_int=true)
  m.set_objective([(x, 1.0), (y, 1.0)])
  m.add_constraint([(x, 2.0), (y, 2.0)], @model.Rel::Equal, 1.0)
  let sol = @moonopt.solve(m)
  assert_eq(sol.status, @moonopt.SolveStatus::Infeasible)
  assert_eq(sol.values.length(), 0)
  assert_true(sol.nodes > 0)
}

///|
/// An unbounded relaxation whose point and ray are integral on the integer variable IS an
/// unbounded model, and the public status says so: `max x + y  s.t.  x - y <= 3` with `x`
/// integer and no upper bound on either variable has the integral ray `(1, 1)` from the
/// integral point `(0, 0)`. The public `Solution` carries no values — a point and a ray are
/// the proof, and they live on the low-level `MipResult` the way the kernel's do — so the
/// status and the message are what the caller reads.
test "an integer model an integral ray proves unbounded is reported unbounded" {
  let m = @model.Model::new(@model.Sense::Maximize)
  let x = m.add_var("x", is_int=true)
  let y = m.add_var("y")
  m.set_objective([(x, 1.0), (y, 1.0)])
  m.add_constraint([(x, 1.0), (y, -1.0)], @model.Rel::LessEqual, 3.0)
  let sol = @moonopt.solve(m)
  assert_eq(sol.status, @moonopt.SolveStatus::Unbounded)
  assert_eq(sol.values.length(), 0)
  assert_true(sol.message.length() > 0)
}

///|
test "an ill-formed model is rejected before any solving happens" {
  let m = @model.Model::new(@model.Sense::Minimize)
  let _ = m.add_var("x", lb=5.0, ub=1.0)
  let sol = @moonopt.solve(m)
  assert_eq(sol.status, @moonopt.SolveStatus::NotSolved)
  assert_true(sol.message.length() > 0)
}

///|
test "iteration limit is surfaced instead of being ignored" {
  let m = @model.Model::new(@model.Sense::Maximize)
  let x = m.add_var("x")
  m.set_objective([(x, 1.0)])
  m.add_constraint([(x, 1.0)], @model.Rel::LessEqual, 1.0)
  // Presolve is on by default and would settle this model by itself (the row is a
  // bound, and the variable then appears in no row), which is a different thing
  // from what this test is about: the kernel's own iteration limit.
  let sol = @moonopt.solve_with(m, {
    max_iterations: 0,
    presolve: false,
    max_nodes: 20000,
    cut_rounds: 2,
  })
  assert_eq(sol.status, @moonopt.SolveStatus::NotSolved)
  assert_true(sol.message.length() > 0)
}

///|
/// Presolve is the default path, and it is allowed to finish a model on its own:
/// `max x` with `x <= 1` is decided by the reduction, so the answer comes back with
/// no kernel run at all — and with the reconstructed values indexed by the
/// original variables.
test "presolve solves a model on its own and reconstructs the value" {
  let m = @model.Model::new(@model.Sense::Maximize)
  let x = m.add_var("x")
  m.set_objective([(x, 1.0)])
  m.add_constraint([(x, 1.0)], @model.Rel::LessEqual, 1.0)
  let sol = @moonopt.solve(m)
  assert_eq(sol.status, @moonopt.SolveStatus::Optimal)
  assert_eq(sol.iterations, 0)
  assert_true(@core.approx_eq(sol.objective, 1.0, rel=1.0e-9, abs=1.0e-9))
  assert_eq(sol.values.length(), 1)
  assert_true(@core.approx_eq(sol.values[0], 1.0, rel=1.0e-9, abs=1.0e-9))
}

///|
/// The two paths have to agree. With a fixed variable the reduced model no longer
/// contains it, so the objective has to pick up its contribution from the offset
/// and the value vector has to put it back in the right place.
test "the presolved path and the kernel-only path agree" {
  // min x + y  s.t.  x + y >= 5,  x = 2,  y >= 0   ->   x = 2, y = 3, objective 5
  let with_presolve = @model.Model::new(@model.Sense::Minimize)
  let x1 = with_presolve.add_var("x", lb=2.0, ub=2.0)
  let y1 = with_presolve.add_var("y", lb=0.0)
  with_presolve.set_objective([(x1, 1.0), (y1, 1.0)])
  with_presolve.add_constraint(
    [(x1, 1.0), (y1, 1.0)],
    @model.Rel::GreaterEqual,
    5.0,
  )
  let presolved = @moonopt.solve(with_presolve)
  assert_eq(presolved.status, @moonopt.SolveStatus::Optimal)
  assert_true(@core.approx_eq(presolved.objective, 5.0, rel=1.0e-9, abs=1.0e-9))
  assert_eq(presolved.values.length(), 2)
  assert_true(@core.approx_eq(presolved.values[0], 2.0, rel=1.0e-9, abs=1.0e-9))
  assert_true(@core.approx_eq(presolved.values[1], 3.0, rel=1.0e-9, abs=1.0e-9))

  let direct = @model.Model::new(@model.Sense::Minimize)
  let x2 = direct.add_var("x", lb=2.0, ub=2.0)
  let y2 = direct.add_var("y", lb=0.0)
  direct.set_objective([(x2, 1.0), (y2, 1.0)])
  direct.add_constraint([(x2, 1.0), (y2, 1.0)], @model.Rel::GreaterEqual, 5.0)
  let kernel_only = @moonopt.solve_with(direct, {
    max_iterations: 10000,
    presolve: false,
    max_nodes: 20000,
    cut_rounds: 2,
  })
  assert_eq(kernel_only.status, @moonopt.SolveStatus::Optimal)
  assert_true(
    @core.approx_eq(kernel_only.objective, 5.0, rel=1.0e-9, abs=1.0e-9),
  )
}

///|
/// A verdict presolve can prove on its own is reported as such, without a kernel
/// run: this variable appears in no row, its cost improves with the value growing,
/// and nothing bounds it above.
test "a model presolve proves unbounded is reported without a kernel run" {
  let m = @model.Model::new(@model.Sense::Minimize)
  let x = m.add_var("x", lb=0.0)
  let y = m.add_var("y", lb=0.0)
  m.set_objective([(x, -1.0), (y, 1.0)])
  m.add_constraint([(y, 1.0)], @model.Rel::GreaterEqual, 4.0)
  let sol = @moonopt.solve(m)
  assert_eq(sol.status, @moonopt.SolveStatus::Unbounded)
  assert_true(sol.message.length() > 0)
}