From a155ebb4b910c73cf84d7fb33af791d7732721da Mon Sep 17 00:00:00 2001
From: Nicola Di Cicco <93935338+nicoladicicco@users.noreply.github.com>
Date: Fri, 19 Jul 2024 10:36:11 +0200
Subject: [PATCH 1/7] Use coordinate descent for continuous variables

---
 Project.toml       |  3 ++-
 src/solver.jl      | 63 +++++++++++++++++++++++++++++++++++++++++++---
 test/raw_solver.jl | 44 ++++++++++++++++++++++++++++++++
 3 files changed, 105 insertions(+), 5 deletions(-)

diff --git a/Project.toml b/Project.toml
index 8cf79d7..16fe459 100644
--- a/Project.toml
+++ b/Project.toml
@@ -30,8 +30,9 @@ julia = "1.6"
 
 [extras]
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
+Intervals = "d8418881-c3e1-53bb-8760-2df7ec849ed5"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 TestItemRunner = "f8b46487-2199-4994-9208-9a1283c18c0a"
 
 [targets]
-test = ["Aqua", "Test", "TestItemRunner"]
+test = ["Aqua", "Intervals", "Test", "TestItemRunner"]
diff --git a/src/solver.jl b/src/solver.jl
index 7aa2ea4..b020818 100644
--- a/src/solver.jl
+++ b/src/solver.jl
@@ -336,6 +336,56 @@ function _move!(s, x::Int, dim::Int = 0)
     return best_values, best_swap, tabu
 end
 
+"""
+    armijo_line_search(f, x, d, fx; α0 = 1.0, β = 0.5, c = 1e-4)
+
+Determines the optimal step size of a line search algorithm via the Armijo condition.
+# Arguments:
+- `f`: a function to minimize
+- `x`: selected variable id
+- `d`: descent direction (e.g., negative gradient)
+- `fx`: value of f at `x`
+- `α0`: initial step size
+- `β`: step size reduction factor
+- `c`: Armijo condition constant
+"""
+function armijo_line_search(f, x, d, fx; α0 = 1.0, β = 0.5, c = 1e-4)
+    α = α0
+    while f(x + α*d) > fx + c*α*d*fx
+        α *= β
+    end
+    return α
+end
+
+"""
+    _coordinate_descent!(s, x)
+
+Runs an iteration of coordinate descent over axis "x".
+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)
+
+    function f(val)
+        _value!(s, x, val)
+        _compute!(s)
+        return get_error(s)
+    end
+
+    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)
+    new_error = f(new_value)
+    
+    if new_error < current_error
+        current_value = new_value
+    end
+end
+
 """
     _step!(s)
 
@@ -345,10 +395,15 @@ function _step!(s)
     # select worst variables
     x = _select_worse(s)
     _verbose(s, "Selected x = $x")
-
-    # Local move (change the value of the selected variable)
-    best_values, best_swap, tabu = _move!(s, x)
-    # _compute!(s)
+    
+    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)
+    end
 
     # If local move is bad (tabu), then try permutation
     if tabu
diff --git a/test/raw_solver.jl b/test/raw_solver.jl
index d20658b..d7e4d63 100644
--- a/test/raw_solver.jl
+++ b/test/raw_solver.jl
@@ -83,6 +83,39 @@ function sudoku(n; start = nothing)
     return m
 end
 
+function chemical_equilibrium(A, B, C)
+    m = model(; kind = :equilibrium)
+
+    N = length(C)
+    M = length(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
+    )
+    
+    # Add constraints
+    for i in 1:M
+        constraint!(m, conserve(i), 1:N)
+    end
+
+    # computes the total energy freed by the reaction
+    free_energy = x -> sum(j -> x[j] * (C[j] + log(x[j] / sum(x))))
+
+    objective!(m, free_energy)
+
+    return m
+end
+
 @testset "Raw solver: internals" begin
     models = [
         sudoku(2)
@@ -203,3 +236,14 @@ end
     @info "Sol (vals): $(!isnothing(best_value(s)) ? best_values(s) : nothing)"
     @info time_info(s)
 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]
