Changing a model does not invalidate stored dual ray

[odow]

I can reproduce this issue in 1.8.1 but not on 1.8.0.

See this example:

julia> using HiGHS

julia> begin
         ret = highs = Highs_create()
         ret = Highs_setStringOptionValue(highs, "presolve", "off")
         ret = Highs_changeObjectiveOffset(highs, 0.0)
         ret = Highs_addCol(highs, 0.0, 0.0, 0.0, 0, C_NULL, C_NULL)
         ret = Highs_addCol(highs, 0.0, 0.0, 0.0, 0, C_NULL, C_NULL)
         ret = Highs_addCol(highs, -1.0, 0.0, Inf, 0, C_NULL, C_NULL)
         ret = Highs_addCol(highs, -1.0, 0.0, Inf, 0, C_NULL, C_NULL)
         ret = Highs_addRow(highs, 0.0, 0.0, 2, Int32[2, 3], [1.0, -1.0])
         ret = Highs_addRow(highs, 1.0, Inf, 2, Int32[2, 3], [1.0, 1.0])
         ret = Highs_addRow(highs, -Inf, 0.0, 2, Int32[0, 2], [-2.0, 1.0])
         ret = Highs_addRow(highs, -Inf, 0.0, 2, Int32[1, 3], [-3.0, 1.0])
         ret = Highs_run(highs)
         has_dual_ray = Base.RefValue{Int32}(0)
         dual_ray_value = [0.0, 0.0, 0.0, 0.0]
         ret = Highs_getDualRay(highs, has_dual_ray, dual_ray_value)
         @show has_dual_ray, dual_ray_value
         ret = Highs_changeColBounds(highs, 1, 1.0, 1.0)
         ret = Highs_run(highs)
         ret = Highs_getDualRay(highs, has_dual_ray, dual_ray_value)
         @show has_dual_ray, dual_ray_value
       end
Running HiGHS 1.9.0 (git hash: 66f735e60): Copyright (c) 2024 HiGHS under MIT licence terms
Coefficient ranges:
  Matrix [1e+00, 3e+00]
  Cost   [1e+00, 1e+00]
  Bound  [0e+00, 0e+00]
  RHS    [1e+00, 1e+00]
Solving LP without presolve, or with basis, or unconstrained
Using EKK dual simplex solver - serial
  Iteration        Objective     Infeasibilities num(sum)
          0    -1.9999968934e+00 Ph1: 3(3); Du: 2(2) 0s
          2     0.0000000000e+00 0s
Model status        : Infeasible
Simplex   iterations: 2
Objective value     :  0.0000000000e+00
Relative P-D gap    :  0.0000000000e+00
HiGHS run time      :          0.00
Solving linear system to compute dual ray
(has_dual_ray, dual_ray_value) = (Base.RefValue{Int32}(1), [0.0, 1.0, -1.0, -1.0])
Coefficient ranges:
  Matrix [1e+00, 3e+00]
  Cost   [1e+00, 1e+00]
  Bound  [1e+00, 1e+00]
  RHS    [1e+00, 1e+00]
Solving LP without presolve, or with basis, or unconstrained
Using EKK dual simplex solver - serial
  Iteration        Objective     Infeasibilities num(sum)
          0    -2.9999961628e+00 Pr: 1(3) 0s
          1     9.7888364081e-13 0s
Model status        : Infeasible
Simplex   iterations: 1
Objective value     :  0.0000000000e+00
Relative P-D gap    :  9.7888364081e-13
HiGHS run time      :          0.00
Copying known dual ray
(has_dual_ray, dual_ray_value) = (Base.RefValue{Int32}(1), [0.0, 1.0, -1.0, -1.0])
(Base.RefValue{Int32}(1), [0.0, 1.0, -1.0, -1.0])

The second call to Highs_getDualRay says Copying known dual ray but the dual ray is not the same because we have modified the problem.

[odow]

Here's a C test:

