From 8ba47f9680d9353f293924941ccf74356d3e911c Mon Sep 17 00:00:00 2001
From: tinatorabi <tntorabii@gmail.com>
Date: Sat, 19 Oct 2024 13:19:38 -0700
Subject: [PATCH 1/5] inplace QR

---
 Project.toml                  |  3 +-
 src/BPDual.jl                 | 44 ++++++++++++----------
 src/helpers.jl                | 69 ++++++++++++++++++++++++++++++++---
 test/test_recover_decaying.jl |  2 +-
 4 files changed, 91 insertions(+), 27 deletions(-)

diff --git a/Project.toml b/Project.toml
index d9c68ce..d099594 100644
--- a/Project.toml
+++ b/Project.toml
@@ -17,11 +17,10 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 DataFrames = "1.6"
 LinearAlgebra = "1.9, 1.10"
 LinearOperators = "2.8"
+Logging = "1.9"
 Printf = "1.9, 1.10"
-QRupdate = "1.0"
 Random = "1.9, 1.10"
 SparseArrays = "1.9, 1.10"
-Logging = "1.9"
 julia = "1.9, 1.10"
 
 [extras]
diff --git a/src/BPDual.jl b/src/BPDual.jl
index e8690f8..99c8b18 100644
--- a/src/BPDual.jl
+++ b/src/BPDual.jl
@@ -76,7 +76,7 @@ function bpdual(
     active::Union{Nothing, Vector{Int}} = nothing, ## maybe write as struct later
     state::Union{Nothing, Vector{Int}} = nothing,
     y::Union{Nothing, Vector{Float64}} = nothing,
-    S = Matrix{Float64}(undef, size(A, 1), 0),
+    S = Matrix{Float64}(undef, size(A, 1), size(A, 2)),
     R::Union{Nothing, Matrix{Float64}} = nothing,
     loglevel::Int = 1,
     coldstart::Bool = true,
@@ -110,8 +110,8 @@ function bpdual(
         active = Vector{Int}([])
         state = Vector{Int}(zeros(Int, n))
         y = Vector{Float64}(zeros(Float64, m))
-        S = Matrix{Float64}(zeros(Float64, m, 0))
-        R = Matrix{Float64}(zeros(Float64, 0, 0))
+        S = Matrix{Float64}(zeros(Float64, m, n))
+        R = Matrix{Float64}(zeros(Float64, n,n))
     end
 
     if homotopy
@@ -171,7 +171,7 @@ function bpdual(
     numtrim = 0
     nprodA = 0
     nprodAt = 0
-
+    current_R_size = 0
     # ------------------------------------------------------------------
     # Cold/warm-start initialization.
     # ------------------------------------------------------------------
@@ -211,12 +211,12 @@ function bpdual(
         sL, sU = infeasibilities(bl, bu, z)
         g = b - λ*y  # Steepest-descent direction
 
-        if isempty(R)
-            condS = 1
-        else
-            rmin = minimum(diag(R))
-            rmax = maximum(diag(R))
+        if current_R_size > 0
+            rmin = minimum(diag(R[1:current_R_size,1:current_R_size]))
+            rmax = maximum(diag(R[1:current_R_size,1:current_R_size]))
             condS = rmax / rmin
+        else
+            condS = 1  # Set to a default condition number if `R` is empty.
         end
 
         if condS > 1e+10
@@ -225,11 +225,10 @@ function bpdual(
             npad = size(S, 2) - size(x, 1)
             x = [x; zeros(npad)]
         else
-            dx, dy = newtonstep(S, R, g, x, λ)
+            dx, dy = newtonstep(S[:,1:current_R_size], R[1:current_R_size,1:current_R_size], g, x, λ)
             x .+= dx
         end
-
-        r = b - S*x
+        r = b - S[:,1:current_R_size]*x
 
         # Print to log.
         yNorm = norm(y, 2)
@@ -273,7 +272,7 @@ function bpdual(
                 @info "\nOptimal solution found. Trimming multipliers..."
             end
             g = b - λin*y
-            trimx(x, S, R, active, state, g, b, λ, feaTol, optTol, loglevel)
+            trimx!(x, S, R, active, state, g, b, λ, feaTol, optTol, current_R_size)
             numtrim = nact - length(active)
             nact = length(active)
         end
@@ -286,9 +285,9 @@ function bpdual(
         # New iteration starts here.
         itn += 1
         p = q = 0
-
+        
         if homotopy
-            x, dy, dz, step, λ, p = htpynewlam(active, state, A, R, S, x, y, sL, sU, λ, lamFinal)
+            x, dy, dz, step, λ, p = htpynewlam(active, state, A, R[1:current_R_size,1:current_R_size], S[:,1:current_R_size], x, y, sL, sU, λ, lamFinal)
             nprodAt += 1
         else
             if norm(dy, Inf) < eps()        
@@ -339,8 +338,10 @@ function bpdual(
             a = A * zerovec
             nprodA += 1
             zerovec[p] = 0
-            R = qraddcol(S, R, a)
-            S = [S a]
+            qraddcol!(S, R, a, current_R_size)
+            # R = qraddcol(S, R, a)
+            current_R_size += 1
+            # S = [S a]
             push!(active, p)
             push!(x, 0)
         else
@@ -358,10 +359,11 @@ function bpdual(
                 _, qa = findmax(abs.(x .* dropa))
                 q = active[qa]
                 state[q] = 0
-                S = S[:, 1:nact .!= qa]
+                S[:, 1:nact][:, 1:nact .== qa] .= 0
                 deleteat!(active, qa)
                 deleteat!(x, qa)
-                R = qrdelcol(R, qa)
+                qrdelcol!(S, R, qa)
+                current_R_size -=1
             else
                 eFlag = :EXIT_OPTIMAL
             end
@@ -391,3 +393,7 @@ function bpdual(
     end
     return tracer
 end
+
+
+
+
diff --git a/src/helpers.jl b/src/helpers.jl
index f13df18..1d9d789 100644
--- a/src/helpers.jl
+++ b/src/helpers.jl
@@ -146,7 +146,61 @@ end
 # products with A. (Implicitly we're assuming that A has a small norm.)
 # ----------------------------------------------------------------------
 
-function trimx(x,S,R,active,state,g,b,λ,featol,opttol,loglevel)
+# function trimx(x,S,R,active,state,g,b,λ,featol,opttol,loglevel)
+
+#     k = 0
+#     nact = length(active)
+#     xabs = abs.(x)
+#     xmin, qa = findmin(xabs)
+#     gNorm = norm(g,Inf)
+
+#     while xmin < opttol
+
+#         e = sign.(x.*(xabs .> opttol)) # Signs of significant multipliers
+#         q = active[qa] # Index of the corresponding constraint
+#         a = S[:,qa]    # Save the col from S in case we need to add it back. 
+#         xsmall = x[qa] # Value of candidate multiplier
+
+#         # Trim quantities related to the small multiplier.
+#         deleteat!(e, qa)
+#         S = [S[:, 1:qa-1] S[:, qa+1:end]]
+#         deleteat!(active, qa)
+#         # R = qrdelcol(R, qa)
+
+#         R = qrdelcol!(S,R, qa)
+
+
+#         # Recompute the remaining multipliers and their signs.
+#         x, dy = csne(R, S, vec(g))           # min ||g - Sx||_2
+#         xabs = abs.(x)
+#         et = sign.(x.*(xabs .> opttol))
+#         dyNorm = norm(dy,Inf)/λ        # dy = (g - Sx) / lambda
+    
+#         # Check if the trimmed active set is still optimal
+#         if any( et .!= e ) || (dyNorm / max(1,gNorm) > featol)
+#             R = qraddcol(S,R,a)
+#             S = [S a]
+#             push!(active, q)
+#             x = csne(R, S,vec(g))
+#             break
+#         end
+
+#         # The trimmed x is still optimal.
+#         k += 1
+#         nact -= 1
+#         state[q] = 0               # Mark this constraint as free.
+#         rNorm = norm(b - S*x)
+#         xNorm = norm(x,1)
+        
+#         # Grab the next canddate multiplier.
+#         xmin, qa = findmin(xabs)
+#     end  
+#     return (x,S,R,active,state)
+# end
+
+
+
+function trimx!(x,S,R,active,state,g,b,λ,featol,opttol, current_R_size)
 
     k = 0
     nact = length(active)
@@ -165,7 +219,10 @@ function trimx(x,S,R,active,state,g,b,λ,featol,opttol,loglevel)
         deleteat!(e, qa)
         S = [S[:, 1:qa-1] S[:, qa+1:end]]
         deleteat!(active, qa)
-        R = qrdelcol(R, qa)
+        # R = qrdelcol(R, qa)
+
+        qrdelcol!(S,R, qa)
+        current_R_size -= 1
 
         # Recompute the remaining multipliers and their signs.
         x, dy = csne(R, S, vec(g))           # min ||g - Sx||_2
@@ -175,7 +232,8 @@ function trimx(x,S,R,active,state,g,b,λ,featol,opttol,loglevel)
     
         # Check if the trimmed active set is still optimal
         if any( et .!= e ) || (dyNorm / max(1,gNorm) > featol)
-            R = qraddcol(S,R,a)
+            qraddcol!(S,R,a, current_R_size)
+            current_R_size +=1
             S = [S a]
             push!(active, q)
             x = csne(R, S,vec(g))
@@ -192,7 +250,7 @@ function trimx(x,S,R,active,state,g,b,λ,featol,opttol,loglevel)
         # Grab the next canddate multiplier.
         xmin, qa = findmin(xabs)
     end  
-    return (x,S,R,active,state)
+    return (x,S,R,active,state, current_R_size)
 end
 
 
@@ -221,7 +279,8 @@ function triminf(active::Vector, state::Vector, S::Matrix, R::Matrix,
         nact = nact - 1 
         S = S[:,1:size(S,2) .!= qa] # Delete column from S
         deleteat!(active, qa)   # Delete index from active set
-        R = qrdelcol(R, qa)     # Recompute new QR factorization
+        # R = qrdelcol(R, qa)     # Recompute new QR factorization
+        R = qrdelcol(R, qa)
         state[q] = 0            # Mark constraint as free
         x = csne(R, S, g)       # Recompute multipliers
 
diff --git a/test/test_recover_decaying.jl b/test/test_recover_decaying.jl
index 687a632..ad43f98 100644
--- a/test/test_recover_decaying.jl
+++ b/test/test_recover_decaying.jl
@@ -22,7 +22,7 @@ function test_recover_decaying()
     bu = +ones(n)
 
     # Solve the basis pursuit problem
-    tracer = asp_homotopy(A, b, min_lambda = 0.0, itnMax = 400, loglevel =0) 
+    tracer = asp_homotopy(A, b, min_lambda = 0.0, actMax = 400, loglevel =0) 
     xs, λ = tracer[end]
 
     cumulative_norm = cumsum(abs.(x[sortperm(abs.(x), rev=true)]))

From 7ade4eb42b29667eb28be8c2101ededa84587b40 Mon Sep 17 00:00:00 2001
From: Tina Torabi <tntorabii@gmail.com>
Date: Sat, 2 Nov 2024 21:40:36 -0700
Subject: [PATCH 2/5] up

---
 src/omp.jl | 27 ++++++++++++---------------
 1 file changed, 12 insertions(+), 15 deletions(-)

diff --git a/src/omp.jl b/src/omp.jl
index 14a7f32..ec9d5eb 100644
--- a/src/omp.jl
+++ b/src/omp.jl
@@ -108,14 +108,10 @@ function asp_omp(
 
     if eFlag != :EXIT_UNKNOWN || active === nothing
         active = Vector{Int}([])
-    end
-    if state === nothing
-        state = zeros(Int, n)
-    end
-    if R === nothing
         R = Matrix{Float64}(undef, 0, 0)
+        S = Matrix{Float64}(undef, m, 0)
+        state = zeros(Int, n)
     end
-
     @info @sprintf("%4s  %8s %12s %12s %12s", "Itn", "Var", "λ", "rNorm", "xNorm")
 
 
@@ -129,7 +125,7 @@ function asp_omp(
             nprodAt += 1
             zmax = norm(z, Inf)
         else
-            x,y = csne(R, S, vec(b))
+            x, r = csne(R, S, b)
             if norm(x, Inf) > 1e12
                 eFlag = :EXIT_SINGULAR_LS
                 break
@@ -181,21 +177,22 @@ function asp_omp(
         zerovec[p] = 0
 
         R = qraddcol(S, R, a)  # Update R
-        S = hcat(S, a)  # Expand S, active
+        S = [S a]
 
         push!(tracer.iteration, itn)
         push!(tracer.lambda, zmax)
-        sparse_x_full = SparseVector(n, copy(active), copy(x))
+        sparse_x_full = spzeros(n)
+        sparse_x_full[copy(active)] = copy(x)  
         push!(tracer.solution, copy(sparse_x_full))
-        
         push!(active, p)
 
     end #while true
 
-    push!(tracer.iteration, itn)
-    push!(tracer.lambda, zmax)
-    sparse_x_full = SparseVector(n, copy(active), copy(x))
-    push!(tracer.solution, copy(sparse_x_full))
+    # push!(tracer.iteration, itn)
+    # push!(tracer.lambda, zmax)
+    # sparse_x_full = spzeros(n)
+    # sparse_x_full[copy(active)] = copy(x)  
+    # push!(tracer.solution, copy(sparse_x_full))
 
     tottime = time() - time0
     if loglevel > 0
@@ -206,4 +203,4 @@ function asp_omp(
         @info "\n"
     end
     return tracer
-end
\ No newline at end of file
+end

From 36d1b17b67c079cd2515b63ab8552fe774f52715 Mon Sep 17 00:00:00 2001
From: Tina Torabi <tntorabii@gmail.com>
Date: Sat, 2 Nov 2024 21:56:08 -0700
Subject: [PATCH 3/5] up

---
 src/omp.jl | 8 ++++++--
 1 file changed, 6 insertions(+), 2 deletions(-)

diff --git a/src/omp.jl b/src/omp.jl
index ec9d5eb..65e04b5 100644
--- a/src/omp.jl
+++ b/src/omp.jl
@@ -31,6 +31,7 @@ Base.lastindex(t::OMPTracer) = lastindex(t.iteration)
     Ax = b.
     ```
 """
+
 function asp_omp(
     A::Union{AbstractMatrix, AbstractLinearOperator},
     b::Vector,
@@ -84,7 +85,8 @@ function asp_omp(
         :EXIT_LAMBDA => "Reached minimum value of lambda",
         :EXIT_RHS_ZERO => "b = 0. The solution is x = 0",
         :EXIT_UNCONSTRAINED => "Unconstrained solution r = b is optimal",
-        :EXIT_UNKNOWN => "unknown exit"
+        :EXIT_UNKNOWN => "unknown exit",
+        :EXIT_ACTMAX => "Max no. of active constraints reached"
     )
 
     itn = 0
@@ -149,8 +151,10 @@ function asp_omp(
             eFlag = :EXIT_OPTIMAL
         elseif itn >= itnMax
             eFlag = :EXIT_TOO_MANY_ITNS
+        elseif eFlag == :EXIT_UNKNOWN && length(active) >= actMax
+            eFlag = :EXIT_ACTMAX
         end
-        
+                
         if eFlag != :EXIT_UNKNOWN
             break
         end

From 5d88cee8c8c0da84ef714cea546827a4b576a61c Mon Sep 17 00:00:00 2001
From: tinatorabi <tntorabii@gmail.com>
Date: Tue, 5 Nov 2024 21:48:24 -0800
Subject: [PATCH 4/5] omp tests + removed in-place

---
 src/BPDual.jl     | 46 +++++++++++++-----------------
 src/helpers.jl    | 71 ++++-------------------------------------------
 src/omp.jl        | 61 +++++++++++++++++++++-------------------
 test/test_bpdn.jl | 31 +++++++++++++--------
 4 files changed, 78 insertions(+), 131 deletions(-)

diff --git a/src/BPDual.jl b/src/BPDual.jl
index 99c8b18..212ef9e 100644
--- a/src/BPDual.jl
+++ b/src/BPDual.jl
@@ -76,7 +76,7 @@ function bpdual(
     active::Union{Nothing, Vector{Int}} = nothing, ## maybe write as struct later
     state::Union{Nothing, Vector{Int}} = nothing,
     y::Union{Nothing, Vector{Float64}} = nothing,
-    S = Matrix{Float64}(undef, size(A, 1), size(A, 2)),
+    S = Matrix{Float64}(undef, size(A, 1), 0),
     R::Union{Nothing, Matrix{Float64}} = nothing,
     loglevel::Int = 1,
     coldstart::Bool = true,
@@ -110,8 +110,8 @@ function bpdual(
         active = Vector{Int}([])
         state = Vector{Int}(zeros(Int, n))
         y = Vector{Float64}(zeros(Float64, m))
-        S = Matrix{Float64}(zeros(Float64, m, n))
-        R = Matrix{Float64}(zeros(Float64, n,n))
+        S = Matrix{Float64}(zeros(Float64, m, 0))
+        R = Matrix{Float64}(zeros(Float64, 0, 0))
     end
 
     if homotopy
@@ -171,7 +171,7 @@ function bpdual(
     numtrim = 0
     nprodA = 0
     nprodAt = 0
-    current_R_size = 0
+
     # ------------------------------------------------------------------
     # Cold/warm-start initialization.
     # ------------------------------------------------------------------
@@ -211,12 +211,12 @@ function bpdual(
         sL, sU = infeasibilities(bl, bu, z)
         g = b - λ*y  # Steepest-descent direction
 
-        if current_R_size > 0
-            rmin = minimum(diag(R[1:current_R_size,1:current_R_size]))
-            rmax = maximum(diag(R[1:current_R_size,1:current_R_size]))
-            condS = rmax / rmin
+        if isempty(R)
+            condS = 1
         else
-            condS = 1  # Set to a default condition number if `R` is empty.
+            rmin = minimum(diag(R))
+            rmax = maximum(diag(R))
+            condS = rmax / rmin
         end
 
         if condS > 1e+10
@@ -225,10 +225,11 @@ function bpdual(
             npad = size(S, 2) - size(x, 1)
             x = [x; zeros(npad)]
         else
-            dx, dy = newtonstep(S[:,1:current_R_size], R[1:current_R_size,1:current_R_size], g, x, λ)
+            dx, dy = newtonstep(S, R, g, x, λ)
             x .+= dx
         end
-        r = b - S[:,1:current_R_size]*x
+
+        r = b - S*x
 
         # Print to log.
         yNorm = norm(y, 2)
@@ -272,7 +273,7 @@ function bpdual(
                 @info "\nOptimal solution found. Trimming multipliers..."
             end
             g = b - λin*y
-            trimx!(x, S, R, active, state, g, b, λ, feaTol, optTol, current_R_size)
+            trimx(x, S, R, active, state, g, b, λ, feaTol, optTol, loglevel)
             numtrim = nact - length(active)
             nact = length(active)
         end
@@ -285,9 +286,9 @@ function bpdual(
         # New iteration starts here.
         itn += 1
         p = q = 0
-        
+
         if homotopy
-            x, dy, dz, step, λ, p = htpynewlam(active, state, A, R[1:current_R_size,1:current_R_size], S[:,1:current_R_size], x, y, sL, sU, λ, lamFinal)
+            x, dy, dz, step, λ, p = htpynewlam(active, state, A, R, S, x, y, sL, sU, λ, lamFinal)
             nprodAt += 1
         else
             if norm(dy, Inf) < eps()        
@@ -338,10 +339,8 @@ function bpdual(
             a = A * zerovec
             nprodA += 1
             zerovec[p] = 0
-            qraddcol!(S, R, a, current_R_size)
-            # R = qraddcol(S, R, a)
-            current_R_size += 1
-            # S = [S a]
+            R = qraddcol(S, R, a)
+            S = [S a]
             push!(active, p)
             push!(x, 0)
         else
@@ -359,11 +358,10 @@ function bpdual(
                 _, qa = findmax(abs.(x .* dropa))
                 q = active[qa]
                 state[q] = 0
-                S[:, 1:nact][:, 1:nact .== qa] .= 0
+                S = S[:, 1:nact .!= qa]
                 deleteat!(active, qa)
                 deleteat!(x, qa)
-                qrdelcol!(S, R, qa)
-                current_R_size -=1
+                R = qrdelcol(R, qa)
             else
                 eFlag = :EXIT_OPTIMAL
             end
@@ -392,8 +390,4 @@ function bpdual(
         @info "Solution time (sec): $tottime"
     end
     return tracer
-end
-
-
-
-
+end
\ No newline at end of file
diff --git a/src/helpers.jl b/src/helpers.jl
index 1d9d789..a71364d 100644
--- a/src/helpers.jl
+++ b/src/helpers.jl
@@ -146,61 +146,7 @@ end
 # products with A. (Implicitly we're assuming that A has a small norm.)
 # ----------------------------------------------------------------------
 
-# function trimx(x,S,R,active,state,g,b,λ,featol,opttol,loglevel)
-
-#     k = 0
-#     nact = length(active)
-#     xabs = abs.(x)
-#     xmin, qa = findmin(xabs)
-#     gNorm = norm(g,Inf)
-
-#     while xmin < opttol
-
-#         e = sign.(x.*(xabs .> opttol)) # Signs of significant multipliers
-#         q = active[qa] # Index of the corresponding constraint
-#         a = S[:,qa]    # Save the col from S in case we need to add it back. 
-#         xsmall = x[qa] # Value of candidate multiplier
-
-#         # Trim quantities related to the small multiplier.
-#         deleteat!(e, qa)
-#         S = [S[:, 1:qa-1] S[:, qa+1:end]]
-#         deleteat!(active, qa)
-#         # R = qrdelcol(R, qa)
-
-#         R = qrdelcol!(S,R, qa)
-
-
-#         # Recompute the remaining multipliers and their signs.
-#         x, dy = csne(R, S, vec(g))           # min ||g - Sx||_2
-#         xabs = abs.(x)
-#         et = sign.(x.*(xabs .> opttol))
-#         dyNorm = norm(dy,Inf)/λ        # dy = (g - Sx) / lambda
-    
-#         # Check if the trimmed active set is still optimal
-#         if any( et .!= e ) || (dyNorm / max(1,gNorm) > featol)
-#             R = qraddcol(S,R,a)
-#             S = [S a]
-#             push!(active, q)
-#             x = csne(R, S,vec(g))
-#             break
-#         end
-
-#         # The trimmed x is still optimal.
-#         k += 1
-#         nact -= 1
-#         state[q] = 0               # Mark this constraint as free.
-#         rNorm = norm(b - S*x)
-#         xNorm = norm(x,1)
-        
-#         # Grab the next canddate multiplier.
-#         xmin, qa = findmin(xabs)
-#     end  
-#     return (x,S,R,active,state)
-# end
-
-
-
-function trimx!(x,S,R,active,state,g,b,λ,featol,opttol, current_R_size)
+function trimx(x,S,R,active,state,g,b,λ,featol,opttol,loglevel)
 
     k = 0
     nact = length(active)
@@ -219,10 +165,7 @@ function trimx!(x,S,R,active,state,g,b,λ,featol,opttol, current_R_size)
         deleteat!(e, qa)
         S = [S[:, 1:qa-1] S[:, qa+1:end]]
         deleteat!(active, qa)
-        # R = qrdelcol(R, qa)
-
-        qrdelcol!(S,R, qa)
-        current_R_size -= 1
+        R = qrdelcol(R, qa)
 
         # Recompute the remaining multipliers and their signs.
         x, dy = csne(R, S, vec(g))           # min ||g - Sx||_2
@@ -232,8 +175,7 @@ function trimx!(x,S,R,active,state,g,b,λ,featol,opttol, current_R_size)
     
         # Check if the trimmed active set is still optimal
         if any( et .!= e ) || (dyNorm / max(1,gNorm) > featol)
-            qraddcol!(S,R,a, current_R_size)
-            current_R_size +=1
+            R = qraddcol(S,R,a)
             S = [S a]
             push!(active, q)
             x = csne(R, S,vec(g))
@@ -250,7 +192,7 @@ function trimx!(x,S,R,active,state,g,b,λ,featol,opttol, current_R_size)
         # Grab the next canddate multiplier.
         xmin, qa = findmin(xabs)
     end  
-    return (x,S,R,active,state, current_R_size)
+    return (x,S,R,active,state)
 end
 
 
@@ -279,8 +221,7 @@ function triminf(active::Vector, state::Vector, S::Matrix, R::Matrix,
         nact = nact - 1 
         S = S[:,1:size(S,2) .!= qa] # Delete column from S
         deleteat!(active, qa)   # Delete index from active set
-        # R = qrdelcol(R, qa)     # Recompute new QR factorization
-        R = qrdelcol(R, qa)
+        R = qrdelcol(R, qa)     # Recompute new QR factorization
         state[q] = 0            # Mark constraint as free
         x = csne(R, S, g)       # Recompute multipliers
 
@@ -446,4 +387,4 @@ function htpynewlam(active, state, A, R, S, x, y, s1, s2, λ, lamFinal)
     end
     
     return x, dy, dz, step, λ, p
-end
+end
\ No newline at end of file
diff --git a/src/omp.jl b/src/omp.jl
index 65e04b5..b1ee1a9 100644
--- a/src/omp.jl
+++ b/src/omp.jl
@@ -31,7 +31,6 @@ Base.lastindex(t::OMPTracer) = lastindex(t.iteration)
     Ax = b.
     ```
 """
-
 function asp_omp(
     A::Union{AbstractMatrix, AbstractLinearOperator},
     b::Vector,
@@ -47,9 +46,8 @@ function asp_omp(
     optTol::Real = 1e-05,
     gapTol::Real = 1e-06,
     pivTol::Real = 1e-12,
-    actMax::Real = Inf)
+    actMax::Union{Real, Nothing} = nothing) 
     
-    # Start the clock and size up the problem.
     time0 = time()
 
     z = A' * b
@@ -59,16 +57,12 @@ function asp_omp(
     nprodA = 0
     nprodAt = 1
 
-    # Initialize the tracer
-
     tracer = OMPTracer(
-        Int[],                  # iteration
-        Float64[],              # lambda
-        Vector{SparseVector{Float64}}() # now stores full sparse solutions
+        Int[],                 
+        Float64[],              
+        Vector{SparseVector{Float64}}() 
     )
 
-    # Print log header.
-
     if loglevel > 0
         @info "-"^124
         @info @sprintf("%-30s : %-10d    %-30s : %-10.4e", "No. rows", m, "λ", λin)
@@ -77,6 +71,7 @@ function asp_omp(
         @info "-"^124
     end
 
+
     # Initialize local variables.
     EXIT_INFO = Dict(
         :EXIT_OPTIMAL => "Optimal solution found -- full Newton step",
@@ -85,8 +80,8 @@ function asp_omp(
         :EXIT_LAMBDA => "Reached minimum value of lambda",
         :EXIT_RHS_ZERO => "b = 0. The solution is x = 0",
         :EXIT_UNCONSTRAINED => "Unconstrained solution r = b is optimal",
-        :EXIT_UNKNOWN => "unknown exit",
-        :EXIT_ACTMAX => "Max no. of active constraints reached"
+        :EXIT_ACTMAX => "Max no. of active constraints reached",
+        :EXIT_UNKNOWN => "unknown exit"
     )
 
     itn = 0
@@ -95,7 +90,6 @@ function asp_omp(
     zerovec = zeros(Float64, n)
     p = 0
 
-    # Quick exit if the RHS is zero.
     if norm(b, Inf) == 0
         r = zeros(m)
         eFlag = :EXIT_RHS_ZERO
@@ -110,12 +104,21 @@ function asp_omp(
 
     if eFlag != :EXIT_UNKNOWN || active === nothing
         active = Vector{Int}([])
-        R = Matrix{Float64}(undef, 0, 0)
-        S = Matrix{Float64}(undef, m, 0)
+    end
+    if state === nothing
         state = zeros(Int, n)
     end
-    @info @sprintf("%4s  %8s %12s %12s %12s", "Itn", "Var", "λ", "rNorm", "xNorm")
+    if R === nothing
+        R = Matrix{Float64}(undef, 0, 0)
+    end
+
+    if actMax === nothing
+        actMax = size(A, 2)
+    end
 
+    if loglevel>0
+        @info @sprintf("%4s  %8s %12s %12s %12s", "Itn", "Var", "λ", "rNorm", "xNorm")
+    end
 
     # Main loop.
     while true
@@ -127,7 +130,7 @@ function asp_omp(
             nprodAt += 1
             zmax = norm(z, Inf)
         else
-            x, r = csne(R, S, b)
+            x,y = csne(R, S, vec(b))
             if norm(x, Inf) > 1e12
                 eFlag = :EXIT_SINGULAR_LS
                 break
@@ -140,7 +143,9 @@ function asp_omp(
         rNorm = norm(r, 2)
         xNorm = norm(x, 1)
 
-        @info @sprintf("%4i  %8i %12.5e %12.5e %12.5e", itn, p, zmax, rNorm, xNorm)
+        if loglevel>0
+            @info @sprintf("%4i  %8i %12.5e %12.5e %12.5e", itn, p, zmax, rNorm, xNorm)
+        end
 
         # Check exit conditions.
         if eFlag != :EXIT_UNKNOWN
@@ -151,10 +156,10 @@ function asp_omp(
             eFlag = :EXIT_OPTIMAL
         elseif itn >= itnMax
             eFlag = :EXIT_TOO_MANY_ITNS
-        elseif eFlag == :EXIT_UNKNOWN && length(active) >= actMax
+        elseif itn == actMax
             eFlag = :EXIT_ACTMAX
         end
-                
+
         if eFlag != :EXIT_UNKNOWN
             break
         end
@@ -181,22 +186,20 @@ function asp_omp(
         zerovec[p] = 0
 
         R = qraddcol(S, R, a)  # Update R
-        S = [S a]
+        S = hcat(S, a)  # Expand S, active
 
         push!(tracer.iteration, itn)
         push!(tracer.lambda, zmax)
-        sparse_x_full = spzeros(n)
-        sparse_x_full[copy(active)] = copy(x)  
+        sparse_x_full = SparseVector(n, copy(active), copy(x))
         push!(tracer.solution, copy(sparse_x_full))
         push!(active, p)
 
     end #while true
 
-    # push!(tracer.iteration, itn)
-    # push!(tracer.lambda, zmax)
-    # sparse_x_full = spzeros(n)
-    # sparse_x_full[copy(active)] = copy(x)  
-    # push!(tracer.solution, copy(sparse_x_full))
+    push!(tracer.iteration, itn)
+    push!(tracer.lambda, zmax)
+    sparse_x_full = SparseVector(n, copy(active), copy(x))
+    push!(tracer.solution, copy(sparse_x_full))
 
     tottime = time() - time0
     if loglevel > 0
@@ -207,4 +210,4 @@ function asp_omp(
         @info "\n"
     end
     return tracer
-end
+end
\ No newline at end of file
diff --git a/test/test_bpdn.jl b/test/test_bpdn.jl
index 106b597..257ba41 100644
--- a/test/test_bpdn.jl
+++ b/test/test_bpdn.jl
@@ -1,11 +1,11 @@
 
 # ------------------------------------------------------------------
-#  Test Basis pursuit
+#  Test BPDN & OMP
 # ------------------------------------------------------------------ 
 
 using Test, LinearAlgebra, Random, SparseArrays, ActiveSetPursuit, LinearOperators
 
-function test_bpdn()
+function test_pursuits()
     m = 600
     n = 2560
     k = 20
@@ -21,12 +21,16 @@ function test_bpdn()
     b = A * x
     
     # Solve the basis pursuit problem
-    tracer = asp_bpdn(A, b, 0.0, loglevel =0);
+    tracer_bpdn = asp_bpdn(A, b, 0.0, loglevel =0);
+    tracer_omp = asp_omp(A, b, 0.0, loglevel =0);
+    xx_bpdn, _ = tracer_bpdn[end]
+    xx_omp, _ = tracer_omp[end]
 
-    xx, λ = tracer[end]
-    pFeas = norm(A * xx - b, Inf) / max(1, norm(xx, Inf))
-    @test pFeas <= 1e-6
+    pFeas_bpdn = norm(A * xx_bpdn - b, Inf) / max(1, norm(xx_bpdn, Inf))
+    pFeas_omp = norm(A * xx_omp - b, Inf) / max(1, norm(xx_omp, Inf))
 
+    @test pFeas_bpdn <= 1e-6
+    @test pFeas_omp <= 1e-6
     # ------------------------
     # Linear operator
     # ------------------------
@@ -35,15 +39,20 @@ function test_bpdn()
     b_op = OP * x
         
     # Solve the basis pursuit problem
-    tracer_op = asp_bpdn(OP, b_op, 0.0, loglevel =0);
+    tracer_op_bpdn = asp_bpdn(OP, b_op, 0.0, loglevel =0);
+    tracer_op_omp = asp_omp(OP, b_op, 0.0, loglevel =0);
 
-    xx_op, λ_op = tracer_op[end]
-    pFeas_op = norm(OP * xx_op - b_op, Inf) / max(1, norm(xx_op, Inf))
-    @test pFeas_op <= 1e-6
+    xx_op_bpdn, _ = tracer_op_bpdn[end]
+    xx_op_omp, _ = tracer_op_omp[end]
+
+    pFeas_bpdn_op = norm(OP * xx_op_bpdn - b_op, Inf) / max(1, norm(xx_op_bpdn, Inf))
+    pFeas_omp_op = norm(OP * xx_op_omp - b_op, Inf) / max(1, norm(xx_op_omp, Inf))
+    @test pFeas_bpdn_op <= 1e-6
+    @test pFeas_omp_op <= 1e-6
 end
 
 Random.seed!(1234)
 
 for ntest = 1:10
-    test_bpdn()
+    test_pursuits()
 end
\ No newline at end of file

From ddc3f80a097a6ed6583d628dcdc645475e150f89 Mon Sep 17 00:00:00 2001
From: tinatorabi <tntorabii@gmail.com>
Date: Fri, 8 Nov 2024 21:25:12 -0800
Subject: [PATCH 5/5] in-place

---
 src/BPDual.jl | 40 ++++++++++++++++++++++++++--------------
 src/omp.jl    | 21 +++++++++++++++------
 2 files changed, 41 insertions(+), 20 deletions(-)

diff --git a/src/BPDual.jl b/src/BPDual.jl
index 212ef9e..ee39ffe 100644
--- a/src/BPDual.jl
+++ b/src/BPDual.jl
@@ -100,6 +100,13 @@ function bpdual(
     # ------------------------------------------------------------------
     m, n = size(A)
 
+
+    work = rand(n)
+    work2 = rand(n)
+    work3 = rand(n)
+    work4 = rand(n)
+    work5 = rand(m)
+
     tracer = ASPTracer(
         Int[],                  # iteration
         Float64[],              # lambda
@@ -110,8 +117,8 @@ function bpdual(
         active = Vector{Int}([])
         state = Vector{Int}(zeros(Int, n))
         y = Vector{Float64}(zeros(Float64, m))
-        S = Matrix{Float64}(zeros(Float64, m, 0))
-        R = Matrix{Float64}(zeros(Float64, 0, 0))
+        S = Matrix{Float64}(zeros(Float64, m, n))
+        R = Matrix{Float64}(zeros(Float64, n, n))
     end
 
     if homotopy
@@ -171,6 +178,8 @@ function bpdual(
     numtrim = 0
     nprodA = 0
     nprodAt = 0
+    cur_r_size = 0
+
 
     # ------------------------------------------------------------------
     # Cold/warm-start initialization.
@@ -211,25 +220,25 @@ function bpdual(
         sL, sU = infeasibilities(bl, bu, z)
         g = b - λ*y  # Steepest-descent direction
 
-        if isempty(R)
+        if isempty(R[1:cur_r_size,1:cur_r_size])
             condS = 1
         else
-            rmin = minimum(diag(R))
-            rmax = maximum(diag(R))
+            rmin = minimum(diag(R[1:cur_r_size, 1:cur_r_size]))
+            rmax = maximum(diag(R[1:cur_r_size, 1:cur_r_size]))
             condS = rmax / rmin
         end
 
         if condS > 1e+10
             eFlag = :EXIT_SINGULAR_LS
             # Pad x with enough zeros to make it compatible with S.
-            npad = size(S, 2) - size(x, 1)
+            npad = size(S[:, 1:cur_r_size], 2) - size(x, 1)
             x = [x; zeros(npad)]
         else
-            dx, dy = newtonstep(S, R, g, x, λ)
+            dx, dy = newtonstep(S[:,1:cur_r_size], R[1:cur_r_size, 1:cur_r_size], g, x, λ)
             x .+= dx
         end
 
-        r = b - S*x
+        r = b - S[:,1:cur_r_size]*x
 
         # Print to log.
         yNorm = norm(y, 2)
@@ -273,7 +282,7 @@ function bpdual(
                 @info "\nOptimal solution found. Trimming multipliers..."
             end
             g = b - λin*y
-            trimx(x, S, R, active, state, g, b, λ, feaTol, optTol, loglevel)
+            trimx(x, S[:, 1:cur_r_size], R[1:cur_r_size, 1:cur_r_size], active, state, g, b, λ, feaTol, optTol, loglevel)
             numtrim = nact - length(active)
             nact = length(active)
         end
@@ -288,7 +297,7 @@ function bpdual(
         p = q = 0
 
         if homotopy
-            x, dy, dz, step, λ, p = htpynewlam(active, state, A, R, S, x, y, sL, sU, λ, lamFinal)
+            x, dy, dz, step, λ, p = htpynewlam(active, state, A, R[1:cur_r_size, 1:cur_r_size], S[:,1:cur_r_size], x, y, sL, sU, λ, lamFinal)
             nprodAt += 1
         else
             if norm(dy, Inf) < eps()        
@@ -339,8 +348,9 @@ function bpdual(
             a = A * zerovec
             nprodA += 1
             zerovec[p] = 0
-            R = qraddcol(S, R, a)
-            S = [S a]
+            qraddcol!(S, R, a, cur_r_size, work, work2, work3, work4, work5)  # Update R
+            cur_r_size +=1             
+            # S = [S a]
             push!(active, p)
             push!(x, 0)
         else
@@ -358,10 +368,12 @@ function bpdual(
                 _, qa = findmax(abs.(x .* dropa))
                 q = active[qa]
                 state[q] = 0
-                S = S[:, 1:nact .!= qa]
+                # S = S[:, 1:nact .!= qa]
                 deleteat!(active, qa)
                 deleteat!(x, qa)
-                R = qrdelcol(R, qa)
+                # R = qrdelcol(R, qa)
+                qrdelcol!(S, R, qa)
+                cur_r_size -=1           
             else
                 eFlag = :EXIT_OPTIMAL
             end
diff --git a/src/omp.jl b/src/omp.jl
index b1ee1a9..d315c5d 100644
--- a/src/omp.jl
+++ b/src/omp.jl
@@ -54,6 +54,13 @@ function asp_omp(
     
     m = length(b)
     n = length(z)
+
+    work = rand(size(A,2))
+    work2 = rand(size(A,2))
+    work3 = rand(size(A,2))
+    work4 = rand(size(A,2))
+    work5 = rand(size(A,1))
+
     nprodA = 0
     nprodAt = 1
 
@@ -89,6 +96,7 @@ function asp_omp(
     x = zeros(Float64, 0)
     zerovec = zeros(Float64, n)
     p = 0
+    cur_r_size = 0
 
     if norm(b, Inf) == 0
         r = zeros(m)
@@ -109,7 +117,8 @@ function asp_omp(
         state = zeros(Int, n)
     end
     if R === nothing
-        R = Matrix{Float64}(undef, 0, 0)
+        R = Matrix{Float64}(undef, size(A,2), size(A,2))
+        S = Matrix{Float64}(undef, size(A,1), size(A,2))
     end
 
     if actMax === nothing
@@ -130,13 +139,13 @@ function asp_omp(
             nprodAt += 1
             zmax = norm(z, Inf)
         else
-            x,y = csne(R, S, vec(b))
+            x,y = csne(R[1:cur_r_size, 1:cur_r_size], S[:,1:cur_r_size], vec(b))
             if norm(x, Inf) > 1e12
                 eFlag = :EXIT_SINGULAR_LS
                 break
             end
 
-            Sx = S * x
+            Sx = S[:,1:cur_r_size] * x
             r = b - Sx
         end
 
@@ -185,9 +194,9 @@ function asp_omp(
         nprodA += 1
         zerovec[p] = 0
 
-        R = qraddcol(S, R, a)  # Update R
-        S = hcat(S, a)  # Expand S, active
-
+        qraddcol!(S, R, a, cur_r_size, work, work2, work3, work4, work5)  # Update R
+        # S = hcat(S, a)  # Expand S, active
+        cur_r_size +=1 
         push!(tracer.iteration, itn)
         push!(tracer.lambda, zmax)
         sparse_x_full = SparseVector(n, copy(active), copy(x))