+    B = [20.0, 30.0, 25.0]
+    C = [-10.0, -8.0, -6.0]
+    m = chemical_equilibrium(A, B, C)
+    s = solver(m; options = Options(print_level = :minimal))
+    solve!(s)
+    display(solution(s))
+    display(s.time_stamps)
+end

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 2/7] 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)

From c718119b0795a97cea55e5f9468bb15452503b62 Mon Sep 17 00:00:00 2001
From: Nicola Di Cicco <93935338+nicoladicicco@users.noreply.github.com>
Date: Fri, 19 Jul 2024 18:57:42 +0200
Subject: [PATCH 3/7] fix issue in mutating solver's state

---
 src/solver.jl | 16 +++++++---------
 1 file changed, 7 insertions(+), 9 deletions(-)

diff --git a/src/solver.jl b/src/solver.jl
index ebb233f..c72c408 100644
--- a/src/solver.jl
+++ b/src/solver.jl
@@ -349,7 +349,7 @@ 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, xmin, xmax; α0 = 1.0, β = 0.9, c = 1e-4)
+function armijo_line_search(f, x, d, fx, xmin, xmax; α0 = 1.0, β = 0.5, c = 1e-4)
     α = α0
     while true
         new_x = x + α * d
@@ -378,8 +378,6 @@ function _coordinate_descent_move!(s, x)
     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)
@@ -395,13 +393,13 @@ function _coordinate_descent_move!(s, x)
     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
-        best_values = [new_value]
-        tabu = false
-    end
+    best_values = [new_value]
+
+    # revert to the original state
+    _value!(s, x, current_value)
+    _compute!(s)
 
-    return best_values, [x], tabu
+    return best_values, [x], false
 end
 
 """

From af5e88447c9a2cbceb789bcda778f4a8e21e4be1 Mon Sep 17 00:00:00 2001
From: Azzaare <jf@baffier.fr>
Date: Mon, 22 Jul 2024 11:52:26 +0000
Subject: [PATCH 4/7] Fix the test for a domain being continuous

---
 src/LocalSearchSolvers.jl |  1 +
 src/solver.jl             | 11 +++++------
 2 files changed, 6 insertions(+), 6 deletions(-)

diff --git a/src/LocalSearchSolvers.jl b/src/LocalSearchSolvers.jl
index 6ae8493..67c2c1a 100644
--- a/src/LocalSearchSolvers.jl
+++ b/src/LocalSearchSolvers.jl
@@ -8,6 +8,7 @@ using Dictionaries
 using Distributed
 using JSON
 using Lazy
+using TestItems
 
 # Exports internal
 export constraint!, variable!, objective!, add!, add_var_to_cons!, add_value!
diff --git a/src/solver.jl b/src/solver.jl
index c72c408..0681df6 100644
--- a/src/solver.jl
+++ b/src/solver.jl
@@ -372,7 +372,7 @@ end
 
 Runs an iteration of coordinate descent over axis "x".
 The derivative is (temporarily?) computed via finite difference.
-The step size is determined via the Armijo condition for line search. 
+The step size is determined via the Armijo condition for line search.
 """
 function _coordinate_descent_move!(s, x)
     current_value = _value(s, x)
@@ -412,13 +412,12 @@ function _step!(s)
     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
+    if get_domain(s, x) isa ContinuousDomain
         # We perform coordinate descent over the variable axis
         best_values, best_swap, tabu = _coordinate_descent_move!(s, x)
+    else
+        # Local move (change the value of the selected variable)
+        best_values, best_swap, tabu = _move!(s, x)
     end
 
     # If local move is bad (tabu), then try permutation

From 1f928dee47c76e111bbc09b18d228a5192b61340 Mon Sep 17 00:00:00 2001
From: Nicola Di Cicco <93935338+nicoladicicco@users.noreply.github.com>
Date: Mon, 22 Jul 2024 21:10:29 +0200
Subject: [PATCH 5/7] add test for a simple quadratic function

---
 src/solver.jl      |  3 ++-
 test/raw_solver.jl | 26 ++++++++++++++++++++++++--
 test/runtests.jl   |  1 +
 3 files changed, 27 insertions(+), 3 deletions(-)