void test_dualRay_twice() {
    void* highs = Highs_create();
    int ret;
    double INF = Highs_getInfinity(highs);
    ret = Highs_changeObjectiveOffset(highs, 0.0);
    assert(ret == 0);
    ret = Highs_setStringOptionValue(highs, "presolve", "off");
    assert(ret == 0);
    ret = Highs_addCol(highs, 0.0, 0.0, 0.0, 0, NULL, NULL);
    assert(ret == 0);
    ret = Highs_addCol(highs, 0.0, 0.0, 0.0, 0, NULL, NULL);
    assert(ret == 0);
    ret = Highs_addCol(highs, -1.0, 0.0, INF, 0, NULL, NULL);
    assert(ret == 0);
    ret = Highs_addCol(highs, -1.0, 0.0, INF, 0, NULL, NULL);
    assert(ret == 0);
    int index[2] = {2, 3};
    double value[2] = {1.0, -1.0};
    ret = Highs_addRow(highs, 0.0, 0.0, 2, index, value);
    assert(ret == 0);
    index[0] = 2;    index[1] = 3;
    value[0] = 1.0;  value[1] = 1.0;
    ret = Highs_addRow(highs, 1.0, INF, 2, index, value);
    assert(ret == 0);
    index[0] = 0;    index[1] = 2;
    value[0] = -2.0; value[1] = 1.0;
    ret = Highs_addRow(highs, -INF, 0.0, 2, index, value);
    assert(ret == 0);
    index[0] = 1;    index[1] = 3;
    value[0] = -3.0; value[1] = 1.0;
    ret = Highs_addRow(highs, -INF, 0.0, 2, index, value);
    assert(ret == 0);
    ret = Highs_run(highs);
    assert(ret == 0);
    int has_dual_ray = 0;
    double dual_ray_value[4] = {0.0, 0.0, 0.0, 0.0};
    ret = Highs_getDualRay(highs, &has_dual_ray, dual_ray_value);
    assert(ret == 0);
    assertIntValuesEqual("has_dual_ray", has_dual_ray, 1);
    assertDoubleValuesEqual("dual_ray_value[0]", dual_ray_value[0], 0.0);
    assertDoubleValuesEqual("dual_ray_value[1]", dual_ray_value[1], 1.0);
    assertDoubleValuesEqual("dual_ray_value[2]", dual_ray_value[2], -1.0);
    assertDoubleValuesEqual("dual_ray_value[3]", dual_ray_value[3], -1.0);
    ret = Highs_changeColBounds(highs, 1, 1.0, 1.0);
    assert(ret == 0);
    ret = Highs_run(highs);
    assert(ret == 0);
    ret = Highs_getDualRay(highs, &has_dual_ray, dual_ray_value);
    assert(ret == 0);
    assertIntValuesEqual("has_dual_ray", has_dual_ray, 1);
    assertDoubleValuesEqual("dual_ray_value[0]", dual_ray_value[0], 1.0);
    assertDoubleValuesEqual("dual_ray_value[1]", dual_ray_value[1], 1.0);
    assertDoubleValuesEqual("dual_ray_value[2]", dual_ray_value[2], -2.0);
    assertDoubleValuesEqual("dual_ray_value[3]", dual_ray_value[3], 0.0);
    Highs_destroy(highs);
    return;
}

[odow]

I don't really understand the implications of what I'm doing, but this fixed it:

diff --git a/src/simplex/HEkk.cpp b/src/simplex/HEkk.cpp
index ff7034006..412af8f5e 100644
--- a/src/simplex/HEkk.cpp
+++ b/src/simplex/HEkk.cpp
@@ -311,6 +311,17 @@ void HEkk::invalidateBasis() {
   this->invalidateBasisArtifacts();
 }
 
+void HEkk::invalidateRays() {
+  this->status_.has_dual_ray = false;
+  this->status_.has_primal_ray = false;
+  this->info_.dual_ray_row_ = -1;
+  this->info_.dual_ray_sign_ = -1;
+  this->dual_ray_.clear();
+  this->info_.primal_ray_col_ = -1;
+  this->info_.primal_ray_sign_ = -1;
+  this->primal_ray_.clear();
+}
+
 void HEkk::invalidateBasisArtifacts() {
   // Invalidate the artifacts of the basis of the simplex LP
   this->status_.has_ar_matrix = false;
@@ -321,14 +332,7 @@ void HEkk::invalidateBasisArtifacts() {
   this->status_.has_dual_objective_value = false;
   this->status_.has_primal_objective_value = false;
   // Invalidate dual and primal ray data
-  this->status_.has_dual_ray = false;
-  this->status_.has_primal_ray = false;
-  this->info_.dual_ray_row_ = -1;
-  this->info_.dual_ray_sign_ = -1;
-  this->dual_ray_.clear();
-  this->info_.primal_ray_col_ = -1;
-  this->info_.primal_ray_sign_ = -1;
-  this->primal_ray_.clear();
+  this->invalidateRays();
 }
 
 void HEkk::updateStatus(LpAction action) {
@@ -343,11 +347,13 @@ void HEkk::updateStatus(LpAction action) {
       this->status_.has_fresh_rebuild = false;
       this->status_.has_dual_objective_value = false;
       this->status_.has_primal_objective_value = false;
+      this->invalidateRays();
       break;
     case LpAction::kNewBounds:
       this->status_.has_fresh_rebuild = false;
       this->status_.has_dual_objective_value = false;
       this->status_.has_primal_objective_value = false;
+      this->invalidateRays();
       break;
     case LpAction::kNewBasis:
       this->invalidateBasis();
diff --git a/src/simplex/HEkk.h b/src/simplex/HEkk.h
index fd270f6f3..179b567b6 100644
--- a/src/simplex/HEkk.h
+++ b/src/simplex/HEkk.h
@@ -77,6 +77,7 @@ class HEkk {
   void invalidate();
   void invalidateBasisMatrix();
   void invalidateBasis();
+  void invalidateRays();
   void invalidateBasisArtifacts();
 
   void updateStatus(LpAction action);