From 584d9692b996051553918e32964e73df8e89ab8e Mon Sep 17 00:00:00 2001
From: Nicola Di Cicco <93935338+nicoladicicco@users.noreply.github.com>
Date: Fri, 19 Jul 2024 15:41:05 +0200
Subject: [PATCH] fix tests and improve algorithm

---
 src/solver.jl      | 38 ++++++++++++++++++++++++++++----------
 test/raw_solver.jl | 19 ++++++++++---------
 2 files changed, 38 insertions(+), 19 deletions(-)

diff --git a/src/solver.jl b/src/solver.jl
index b020818..ebb233f 100644
--- a/src/solver.jl
+++ b/src/solver.jl
@@ -349,9 +349,19 @@ Determines the optimal step size of a line search algorithm via the Armijo condi
 - `β`: step size reduction factor
 - `c`: Armijo condition constant
 """
-function armijo_line_search(f, x, d, fx; α0 = 1.0, β = 0.5, c = 1e-4)
+function armijo_line_search(f, x, d, fx, xmin, xmax; α0 = 1.0, β = 0.9, c = 1e-4)
     α = α0
-    while f(x + α*d) > fx + c*α*d*fx
+    while true
+        new_x = x + α * d
+        # Ensure new_x is within [xmin, xmax]
+        if new_x < xmin
+            new_x = xmin
+        elseif new_x > xmax
+            new_x = xmax
+        end
+        if f(new_x) <= fx + c * α * d * fx || α < 1e-8
+            break
+        end
         α *= β
     end
     return α
@@ -365,8 +375,11 @@ The derivative is (temporarily?) computed via finite difference.
 The step size is determined via the Armijo condition for line search. 
 """
 function _coordinate_descent_move!(s, x)
-    domain = get_variable(s, x).domain
     current_value = _value(s, x)
+    xmin = minimum(first, get_domain(s, x))
+    xmax = maximum(last, get_domain(s, x))
+    best_values = [current_value]
+    tabu = true
 
     function f(val)
         _value!(s, x, val)
@@ -376,14 +389,19 @@ function _coordinate_descent_move!(s, x)
 
     current_error = f(current_value)
     grad = (f(current_value + 1e-6) - f(current_value - 1e-6)) / (2e-6)
-    
-    α = armijo_line_search(f, current_value, -grad, current_error)
-    new_value = clamp(current_value - α * grad, domain.lb, domain.ub)
+
+    α = armijo_line_search(f, current_value, -grad, current_error, xmin, xmax)
+
+    new_value = clamp(current_value - α * grad,
+        minimum(first, get_domain(s, x)), maximum(last, get_domain(s, x)))
+
     new_error = f(new_value)
-    
     if new_error < current_error
-        current_value = new_value
+        best_values = [new_value]
+        tabu = false
     end
+
+    return best_values, [x], tabu
 end
 
 """
@@ -395,14 +413,14 @@ function _step!(s)
     # select worst variables
     x = _select_worse(s)
     _verbose(s, "Selected x = $x")
-    
+
     if _value(s, x) isa Int
         # Local move (change the value of the selected variable)
         best_values, best_swap, tabu = _move!(s, x)
         # _compute!(s)
     else
         # We perform coordinate descent over the variable axis
-        _coordinate_descent_move!(s, x)
+        best_values, best_swap, tabu = _coordinate_descent_move!(s, x)
     end
 
     # If local move is bad (tabu), then try permutation
diff --git a/test/raw_solver.jl b/test/raw_solver.jl
index d7e4d63..ca1e3af 100644
--- a/test/raw_solver.jl
+++ b/test/raw_solver.jl
@@ -1,3 +1,5 @@
+using Intervals
+
 function mincut(graph; source, sink, interdiction = 0)
     m = model(; kind = :cut)
     n = size(graph, 1)
@@ -89,20 +91,19 @@ function chemical_equilibrium(A, B, C)
     N = length(C)
     M = length(B)
 
-    d = domain(0..maximum(B))
-    
+    d = domain(0 .. maximum(B))
+
     # Add variables, number of moles per compound
 
     foreach(_ -> variable!(m, d), 1:N)
 
     # mass_conservation function
-    conserve = i -> (x ->
-        begin
-            δ = abs(sum(A[:, i] .* x) - B[i])
-            return δ ≤ 1.e-6 ? 0. : δ
-        end
+    conserve = i -> (x -> begin
+        δ = abs(sum(A[:, i] .* x) - B[i])
+        return δ ≤ 1.e-6 ? 0.0 : δ
+    end
     )
-    
+
     # Add constraints
     for i in 1:M
         constraint!(m, conserve(i), 1:N)
@@ -238,7 +239,7 @@ end
 end
 
 @testset "Raw solver: chemical equilibrium" begin
-    A = [2.0 1.0 0.0; 6.0 2.0 1.0; 1.0 2.0 4.0]
+    A = [2.0 1.0 3.0; 6.0 2.0 1.0; 1.0 2.0 4.0]
     B = [20.0, 30.0, 25.0]
     C = [-10.0, -8.0, -6.0]
     m = chemical_equilibrium(A, B, C)