diff --git a/src/solver.jl b/src/solver.jl
index 0681df6..65de885 100644
--- a/src/solver.jl
+++ b/src/solver.jl
@@ -368,7 +368,7 @@ function armijo_line_search(f, x, d, fx, xmin, xmax; α0 = 1.0, β = 0.5, c = 1e
 end
 
 """
-    _coordinate_descent!(s, x)
+    _coordinate_descent_move!(s, x)
 
 Runs an iteration of coordinate descent over axis "x".
 The derivative is (temporarily?) computed via finite difference.
@@ -412,6 +412,7 @@ function _step!(s)
     x = _select_worse(s)
     _verbose(s, "Selected x = $x")
 
+    # TODO: this if statement is currently not working, every variable is treated as integer 
     if get_domain(s, x) isa ContinuousDomain
         # We perform coordinate descent over the variable axis
         best_values, best_swap, tabu = _coordinate_descent_move!(s, x)
diff --git a/test/raw_solver.jl b/test/raw_solver.jl
index ca1e3af..02002c5 100644
--- a/test/raw_solver.jl
+++ b/test/raw_solver.jl
@@ -1,5 +1,3 @@
-using Intervals
-
 function mincut(graph; source, sink, interdiction = 0)
     m = model(; kind = :cut)
     n = size(graph, 1)
@@ -117,6 +115,20 @@ function chemical_equilibrium(A, B, C)
     return m
 end
 
+function sum_squares(n)
+    m = model(; kind = :sum_squares)
+
+    d = domain(-10.0 .. 10.0)
+
+    foreach(_ -> variable!(m, d), 1:n)
+
+    ss = x -> sum(j -> j * x[j] * x[j], 1:n)
+
+    objective!(m, ss)
+
+    return m
+end
+
 @testset "Raw solver: internals" begin
     models = [
         sudoku(2)
@@ -248,3 +260,13 @@ end
     display(solution(s))
     display(s.time_stamps)
 end
+
+@testset "Raw solver: sum squares" begin
+    # Sanity check: simple quadratic function with a trivial minimum at 0
+    n = 10
+    m = sum_squares(n)
+    s = solver(m; options = Options(print_level = :minimal))
+    solve!(s)
+    display(solution(s))
+    display(s.time_stamps)
+end
\ No newline at end of file
diff --git a/test/runtests.jl b/test/runtests.jl
index 1e5a532..1ee8e2f 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -4,6 +4,7 @@ import ConstraintDomains
 import CompositionalNetworks
 @everywhere using Constraints
 using Dictionaries
+using Intervals
 @everywhere using LocalSearchSolvers
 using Test
 using TestItemRunner

From bfe33f5ffc94fbec0e9b95005c19019666838f82 Mon Sep 17 00:00:00 2001
From: Nicola Di Cicco <93935338+nicoladicicco@users.noreply.github.com>
Date: Tue, 23 Jul 2024 13:05:20 +0200
Subject: [PATCH 6/7] fix domain check not triggering

---
 src/solver.jl | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/src/solver.jl b/src/solver.jl
index 65de885..1f929ab 100644
--- a/src/solver.jl
+++ b/src/solver.jl
@@ -412,8 +412,7 @@ function _step!(s)
     x = _select_worse(s)
     _verbose(s, "Selected x = $x")
 
-    # TODO: this if statement is currently not working, every variable is treated as integer 
-    if get_domain(s, x) isa ContinuousDomain
+    if typeof(get_variable(s, x).domain) <: ContinuousDomain
         # We perform coordinate descent over the variable axis
         best_values, best_swap, tabu = _coordinate_descent_move!(s, x)
     else

From 2fe68121bf7a430f73331b621daad0814c0a659f Mon Sep 17 00:00:00 2001
From: Azzaare <jf@baffier.fr>
Date: Thu, 1 Aug 2024 11:15:05 +0900
Subject: [PATCH 7/7] Fix format in tests

---
 test/raw_solver.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test/raw_solver.jl b/test/raw_solver.jl
index 02002c5..2d9dc0e 100644
--- a/test/raw_solver.jl
+++ b/test/raw_solver.jl
@@ -269,4 +269,4 @@ end
     solve!(s)
     display(solution(s))
     display(s.time_stamps)
-end
\ No newline at end of file
+end