From 002b3927fbd461e36c6ec234b928a6a3f4dab25a Mon Sep 17 00:00:00 2001
From: LJS42 <1457928939@qq.com>
Date: Wed, 26 Nov 2025 12:24:14 +0100
Subject: [PATCH 01/20] start

---
 src/polytope_blmos.jl | 5 +++++
 1 file changed, 5 insertions(+)

diff --git a/src/polytope_blmos.jl b/src/polytope_blmos.jl
index 1656ddb21..e533cf226 100644
--- a/src/polytope_blmos.jl
+++ b/src/polytope_blmos.jl
@@ -716,3 +716,8 @@ function rounding_hyperplane_heuristic(
     end
     return [z], false
 end
+
+struct ReverseKnapsackLMO <: FrankWolfe.LinearMinimizationOracle
+    N::Float64
+    upper_bounds::Vector{Float64}
+end
\ No newline at end of file

From fcdf67eb7d65d162b9bfb80e9bea18ddd37daeb4 Mon Sep 17 00:00:00 2001
From: LJS42 <1457928939@qq.com>
Date: Wed, 26 Nov 2025 17:50:41 +0100
Subject: [PATCH 02/20] testing

---
 src/polytope_blmos.jl | 123 ++++++++++++++++++++++++++++++++++++++++--
 test/LMO_test.jl      | 103 +++++++++++++++++++++++++++++++++++
 2 files changed, 223 insertions(+), 3 deletions(-)

diff --git a/src/polytope_blmos.jl b/src/polytope_blmos.jl
index e533cf226..1ae32e51f 100644
--- a/src/polytope_blmos.jl
+++ b/src/polytope_blmos.jl
@@ -717,7 +717,124 @@ function rounding_hyperplane_heuristic(
     return [z], false
 end
 
-struct ReverseKnapsackLMO <: FrankWolfe.LinearMinimizationOracle
-    N::Float64
-    upper_bounds::Vector{Float64}
+"""
+    2normBallBLMO()
+
+BLMO denotes the L2normBall, It is unit ball which means R = 1
+"""
+
+struct L2normBallBLMO <: FrankWolfe.LinearMinimizationOracle
+end
+
+function bounded_compute_extreme_point(blmo::L2normBallBLMO, d, lb, ub, int_vars; kwargs...)
+    max = abs(d[1])
+    max_idx = 1
+    idx = 0
+    for idx in int_vars
+        if abs(d[idx]) > max
+            max = -d[idx]
+            max_idx = idx
+        end
+        i = findfirst(x -> x == idx, int_vars)
+        if lb[i] > 0 
+            v = zeros(length(d))
+            v[idx] = 1
+            return v
+        elseif ub[i] < 0 
+            v = zeros(length(d))
+            v[idx] = -1
+            return v
+        end
+    end
+    v = -d
+    v[int_vars] .= 0
+    v = v ./ norm(v)
+    prod = dot(v,d)
+    if max < prod 
+        v = zeros(length(d))
+        v[max_idx] = 1
+    end
+    return v
+end
+
+function is_simple_linear_feasible(lmo::L2normBallBLMO, v)
+    if norm(v) > 1 + 1e-6
+        @debug "norm(v) : $(norm(v)) > 1"
+        return false
+    end
+    return true
+end
+
+function check_feasibility(lmo::L2normBallBLMO, lb, ub, int_vars, n)
+    if any(lb .> 1) || any(ub .< -1)
+        return INFEASIBLE
+    end
+    count = 0
+    for idx in 1:length(int_vars)
+        if !(lb[idx] <= 0 <= ub[idx])
+            count = count + 1
+        end
+        if count > 1 
+            return INFEASIBLE
+        end
+    end
+    return OPTIMAL
+end 
+
+"""
+    ConvexHullLMO(Vector{Vector{Float64}})
+
+Linear Minimization Oracle over the convex hull of a finite set of integer points.
+The extreme points are the integer vertices provided at construction.
+Primarily used for integer programming relaxations.
+"""
+
+struct ConvexHullLMO{T, V <: AbstractVector{T}} <: FrankWolfe.LinearMinimizationOracle
+    vertices::Vector{V}
+end
+
+function bounded_compute_extreme_point(
+    lmo::ConvexHullLMO,
+    direction::AbstractVector{T},
+    lb,
+    ub, 
+    int_vars;
+    kwargs...,
+) where {T}
+    # Find the integer vertex that minimizes the linear objective
+    min_val = Inf
+    best_vertex = lmo.vertices[1]
+    
+    @inbounds for vertex in lmo.vertices
+        val = dot(direction, vertex)
+        if val < min_val
+            min_val = val
+            best_vertex = vertex
+        end
+    end
+    
+    # Return as Float64 for compatibility with continuous optimization
+    return convert(Vector{Float64}, best_vertex)
+end
+
+function is_simple_linear_feasible(lmo::ConvexHullLMO, v)
+    return true
+end
+
+
+function check_feasibility(lmo::ConvexHullLMO, lb, ub, int_vars, n)
+    for vertex in lmo.vertices
+        is_feasible_vertex = true
+        
+        for i in 1:n 
+            if vertex[i] < lb[i] - eps(Float64) || vertex[i] > ub[i] + eps(Float64)
+                is_feasible_vertex = false
+                break
+            end
+        end
+        if is_feasible_vertex
+            return OPTIMAL
+        end
+    end
+    return INFEASIBLE
 end
\ No newline at end of file
diff --git a/test/LMO_test.jl b/test/LMO_test.jl
index f44f4415c..d0bd1b29b 100644
--- a/test/LMO_test.jl
+++ b/test/LMO_test.jl
@@ -8,6 +8,7 @@ import Bonobo
 using HiGHS
 using Printf
 using Dates
+using LinearAlgebra
 const MOI = MathOptInterface
 const MOIU = MOI.Utilities
 using StableRNGs
@@ -217,3 +218,105 @@ diffi = x_sol + 0.3 * rand([-1, 1], n)
     @test sum(isapprox.(x, x_sol, atol=1e-6, rtol=1e-2)) == n
     @test isapprox(f(x), f(result[:raw_solution]), atol=1e-6, rtol=1e-3)
 end
+
+# Test for L2normBallBLMO
+@testset "L2normBall BLMO" begin
+    n = 20
+    
+    # Generate a solution inside the L2 ball (||x|| <= 1)
+    x_sol = randn(rng, n)
+    x_sol = x_sol ./ (norm(x_sol) * 1.5)  
+    
+    # Round some coordinates to integers for testing
+    num_int = 5
+    int_indices = sort(rand(rng, 1:n, num_int))
+    for idx in int_indices
+        x_sol[idx] = 0
+    end
+    
+    function f(x)
+        return 0.5 * sum((x[i] - x_sol[i])^2 for i in eachindex(x))
+    end
+    
+    function grad!(storage, x)
+        @. storage = x - x_sol
+    end
+    
+    # Create L2normBallBLMO
+    blmo = Boscia.L2normBallBLMO()
+    
+    # Bounds: each variable in [-1, 1] due to L2 ball constraint
+    lower_bounds = fill(-2.0, num_int)
+    upper_bounds = fill(2.0, num_int)
+    
+    # Some variables are integer
+    int_vars = collect(int_indices)
+    
+    x, _, result = Boscia.solve(
+        f, 
+        grad!, 
+        blmo, 
+        lower_bounds, 
+        upper_bounds, 
+        int_vars, 
+        n
+    )
+    
+    # Check objective value
+    @test sum(isapprox.(x, x_sol, atol=1e-6, rtol=1e-2)) == n  # Should improve or stay similar
+    @test isapprox(f(x), f(result[:raw_solution]), atol=1e-6, rtol=1e-3)
+end
+
+
+# Test for ConvexHullBLMO
+@testset "ConvexHull BLMO" begin
+    n = 20
+    
+    # Generate random integer vertices for the convex hull
+    num_vertices = 30
+    vertices = [rand(rng, 0:10, n) for _ in 1:num_vertices]
+    
+    # Create a solution that is a convex combination of vertices
+    # Pick a few vertices randomly
+    num_active = 3
+    active_indices = rand(rng, 1:num_vertices, num_active)
+    weights = rand(rng, num_active)
+    weights = weights ./ sum(weights)  # normalize to sum to 1
+    
+    x_sol = zeros(n)
+    for (idx, w) in zip(active_indices, weights)
+        x_sol .+= w .* vertices[idx]
+    end
+    
+    # Add small perturbation
+    diffi = x_sol + 0.3 * rand(rng, n) .* randn(rng, n)
+    
+    function f(x)
+        return 0.5 * sum((x[i] - diffi[i])^2 for i in eachindex(x))
+    end
+    
+    function grad!(storage, x)
+        @. storage = x - diffi
+    end
+    
+    # Create ConvexHullBLMO with the vertices
+    blmo = Boscia.ConvexHullLMO(vertices)
+    
+    # All variables are integer in the original problem
+    int_vars = collect(1:n)
+    
+    x, _, result = Boscia.solve(
+        f, 
+        grad!, 
+        blmo, 
+        fill(0.0, n),
+        fill(1.0 * N, n),
+        int_vars, 
+        n
+    )
+    
+    # Check that solution is close to target
+    @test sum(isapprox.(x, x_sol, atol=1e-6, rtol=1e-2)) == n
+    @test isapprox(f(x), f(result[:raw_solution]), atol=1e-6, rtol=1e-3)
+end
+

From 9321f79e87043cd1a9bc07a747c2237ea04d1e4f Mon Sep 17 00:00:00 2001
From: LJS42 <1457928939@qq.com>
Date: Sat, 29 Nov 2025 22:16:24 +0100
Subject: [PATCH 03/20] pass all the test

---
 src/polytope_blmos.jl | 171 ++++++++++++++++++++++++++++++++++++------
 test/LMO_test.jl      |  12 +--
 2 files changed, 151 insertions(+), 32 deletions(-)

diff --git a/src/polytope_blmos.jl b/src/polytope_blmos.jl
index 1ae32e51f..2d2ec1485 100644
--- a/src/polytope_blmos.jl
+++ b/src/polytope_blmos.jl
@@ -793,48 +793,173 @@ struct ConvexHullLMO{T, V <: AbstractVector{T}} <: FrankWolfe.LinearMinimization
     vertices::Vector{V}
 end
 
+function proj_simplex!(λ::AbstractVector{T}) where T<:Real
+    n = length(λ)
+    u = sort(collect(λ); rev=true)
+    cssv = cumsum(u)
+    rho = findlast(i -> u[i] + (1 - cssv[i]) / i > 0, 1:n)
+    if rho === nothing
+        fill!(λ, 1/n)
+        return λ
+    end
+    theta = (cssv[rho] - 1) / rho
+    @inbounds for i in eachindex(λ)
+        λ[i] = max(λ[i] - theta, zero(T))
+    end
+    return λ
+end
+
+function vertices_matrix(vertices::Vector{<:AbstractVector{T}}) where T
+    m = length(vertices)
+    n = length(vertices[1])
+    V = zeros(Float64, n, m)
+    @inbounds for j in 1:m
+        V[:, j] .= Float64.(vertices[j])
+    end
+    return V
+end
+
+function try_convex_combination_in_box(vertices::Vector{<:AbstractVector}, lb::AbstractVector, ub::AbstractVector; tol=1e-8, max_iters=200)
+    m = length(vertices)
+    if m == 0
+        return nothing
+    end
+    V = vertices_matrix(vertices)  # n x m
+    # initial λ: uniform
+    λ = fill(1.0 / m, m)
+    # target mid point
+    mid = (Float64.(lb) .+ Float64.(ub)) ./ 2.0
+    # gradient descent w/ projection onto simplex
+    α = 1.0
+    for iter in 1:max_iters
+        x = V * λ
+        grad = 2.0 * (V' * (x - mid))
+        # step
+        λ -= α * grad
+        proj_simplex!(λ)
+        # check residual
+        x = V * λ
+        if all(x .>= Float64.(lb) .- 1e-8) && all(x .<= Float64.(ub) .+ 1e-8)
+            return x  # feasible convex combination found
+        end
+        # optionally reduce step size
+        α *= 0.99
+    end
+    return nothing
+end
+
 function bounded_compute_extreme_point(
     lmo::ConvexHullLMO,
     direction::AbstractVector{T},
     lb,
-    ub, 
+    ub,
     int_vars;
-    kwargs...,
+    kwargs...
 ) where {T}
-    # Find the integer vertex that minimizes the linear objective
-    min_val = Inf
-    best_vertex = lmo.vertices[1]
-    
-    @inbounds for vertex in lmo.vertices
-        val = dot(direction, vertex)
-        if val < min_val
-            min_val = val
-            best_vertex = vertex
+
+    verts = lmo.vertices
+    nverts = length(verts)
+    n = length(direction)
+    feasible_vs = Vector{Vector{Float64}}()
+    @inbounds for v in verts
+        vfloat = Float64.(v)
+        ok = true
+        for i in 1:n
+            if vfloat[i] < Float64(lb[i]) - eps(Float64) || vfloat[i] > Float64(ub[i]) + eps(Float64)
+                ok = false
+                break
+            end
+        end
+        if ok
+            push!(feasible_vs, vfloat)
         end
     end
-    
-    # Return as Float64 for compatibility with continuous optimization
-    return convert(Vector{Float64}, best_vertex)
+
+    if !isempty(feasible_vs)
+        best_val = Inf
+        best_v = feasible_vs[1]
+        @inbounds for v in feasible_vs
+            val = dot(Float64.(direction), v)
+            if val < best_val
+                best_val = val
+                best_v = v
+            end
+        end
+        return best_v
+    end
+
+    x_comb = try_convex_combination_in_box(verts, lb, ub)
+    if x_comb !== nothing
+        return x_comb
+    end
+
+    best_val = Inf
+    best_v = Float64.(verts[1])
+    dirf = Float64.(direction)
+    @inbounds for v in verts
+        vfloat = Float64.(v)
+        val = dot(dirf, vfloat)
+        if val < best_val
+            best_val = val
+            best_v = vfloat
+        end
+    end
+    # clip
+    clipped = clamp.(best_v, Float64.(lb), Float64.(ub))
+    @warn "ConvexHullLMO: no single vertex or convex combination found inside [lb,ub]. Returning clipped best-vertex. This may cause degenerate branching if many vertices share same int value." clipped=clipped
+    return clipped
 end
 
+# -----------------------
+# is_simple_linear_feasible: quick check if any vertex individually satisfies bounds.
 function is_simple_linear_feasible(lmo::ConvexHullLMO, v)
-    return true
+    # v is candidate point; simply check if in box
+    return all(Float64.(v) .>= Float64.(v) .- 0.0)  # trivial placeholder: always true in user code previously
 end
 
-
 function check_feasibility(lmo::ConvexHullLMO, lb, ub, int_vars, n)
-    for vertex in lmo.vertices
-        is_feasible_vertex = true
-        
-        for i in 1:n 
-            if vertex[i] < lb[i] - eps(Float64) || vertex[i] > ub[i] + eps(Float64)
-                is_feasible_vertex = false
+    verts = lmo.vertices
+    if isempty(verts)
+        return INFEASIBLE
+    end
+    @inbounds for v in verts
+        vfloat = Float64.(v)
+        ok = true
+        for i in 1:n
+            if vfloat[i] < Float64(lb[i]) - eps(Float64) || vfloat[i] > Float64(ub[i]) + eps(Float64)
+                ok = false
                 break
             end
         end
-        if is_feasible_vertex
+        if ok
             return OPTIMAL
         end
     end
+
+    nverts = length(verts)
+    mins = fill(Inf, n)
+    maxs = fill(-Inf, n)
+    @inbounds for v in verts
+        vfloat = Float64.(v)
+        for i in 1:n
+            if vfloat[i] < mins[i]
+                mins[i] = vfloat[i]
+            end
+            if vfloat[i] > maxs[i]
+                maxs[i] = vfloat[i]
+            end
+        end
+    end
+    for i in 1:n
+        if maxs[i] < Float64(lb[i]) - 1e-12 || mins[i] > Float64(ub[i]) + 1e-12
+            return INFEASIBLE
+        end
+    end
+
+    x = try_convex_combination_in_box(verts, lb, ub)
+    if x !== nothing
+        return OPTIMAL
+    end
+
     return INFEASIBLE
 end
\ No newline at end of file
diff --git a/test/LMO_test.jl b/test/LMO_test.jl
index d0bd1b29b..490d15bef 100644
--- a/test/LMO_test.jl
+++ b/test/LMO_test.jl
@@ -279,14 +279,8 @@ end
     # Create a solution that is a convex combination of vertices
     # Pick a few vertices randomly
     num_active = 3
-    active_indices = rand(rng, 1:num_vertices, num_active)
-    weights = rand(rng, num_active)
-    weights = weights ./ sum(weights)  # normalize to sum to 1
-    
-    x_sol = zeros(n)
-    for (idx, w) in zip(active_indices, weights)
-        x_sol .+= w .* vertices[idx]
-    end
+
+    x_sol = vertices[3]
     
     # Add small perturbation
     diffi = x_sol + 0.3 * rand(rng, n) .* randn(rng, n)
@@ -314,7 +308,7 @@ end
         int_vars, 
         n
     )
-    
+
     # Check that solution is close to target
     @test sum(isapprox.(x, x_sol, atol=1e-6, rtol=1e-2)) == n
     @test isapprox(f(x), f(result[:raw_solution]), atol=1e-6, rtol=1e-3)

From c467d1633b4d5d7036ed64516b6fe7d5d144e6ff Mon Sep 17 00:00:00 2001
From: LJS42 <1457928939@qq.com>
Date: Sat, 29 Nov 2025 22:30:15 +0100
Subject: [PATCH 04/20] format check

---
 src/polytope_blmos.jl | 44 +++++++++++++++++++--------------
 test/LMO_test.jl      | 57 +++++++++++++++----------------------------
 2 files changed, 46 insertions(+), 55 deletions(-)

diff --git a/src/polytope_blmos.jl b/src/polytope_blmos.jl
index 2d2ec1485..4642e7ba6 100644
--- a/src/polytope_blmos.jl
+++ b/src/polytope_blmos.jl
@@ -723,8 +723,7 @@ end
 BLMO denotes the L2normBall, It is unit ball which means R = 1
 """
 
-struct L2normBallBLMO <: FrankWolfe.LinearMinimizationOracle
-end
+struct L2normBallBLMO <: FrankWolfe.LinearMinimizationOracle end
 
 function bounded_compute_extreme_point(blmo::L2normBallBLMO, d, lb, ub, int_vars; kwargs...)
     max = abs(d[1])
@@ -736,11 +735,11 @@ function bounded_compute_extreme_point(blmo::L2normBallBLMO, d, lb, ub, int_vars
             max_idx = idx
         end
         i = findfirst(x -> x == idx, int_vars)
-        if lb[i] > 0 
+        if lb[i] > 0
             v = zeros(length(d))
             v[idx] = 1
             return v
-        elseif ub[i] < 0 
+        elseif ub[i] < 0
             v = zeros(length(d))
             v[idx] = -1
             return v
@@ -749,8 +748,8 @@ function bounded_compute_extreme_point(blmo::L2normBallBLMO, d, lb, ub, int_vars
     v = -d
     v[int_vars] .= 0
     v = v ./ norm(v)
-    prod = dot(v,d)
-    if max < prod 
+    prod = dot(v, d)
+    if max < prod
         v = zeros(length(d))
         v[max_idx] = 1
     end
@@ -774,12 +773,12 @@ function check_feasibility(lmo::L2normBallBLMO, lb, ub, int_vars, n)
         if !(lb[idx] <= 0 <= ub[idx])
             count = count + 1
         end
-        if count > 1 
+        if count > 1
             return INFEASIBLE
         end
     end
     return OPTIMAL
-end 
+end
 
 """
     ConvexHullLMO(Vector{Vector{Float64}})
@@ -789,17 +788,17 @@ The extreme points are the integer vertices provided at construction.
 Primarily used for integer programming relaxations.
 """
 
-struct ConvexHullLMO{T, V <: AbstractVector{T}} <: FrankWolfe.LinearMinimizationOracle
+struct ConvexHullLMO{T,V<:AbstractVector{T}} <: FrankWolfe.LinearMinimizationOracle
     vertices::Vector{V}
 end
 
-function proj_simplex!(λ::AbstractVector{T}) where T<:Real
+function proj_simplex!(λ::AbstractVector{T}) where {T<:Real}
     n = length(λ)
     u = sort(collect(λ); rev=true)
     cssv = cumsum(u)
     rho = findlast(i -> u[i] + (1 - cssv[i]) / i > 0, 1:n)
     if rho === nothing
-        fill!(λ, 1/n)
+        fill!(λ, 1 / n)
         return λ
     end
     theta = (cssv[rho] - 1) / rho
@@ -809,7 +808,7 @@ function proj_simplex!(λ::AbstractVector{T}) where T<:Real
     return λ
 end
 
-function vertices_matrix(vertices::Vector{<:AbstractVector{T}}) where T
+function vertices_matrix(vertices::Vector{<:AbstractVector{T}}) where {T}
     m = length(vertices)
     n = length(vertices[1])
     V = zeros(Float64, n, m)
@@ -819,7 +818,13 @@ function vertices_matrix(vertices::Vector{<:AbstractVector{T}}) where T
     return V
 end
 
-function try_convex_combination_in_box(vertices::Vector{<:AbstractVector}, lb::AbstractVector, ub::AbstractVector; tol=1e-8, max_iters=200)
+function try_convex_combination_in_box(
+    vertices::Vector{<:AbstractVector},
+    lb::AbstractVector,
+    ub::AbstractVector;
+    tol=1e-8,
+    max_iters=200,
+)
     m = length(vertices)
     if m == 0
         return nothing
@@ -854,7 +859,7 @@ function bounded_compute_extreme_point(
     lb,
     ub,
     int_vars;
-    kwargs...
+    kwargs...,
 ) where {T}
 
     verts = lmo.vertices
@@ -865,7 +870,8 @@ function bounded_compute_extreme_point(
         vfloat = Float64.(v)
         ok = true
         for i in 1:n
-            if vfloat[i] < Float64(lb[i]) - eps(Float64) || vfloat[i] > Float64(ub[i]) + eps(Float64)
+            if vfloat[i] < Float64(lb[i]) - eps(Float64) ||
+               vfloat[i] > Float64(ub[i]) + eps(Float64)
                 ok = false
                 break
             end
@@ -906,7 +912,8 @@ function bounded_compute_extreme_point(
     end
     # clip
     clipped = clamp.(best_v, Float64.(lb), Float64.(ub))
-    @warn "ConvexHullLMO: no single vertex or convex combination found inside [lb,ub]. Returning clipped best-vertex. This may cause degenerate branching if many vertices share same int value." clipped=clipped
+    @warn "ConvexHullLMO: no single vertex or convex combination found inside [lb,ub]. Returning clipped best-vertex. This may cause degenerate branching if many vertices share same int value." clipped =
+        clipped
     return clipped
 end
 
@@ -926,7 +933,8 @@ function check_feasibility(lmo::ConvexHullLMO, lb, ub, int_vars, n)
         vfloat = Float64.(v)
         ok = true
         for i in 1:n
-            if vfloat[i] < Float64(lb[i]) - eps(Float64) || vfloat[i] > Float64(ub[i]) + eps(Float64)
+            if vfloat[i] < Float64(lb[i]) - eps(Float64) ||
+               vfloat[i] > Float64(ub[i]) + eps(Float64)
                 ok = false
                 break
             end
@@ -962,4 +970,4 @@ function check_feasibility(lmo::ConvexHullLMO, lb, ub, int_vars, n)
     end
 
     return INFEASIBLE
-end
\ No newline at end of file
+end
diff --git a/test/LMO_test.jl b/test/LMO_test.jl
index 490d15bef..3fa4996cf 100644
--- a/test/LMO_test.jl
+++ b/test/LMO_test.jl
@@ -222,46 +222,38 @@ end
 # Test for L2normBallBLMO
 @testset "L2normBall BLMO" begin
     n = 20
-    
+
     # Generate a solution inside the L2 ball (||x|| <= 1)
     x_sol = randn(rng, n)
-    x_sol = x_sol ./ (norm(x_sol) * 1.5)  
-    
+    x_sol = x_sol ./ (norm(x_sol) * 1.5)
+
     # Round some coordinates to integers for testing
     num_int = 5
     int_indices = sort(rand(rng, 1:n, num_int))
     for idx in int_indices
         x_sol[idx] = 0
     end
-    
+
     function f(x)
         return 0.5 * sum((x[i] - x_sol[i])^2 for i in eachindex(x))
     end
-    
+
     function grad!(storage, x)
         @. storage = x - x_sol
     end
-    
+
     # Create L2normBallBLMO
     blmo = Boscia.L2normBallBLMO()
-    
+
     # Bounds: each variable in [-1, 1] due to L2 ball constraint
     lower_bounds = fill(-2.0, num_int)
     upper_bounds = fill(2.0, num_int)
-    
+
     # Some variables are integer
     int_vars = collect(int_indices)
-    
-    x, _, result = Boscia.solve(
-        f, 
-        grad!, 
-        blmo, 
-        lower_bounds, 
-        upper_bounds, 
-        int_vars, 
-        n
-    )
-    
+
+    x, _, result = Boscia.solve(f, grad!, blmo, lower_bounds, upper_bounds, int_vars, n)
+
     # Check objective value
     @test sum(isapprox.(x, x_sol, atol=1e-6, rtol=1e-2)) == n  # Should improve or stay similar
     @test isapprox(f(x), f(result[:raw_solution]), atol=1e-6, rtol=1e-3)
@@ -271,46 +263,37 @@ end
 # Test for ConvexHullBLMO
 @testset "ConvexHull BLMO" begin
     n = 20
-    
+
     # Generate random integer vertices for the convex hull
     num_vertices = 30
     vertices = [rand(rng, 0:10, n) for _ in 1:num_vertices]
-    
+
     # Create a solution that is a convex combination of vertices
     # Pick a few vertices randomly
     num_active = 3
 
     x_sol = vertices[3]
-    
+
     # Add small perturbation
     diffi = x_sol + 0.3 * rand(rng, n) .* randn(rng, n)
-    
+
     function f(x)
         return 0.5 * sum((x[i] - diffi[i])^2 for i in eachindex(x))
     end
-    
+
     function grad!(storage, x)
         @. storage = x - diffi
     end
-    
+
     # Create ConvexHullBLMO with the vertices
     blmo = Boscia.ConvexHullLMO(vertices)
-    
+
     # All variables are integer in the original problem
     int_vars = collect(1:n)
-    
-    x, _, result = Boscia.solve(
-        f, 
-        grad!, 
-        blmo, 
-        fill(0.0, n),
-        fill(1.0 * N, n),
-        int_vars, 
-        n
-    )
+
+    x, _, result = Boscia.solve(f, grad!, blmo, fill(0.0, n), fill(1.0 * N, n), int_vars, n)
 
     # Check that solution is close to target
     @test sum(isapprox.(x, x_sol, atol=1e-6, rtol=1e-2)) == n
     @test isapprox(f(x), f(result[:raw_solution]), atol=1e-6, rtol=1e-3)
 end
-

From 12ad3ce4f001726556e881968efbb29c97a9040b Mon Sep 17 00:00:00 2001
From: LJS42 <1457928939@qq.com>
Date: Wed, 3 Dec 2025 09:53:05 +0100
Subject: [PATCH 05/20] save

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

diff --git a/test/LMO_test.jl b/test/LMO_test.jl
index 3fa4996cf..bb108d0db 100644
--- a/test/LMO_test.jl
+++ b/test/LMO_test.jl
@@ -291,7 +291,7 @@ end
     # All variables are integer in the original problem
     int_vars = collect(1:n)
 
-    x, _, result = Boscia.solve(f, grad!, blmo, fill(0.0, n), fill(1.0 * N, n), int_vars, n)
+    x, _, result = Boscia.solve(f, grad!, blmo, fill(0.0, n), fill(10.0, n), int_vars, n)
 
     # Check that solution is close to target
     @test sum(isapprox.(x, x_sol, atol=1e-6, rtol=1e-2)) == n

From 23cb9a2b357b7543787c504d5ab7a8817dcfdf14 Mon Sep 17 00:00:00 2001
From: LJS42 <1457928939@qq.com>
Date: Wed, 3 Dec 2025 14:29:25 +0100
Subject: [PATCH 06/20] seperate

---
 src/polytope_blmos.jl | 194 +-----------------------------------------
 test/LMO_test.jl      |  39 ---------
 2 files changed, 1 insertion(+), 232 deletions(-)

diff --git a/src/polytope_blmos.jl b/src/polytope_blmos.jl
index 4642e7ba6..a0ba57a6a 100644
--- a/src/polytope_blmos.jl
+++ b/src/polytope_blmos.jl
@@ -778,196 +778,4 @@ function check_feasibility(lmo::L2normBallBLMO, lb, ub, int_vars, n)
         end
     end
     return OPTIMAL
-end
-
-"""
-    ConvexHullLMO(Vector{Vector{Float64}})
-
-Linear Minimization Oracle over the convex hull of a finite set of integer points.
-The extreme points are the integer vertices provided at construction.
-Primarily used for integer programming relaxations.
-"""
-
-struct ConvexHullLMO{T,V<:AbstractVector{T}} <: FrankWolfe.LinearMinimizationOracle
-    vertices::Vector{V}
-end
-
-function proj_simplex!(λ::AbstractVector{T}) where {T<:Real}
-    n = length(λ)
-    u = sort(collect(λ); rev=true)
-    cssv = cumsum(u)
-    rho = findlast(i -> u[i] + (1 - cssv[i]) / i > 0, 1:n)
-    if rho === nothing
-        fill!(λ, 1 / n)
-        return λ
-    end
-    theta = (cssv[rho] - 1) / rho
-    @inbounds for i in eachindex(λ)
-        λ[i] = max(λ[i] - theta, zero(T))
-    end
-    return λ
-end
-
-function vertices_matrix(vertices::Vector{<:AbstractVector{T}}) where {T}
-    m = length(vertices)
-    n = length(vertices[1])
-    V = zeros(Float64, n, m)
-    @inbounds for j in 1:m
-        V[:, j] .= Float64.(vertices[j])
-    end
-    return V
-end
-
-function try_convex_combination_in_box(
-    vertices::Vector{<:AbstractVector},
-    lb::AbstractVector,
-    ub::AbstractVector;
-    tol=1e-8,
-    max_iters=200,
-)
-    m = length(vertices)
-    if m == 0
-        return nothing
-    end
-    V = vertices_matrix(vertices)  # n x m
-    # initial λ: uniform
-    λ = fill(1.0 / m, m)
-    # target mid point
-    mid = (Float64.(lb) .+ Float64.(ub)) ./ 2.0
-    # gradient descent w/ projection onto simplex
-    α = 1.0
-    for iter in 1:max_iters
-        x = V * λ
-        grad = 2.0 * (V' * (x - mid))
-        # step
-        λ -= α * grad
-        proj_simplex!(λ)
-        # check residual
-        x = V * λ
-        if all(x .>= Float64.(lb) .- 1e-8) && all(x .<= Float64.(ub) .+ 1e-8)
-            return x  # feasible convex combination found
-        end
-        # optionally reduce step size
-        α *= 0.99
-    end
-    return nothing
-end
-
-function bounded_compute_extreme_point(
-    lmo::ConvexHullLMO,
-    direction::AbstractVector{T},
-    lb,
-    ub,
-    int_vars;
-    kwargs...,
-) where {T}
-
-    verts = lmo.vertices
-    nverts = length(verts)
-    n = length(direction)
-    feasible_vs = Vector{Vector{Float64}}()
-    @inbounds for v in verts
-        vfloat = Float64.(v)
-        ok = true
-        for i in 1:n
-            if vfloat[i] < Float64(lb[i]) - eps(Float64) ||
-               vfloat[i] > Float64(ub[i]) + eps(Float64)
-                ok = false
-                break
-            end
-        end
-        if ok
-            push!(feasible_vs, vfloat)
-        end
-    end
-
-    if !isempty(feasible_vs)
-        best_val = Inf
-        best_v = feasible_vs[1]
-        @inbounds for v in feasible_vs
-            val = dot(Float64.(direction), v)
-            if val < best_val
-                best_val = val
-                best_v = v
-            end
-        end
-        return best_v
-    end
-
-    x_comb = try_convex_combination_in_box(verts, lb, ub)
-    if x_comb !== nothing
-        return x_comb
-    end
-
-    best_val = Inf
-    best_v = Float64.(verts[1])
-    dirf = Float64.(direction)
-    @inbounds for v in verts
-        vfloat = Float64.(v)
-        val = dot(dirf, vfloat)
-        if val < best_val
-            best_val = val
-            best_v = vfloat
-        end
-    end
-    # clip
-    clipped = clamp.(best_v, Float64.(lb), Float64.(ub))
-    @warn "ConvexHullLMO: no single vertex or convex combination found inside [lb,ub]. Returning clipped best-vertex. This may cause degenerate branching if many vertices share same int value." clipped =
-        clipped
-    return clipped
-end
-
-# -----------------------
-# is_simple_linear_feasible: quick check if any vertex individually satisfies bounds.
-function is_simple_linear_feasible(lmo::ConvexHullLMO, v)
-    # v is candidate point; simply check if in box
-    return all(Float64.(v) .>= Float64.(v) .- 0.0)  # trivial placeholder: always true in user code previously
-end
-
-function check_feasibility(lmo::ConvexHullLMO, lb, ub, int_vars, n)
-    verts = lmo.vertices
-    if isempty(verts)
-        return INFEASIBLE
-    end
-    @inbounds for v in verts
-        vfloat = Float64.(v)
-        ok = true
-        for i in 1:n
-            if vfloat[i] < Float64(lb[i]) - eps(Float64) ||
-               vfloat[i] > Float64(ub[i]) + eps(Float64)
-                ok = false
-                break
-            end
-        end
-        if ok
-            return OPTIMAL
-        end
-    end
-
-    nverts = length(verts)
-    mins = fill(Inf, n)
-    maxs = fill(-Inf, n)
-    @inbounds for v in verts
-        vfloat = Float64.(v)
-        for i in 1:n
-            if vfloat[i] < mins[i]
-                mins[i] = vfloat[i]
-            end
-            if vfloat[i] > maxs[i]
-                maxs[i] = vfloat[i]
-            end
-        end
-    end
-    for i in 1:n
-        if maxs[i] < Float64(lb[i]) - 1e-12 || mins[i] > Float64(ub[i]) + 1e-12
-            return INFEASIBLE
-        end
-    end
-
-    x = try_convex_combination_in_box(verts, lb, ub)
-    if x !== nothing
-        return OPTIMAL
-    end
-
-    return INFEASIBLE
-end
+end
\ No newline at end of file
diff --git a/test/LMO_test.jl b/test/LMO_test.jl
index bb108d0db..d7da828d2 100644
--- a/test/LMO_test.jl
+++ b/test/LMO_test.jl
@@ -258,42 +258,3 @@ end
     @test sum(isapprox.(x, x_sol, atol=1e-6, rtol=1e-2)) == n  # Should improve or stay similar
     @test isapprox(f(x), f(result[:raw_solution]), atol=1e-6, rtol=1e-3)
 end
-
-
-# Test for ConvexHullBLMO
-@testset "ConvexHull BLMO" begin
-    n = 20
-
-    # Generate random integer vertices for the convex hull
-    num_vertices = 30
-    vertices = [rand(rng, 0:10, n) for _ in 1:num_vertices]
-
-    # Create a solution that is a convex combination of vertices
-    # Pick a few vertices randomly
-    num_active = 3
-
-    x_sol = vertices[3]
-
-    # Add small perturbation
-    diffi = x_sol + 0.3 * rand(rng, n) .* randn(rng, n)
-
-    function f(x)
-        return 0.5 * sum((x[i] - diffi[i])^2 for i in eachindex(x))
-    end
-
-    function grad!(storage, x)
-        @. storage = x - diffi
-    end
-
-    # Create ConvexHullBLMO with the vertices
-    blmo = Boscia.ConvexHullLMO(vertices)
-
-    # All variables are integer in the original problem
-    int_vars = collect(1:n)
-
-    x, _, result = Boscia.solve(f, grad!, blmo, fill(0.0, n), fill(10.0, n), int_vars, n)
-
-    # Check that solution is close to target
-    @test sum(isapprox.(x, x_sol, atol=1e-6, rtol=1e-2)) == n
-    @test isapprox(f(x), f(result[:raw_solution]), atol=1e-6, rtol=1e-3)
-end

From 7ffd9164e7faee6a5c34fc37abcbfeef72271996 Mon Sep 17 00:00:00 2001
From: LJS42 <1457928939@qq.com>
Date: Wed, 3 Dec 2025 15:02:08 +0100
Subject: [PATCH 07/20] test

---
 src/polytope_blmos.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/polytope_blmos.jl b/src/polytope_blmos.jl
index a0ba57a6a..265423f79 100644
--- a/src/polytope_blmos.jl
+++ b/src/polytope_blmos.jl
@@ -778,4 +778,4 @@ function check_feasibility(lmo::L2normBallBLMO, lb, ub, int_vars, n)
         end
     end
     return OPTIMAL
-end
\ No newline at end of file
+end

From 94cdaf0aacfb06ca2c4657359399c3465969cc8c Mon Sep 17 00:00:00 2001
From: LJS42 <1457928939@qq.com>
Date: Wed, 3 Dec 2025 15:28:23 +0100
Subject: [PATCH 08/20] solve the format problem

---
 test/LMO_test.jl          |  6 +++---
 test/heuristics.jl        |  4 ++--
 test/indicator_test.jl    |  2 +-
 test/interface_test.jl    |  6 +++---
 test/poisson.jl           | 18 +++++++++---------
 test/sparse_regression.jl | 20 ++++++++++----------
 6 files changed, 28 insertions(+), 28 deletions(-)

diff --git a/test/LMO_test.jl b/test/LMO_test.jl
index d7da828d2..e1fcbdcb3 100644
--- a/test/LMO_test.jl
+++ b/test/LMO_test.jl
@@ -140,7 +140,7 @@ end
 end
 
 n = 20
-x_sol = rand(rng, 1:floor(Int, n / 4), n)
+x_sol = rand(rng, 1:floor(Int, n/4), n)
 N = sum(x_sol)
 dir = vcat(fill(1, floor(Int, n / 2)), fill(-1, floor(Int, n / 2)), fill(0, mod(n, 2)))
 diffi = x_sol + 0.3 * dir
@@ -178,7 +178,7 @@ diffi = x_sol + 0.3 * dir
 end
 
 n = 20
-x_sol = rand(rng, 1:floor(Int, n / 4), n)
+x_sol = rand(rng, 1:floor(Int, n/4), n)
 diffi = x_sol + 0.3 * rand(rng, [-1, 1], n)
 
 @testset "Unit Simplex LMO" begin
@@ -199,7 +199,7 @@ diffi = x_sol + 0.3 * rand(rng, [-1, 1], n)
 end
 
 n = 20
-x_sol = rand(1:floor(Int, n / 4), n)
+x_sol = rand(1:floor(Int, n/4), n)
 diffi = x_sol + 0.3 * rand([-1, 1], n)
 
 @testset "Reverse Knapsack LMO" begin
diff --git a/test/heuristics.jl b/test/heuristics.jl
index 705ba72eb..73ef06087 100644
--- a/test/heuristics.jl
+++ b/test/heuristics.jl
@@ -19,7 +19,7 @@ seed = rand(UInt64)
 rng = StableRNG(seed)
 
 n = 20
-x_sol = rand(rng, 1:floor(Int, n / 4), n)
+x_sol = rand(rng, 1:floor(Int, n/4), n)
 N = sum(x_sol)
 dir = vcat(fill(1, floor(Int, n / 2)), fill(-1, floor(Int, n / 2)), fill(0, mod(n, 2)))
 diffi = x_sol + 0.3 * dir
@@ -52,7 +52,7 @@ diffi = x_sol + 0.3 * dir
 end
 
 n = 20
-x_sol = rand(rng, 1:floor(Int, n / 4), n)
+x_sol = rand(rng, 1:floor(Int, n/4), n)
 diffi = x_sol + 0.3 * rand(rng, [-1, 1], n)
 
 @testset "Hyperplane Aware Rounding - Unit Simplex" begin
diff --git a/test/indicator_test.jl b/test/indicator_test.jl
index 0d9721624..d037982b0 100644
--- a/test/indicator_test.jl
+++ b/test/indicator_test.jl
@@ -53,7 +53,7 @@ println("\nIndicator Tests")
     @test Boscia.indicator_present(blmo) == true
 
     function ind_rounding(x)
-        round.(x[n+1:2n])
+        round.(x[(n+1):2n])
         for i in 1:n
             if isapprox(x[n+i], 1.0)
                 x[i] = 0.0
diff --git a/test/interface_test.jl b/test/interface_test.jl
index aa11a17f1..2821cc0a7 100644
--- a/test/interface_test.jl
+++ b/test/interface_test.jl
@@ -332,7 +332,7 @@ Ns = 0.1
         w = @view(θ[1:p])
         b = θ[end]
         storage[1:p] .= 2α .* w
-        storage[p+1:2p] .= 0
+        storage[(p+1):2p] .= 0
         storage[end] = 0
         for i in 1:n
             xi = @view(Xs[:, i])
@@ -348,13 +348,13 @@ Ns = 0.1
     settings.branch_and_bound[:verbose] = false
     x, _, result = Boscia.solve(f, grad!, lmo, settings=settings)
 
-    @test sum(x[p+1:2p]) <= k
+    @test sum(x[(p+1):2p]) <= k
     @test f(x) <= f(result[:raw_solution])
 
     settings = Boscia.create_default_settings()
     settings.branch_and_bound[:start_solution] = x
     x2, _, result = Boscia.solve(f, grad!, lmo, settings=settings)
-    @test sum(x2[p+1:2p]) <= k
+    @test sum(x2[(p+1):2p]) <= k
     @test f(x2) == f(x)
 end
 
diff --git a/test/poisson.jl b/test/poisson.jl
index 47831e13f..8579b9393 100644
--- a/test/poisson.jl
+++ b/test/poisson.jl
@@ -79,7 +79,7 @@ N = 1.0
         w = @view(θ[1:p])
         b = θ[end]
         storage[1:p] .= 2α .* w
-        storage[p+1:2p] .= 0
+        storage[(p+1):2p] .= 0
         storage[end] = 0
         for i in 1:n0
             xi = @view(X0[:, i])
@@ -97,7 +97,7 @@ N = 1.0
     x, _, result = Boscia.solve(f, grad!, lmo, settings=settings)
 
     @test f(x) <= f(result[:raw_solution]) + 1e-6
-    @test sum(x[p+1:2p]) <= k
+    @test sum(x[(p+1):2p]) <= k
 end
 
 @testset "Hybrid branching poisson sparse regression" begin
@@ -151,7 +151,7 @@ end
         w = @view(θ[1:p])
         b = θ[end]
         storage[1:p] .= 2α .* w
-        storage[p+1:2p] .= 0
+        storage[(p+1):2p] .= 0
         storage[end] = 0
         for i in 1:n0
             xi = @view(X0[:, i])
@@ -172,9 +172,9 @@ end
     settings.branch_and_bound[:branching_strategy] = branching_strategy
     settings.tolerances[:fw_epsilon] = 1e-3
     x, _, result = Boscia.solve(f, grad!, lmo, settings=settings)
-    @test sum(x[p+1:2p]) <= k
+    @test sum(x[(p+1):2p]) <= k
     @test f(x) <= f(result[:raw_solution]) + 1e-6
-    @test sum(x[p+1:2p]) <= k
+    @test sum(x[(p+1):2p]) <= k
 end
 
 n0g = 20
@@ -257,7 +257,7 @@ push!(groups, ((k-1)*group_size+1):pg)
         w = @view(θ[1:pg])
         b = θ[end]
         storage[1:pg] .= 2α .* w
-        storage[pg+1:2pg] .= 0
+        storage[(pg+1):2pg] .= 0
         storage[end] = 0
         for i in 1:n0g
             xi = @view(X0g[:, i])
@@ -273,7 +273,7 @@ push!(groups, ((k-1)*group_size+1):pg)
     settings.branch_and_bound[:verbose] = true
     x, _, result = Boscia.solve(f, grad!, lmo, settings=settings)
     @test f(x) <= f(result[:raw_solution]) + 1e-6
-    @test sum(x[p+1:2pg]) <= k
+    @test sum(x[(p+1):2pg]) <= k
 end
 
 
@@ -331,7 +331,7 @@ end
         w = @view(θ[1:pg])
         b = θ[end]
         storage[1:pg] .= 2α .* w
-        storage[pg+1:2pg] .= 0
+        storage[(pg+1):2pg] .= 0
         storage[end] = 0
         for i in 1:n0g
             xi = @view(X0g[:, i])
@@ -353,5 +353,5 @@ end
     x, _, result = Boscia.solve(f, grad!, lmo, settings=settings)
 
     @test f(x) <= f(result[:raw_solution]) + 1e-6
-    @test sum(x[p+1:2pg]) <= k
+    @test sum(x[(p+1):2pg]) <= k
 end
diff --git a/test/sparse_regression.jl b/test/sparse_regression.jl
index 89137ab6a..56895722c 100644
--- a/test/sparse_regression.jl
+++ b/test/sparse_regression.jl
@@ -36,7 +36,7 @@ const M = 2 * var(A)
     MOI.set(o, MOI.Silent(), true)
     MOI.empty!(o)
     x = MOI.add_variables(o, 2p)
-    for i in p+1:2p
+    for i in (p+1):2p
         MOI.add_constraint(o, x[i], MOI.GreaterThan(0.0))
         MOI.add_constraint(o, x[i], MOI.LessThan(1.0))
         MOI.add_constraint(o, x[i], MOI.ZeroOne()) # or MOI.Integer()
@@ -55,13 +55,13 @@ const M = 2 * var(A)
     end
     MOI.add_constraint(
         o,
-        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(ones(p), x[p+1:2p]), 0.0),
+        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(ones(p), x[(p+1):2p]), 0.0),
         MOI.LessThan(k),
     )
     lmo = FrankWolfe.MathOptLMO(o)
 
     function f(x)
-        return sum((y - A * x[1:p]) .^ 2) + lambda_0 * sum(x[p+1:2p]) + lambda_2 * norm(x[1:p])^2
+        return sum((y - A * x[1:p]) .^ 2) + lambda_0 * sum(x[(p+1):2p]) + lambda_2 * norm(x[1:p])^2
     end
     function grad!(storage, x)
         storage .= vcat(
@@ -76,7 +76,7 @@ const M = 2 * var(A)
     settings.branch_and_bound[:time_limit] = 100
     x, _, result = Boscia.solve(f, grad!, lmo, settings=settings)
     # println("Solution: $(x[1:p])")
-    @test sum(x[1+p:2p]) <= k
+    @test sum(x[(1+p):2p]) <= k
     @test f(x) <= f(result[:raw_solution]) + 1e-6
 end
 
@@ -103,7 +103,7 @@ function build_sparse_lmo_grouped()
     MOI.set(o, MOI.Silent(), true)
     MOI.empty!(o)
     x = MOI.add_variables(o, 2p)
-    for i in p+1:2p
+    for i in (p+1):2p
         MOI.add_constraint(o, x[i], MOI.GreaterThan(0.0))
         MOI.add_constraint(o, x[i], MOI.LessThan(1.0))
         MOI.add_constraint(o, x[i], MOI.ZeroOne()) # or MOI.Integer()
@@ -122,7 +122,7 @@ function build_sparse_lmo_grouped()
     end
     MOI.add_constraint(
         o,
-        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(ones(p), x[p+1:2p]), 0.0),
+        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(ones(p), x[(p+1):2p]), 0.0),
         MOI.LessThan(k),
     )
     for i in 1:k_int
@@ -139,13 +139,13 @@ end
 @testset "Sparse Regression Group" begin
     function f(x)
         return sum(abs2, y_g - A_g * x[1:p]) +
-               lambda_0_g * norm(x[p+1:2p])^2 +
+               lambda_0_g * norm(x[(p+1):2p])^2 +
                lambda_2_g * norm(x[1:p])^2
     end
     function grad!(storage, x)
         storage .= vcat(
             2 * (transpose(A_g) * A_g * x[1:p] - transpose(A_g) * y_g + lambda_2_g * x[1:p]),
-            2 * lambda_0_g * x[p+1:2p],
+            2 * lambda_0_g * x[(p+1):2p],
         )
         return storage
     end
@@ -156,7 +156,7 @@ end
     settings.tolerances[:fw_epsilon] = 1e-3
     x, _, result = Boscia.solve(f, grad!, lmo, settings=settings)
 
-    @test sum(x[p+1:2p]) <= k
+    @test sum(x[(p+1):2p]) <= k
     for i in 1:k_int
         @test sum(x[groups[i]]) >= 1
     end
@@ -176,7 +176,7 @@ end
     settings.tolerances[:fw_epsilon] = 1e-3
     settings.tightening[:strong_convexity] = μ
     x2, _, result2 = Boscia.solve(f, grad!, lmo, settings=settings)
-    @test sum(x2[p+1:2p]) <= k
+    @test sum(x2[(p+1):2p]) <= k
     for i in 1:k_int
         @test sum(x2[groups[i]]) >= 1
     end

From 5130785d01789850144a9d02214719cbb1b433bf Mon Sep 17 00:00:00 2001
From: LJS42 <1457928939@qq.com>
Date: Wed, 3 Dec 2025 17:58:37 +0100
Subject: [PATCH 09/20] coverage and format

---
 src/polytope_blmos.jl | 31 ++++++++++++++++++++-----------
 test/LMO_test.jl      | 37 ++++++++++++++++++++++++++++++++++++-
 2 files changed, 56 insertions(+), 12 deletions(-)

diff --git a/src/polytope_blmos.jl b/src/polytope_blmos.jl
index 265423f79..a3616ae1e 100644
--- a/src/polytope_blmos.jl
+++ b/src/polytope_blmos.jl
@@ -726,14 +726,9 @@ BLMO denotes the L2normBall, It is unit ball which means R = 1
 struct L2normBallBLMO <: FrankWolfe.LinearMinimizationOracle end
 
 function bounded_compute_extreme_point(blmo::L2normBallBLMO, d, lb, ub, int_vars; kwargs...)
-    max = abs(d[1])
-    max_idx = 1
-    idx = 0
+    max = 0
+    max_idx = 0
     for idx in int_vars
-        if abs(d[idx]) > max
-            max = -d[idx]
-            max_idx = idx
-        end
         i = findfirst(x -> x == idx, int_vars)
         if lb[i] > 0
             v = zeros(length(d))
@@ -743,15 +738,29 @@ function bounded_compute_extreme_point(blmo::L2normBallBLMO, d, lb, ub, int_vars
             v = zeros(length(d))
             v[idx] = -1
             return v
+        elseif lb[i] == 0 && ub[i] == 0
+            continue
+        end
+        if abs(d[idx]) > max
+            max = abs(d[idx])
+            max_idx = idx
         end
     end
     v = -d
     v[int_vars] .= 0
-    v = v ./ norm(v)
-    prod = dot(v, d)
-    if max < prod
+    if isapprox(norm(v), 0.0; rtol=1e-6)
+        prod = 0
+    else
+        v = v ./ norm(v)
+        prod = dot(v, d)
+    end
+    if -max < prod
         v = zeros(length(d))
-        v[max_idx] = 1
+        if d[max_idx] < 0
+            v[max_idx] = 1
+        elseif d[max_idx] > 0
+            v[max_idx] = -1
+        end
     end
     return v
 end
diff --git a/test/LMO_test.jl b/test/LMO_test.jl
index e1fcbdcb3..2a665a615 100644
--- a/test/LMO_test.jl
+++ b/test/LMO_test.jl
@@ -220,7 +220,7 @@ diffi = x_sol + 0.3 * rand([-1, 1], n)
 end
 
 # Test for L2normBallBLMO
-@testset "L2normBall BLMO" begin
+@testset "L2normBall BLMO continuous" begin
     n = 20
 
     # Generate a solution inside the L2 ball (||x|| <= 1)
@@ -258,3 +258,38 @@ end
     @test sum(isapprox.(x, x_sol, atol=1e-6, rtol=1e-2)) == n  # Should improve or stay similar
     @test isapprox(f(x), f(result[:raw_solution]), atol=1e-6, rtol=1e-3)
 end
+
+@testset "L2normBall BLMO integer" begin
+    n = 20
+
+    # Generate a solution inside the L2 ball (||x|| <= 1)
+    # Round some coordinates to integers for testing
+    num_int = 5
+    int_indices = sort(rand(rng, 1:n, num_int))
+    x_sol = zeros(n)
+    x_sol[int_indices[2]] = 1
+
+    function f(x)
+        return 0.5 * sum((x[i] - x_sol[i])^2 for i in eachindex(x))
+    end
+
+    function grad!(storage, x)
+        @. storage = x - x_sol
+    end
+
+    # Create L2normBallBLMO
+    blmo = Boscia.L2normBallBLMO()
+
+    # Bounds: each variable in [-1, 1] due to L2 ball constraint
+    lower_bounds = fill(-2.0, num_int)
+    upper_bounds = fill(3.0, num_int)
+
+    # Some variables are integer
+    int_vars = collect(int_indices)
+
+    x, _, result = Boscia.solve(f, grad!, blmo, lower_bounds, upper_bounds, int_vars, n)
+
+    # Check objective value
+    @test sum(isapprox.(x, x_sol, atol=1e-6, rtol=1e-2)) == n  # Should improve or stay similar
+    @test isapprox(f(x), f(result[:raw_solution]), atol=1e-6, rtol=1e-3)
+end

From 416e82648a57309414c98858f59e9e40ebdd4aa9 Mon Sep 17 00:00:00 2001
From: LJS42 <1457928939@qq.com>
Date: Wed, 10 Dec 2025 10:07:20 +0100
Subject: [PATCH 10/20] add test and format

---
 test/LMO_test.jl | 37 ++++++++++++++++++++++++++++++++++++-
 1 file changed, 36 insertions(+), 1 deletion(-)

diff --git a/test/LMO_test.jl b/test/LMO_test.jl
index 2a665a615..89ecb21a6 100644
--- a/test/LMO_test.jl
+++ b/test/LMO_test.jl
@@ -259,7 +259,7 @@ end
     @test isapprox(f(x), f(result[:raw_solution]), atol=1e-6, rtol=1e-3)
 end
 
-@testset "L2normBall BLMO integer" begin
+@testset "L2normBall BLMO integer 1" begin
     n = 20
 
     # Generate a solution inside the L2 ball (||x|| <= 1)
@@ -293,3 +293,38 @@ end
     @test sum(isapprox.(x, x_sol, atol=1e-6, rtol=1e-2)) == n  # Should improve or stay similar
     @test isapprox(f(x), f(result[:raw_solution]), atol=1e-6, rtol=1e-3)
 end
+
+@testset "L2normBall BLMO integer 2" begin
+    n = 20
+
+    # Generate a solution inside the L2 ball (||x|| <= 1)
+    # Round some coordinates to integers for testing
+    num_int = 5
+    int_indices = sort(rand(rng, 1:n, num_int))
+    x_sol = zeros(n)
+    x_sol[int_indices[2]] = -1
+
+    function f(x)
+        return 0.5 * sum((x[i] - x_sol[i])^2 for i in eachindex(x))
+    end
+
+    function grad!(storage, x)
+        @. storage = x - x_sol
+    end
+
+    # Create L2normBallBLMO
+    blmo = Boscia.L2normBallBLMO()
+
+    # Bounds: each variable in [-1, 1] due to L2 ball constraint
+    lower_bounds = fill(-2.0, num_int)
+    upper_bounds = fill(3.0, num_int)
+
+    # Some variables are integer
+    int_vars = collect(int_indices)
+
+    x, _, result = Boscia.solve(f, grad!, blmo, lower_bounds, upper_bounds, int_vars, n)
+
+    # Check objective value
+    @test sum(isapprox.(x, x_sol, atol=1e-6, rtol=1e-2)) == n  # Should improve or stay similar
+    @test isapprox(f(x), f(result[:raw_solution]), atol=1e-6, rtol=1e-3)
+end

From 25e0fb63c992fd59725c8b5b1f861830a139cff9 Mon Sep 17 00:00:00 2001
From: LJS42 <1457928939@qq.com>
Date: Wed, 10 Dec 2025 11:18:48 +0100
Subject: [PATCH 11/20] format

---
 Project.toml                          |   2 +
 examples/birkhoff_decomposition.jl    |   4 +-
 examples/docs-01-network-design.jl    | 158 ++++++++++++++++----------
 examples/docs-02-graph-isomorphism.jl |   2 +-
 examples/docs-03-optimal-design.jl    |  46 ++++----
 examples/lasso.jl                     |   6 +-
 examples/oed_utils.jl                 |   1 -
 examples/plot_utilities.jl            |  39 +++----
 examples/poisson_reg.jl               |   4 +-
 examples/quadratic_over_birkhoff.jl   |   6 +-
 examples/sparse_reg.jl                |   8 +-
 test/LMO_test.jl                      |   6 +-
 test/heuristics.jl                    |   4 +-
 13 files changed, 161 insertions(+), 125 deletions(-)

diff --git a/Project.toml b/Project.toml
index d77196ac5..4f9af0ed7 100644
--- a/Project.toml
+++ b/Project.toml
@@ -7,6 +7,7 @@ authors = ["ZIB IOL"]
 Bonobo = "f7b14807-3d4d-461a-888a-05dd4bca8bc3"
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
 FrankWolfe = "f55ce6ea-fdc5-4628-88c5-0087fe54bd30"
+JuliaFormatter = "98e50ef6-434e-11e9-1051-2b60c6c9e899"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee"
 MathOptSetDistances = "3b969827-a86c-476c-9527-bb6f1a8fbad5"
@@ -34,6 +35,7 @@ DoubleFloats = "1"
 FrankWolfe = "0.6"
 Graphs = "1.13"
 HiGHS = "1"
+JuliaFormatter = "1.0.62"
 LinearAlgebra = "1.6"
 MathOptInterface = "1"
 MathOptSetDistances = "0.2"
diff --git a/examples/birkhoff_decomposition.jl b/examples/birkhoff_decomposition.jl
index 291fba102..d43575771 100644
--- a/examples/birkhoff_decomposition.jl
+++ b/examples/birkhoff_decomposition.jl
@@ -56,10 +56,10 @@ end
 function grad!(storage, x)
     storage .= 0
     for j in 1:k
-        Sk = reshape(@view(storage[(j-1)*n^2+1:j*n^2]), n, n)
+        Sk = reshape(@view(storage[((j-1)*n^2+1):(j*n^2)]), n, n)
         @. Sk = -Xstar
         for m in 1:k
-            Yk = reshape(@view(x[(m-1)*n^2+1:m*n^2]), n, n)
+            Yk = reshape(@view(x[((m-1)*n^2+1):(m*n^2)]), n, n)
             @. Sk += Yk
         end
     end
diff --git a/examples/docs-01-network-design.jl b/examples/docs-01-network-design.jl
index 8c1bbbf3b..c84c57661 100644
--- a/examples/docs-01-network-design.jl
+++ b/examples/docs-01-network-design.jl
@@ -41,14 +41,14 @@ using SparseArrays
 using LinearAlgebra
 import MathOptInterface
 const MOI = MathOptInterface
-using HiGHS 
+using HiGHS
 
 println("\nDocumentation Example 01: Network Design Problem")
 
 # The graph structure is shown below.
 mutable struct NetworkData
     num_nodes::Int
-    num_edges::Int  
+    num_edges::Int
     init_nodes::Vector{Int}
     term_nodes::Vector{Int}
     free_flow_time::Vector{Float64}
@@ -73,8 +73,18 @@ function load_braess_network()
     b = [0.1, 0.1, 0.1, 0.1, 3.0, 0.1, 0.1, 0.1, 0.1, 0.1, 3.0, 0.1]
     power = [2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0]
     travel_demand = [0.0 0.0 1.0; 0.0 0.0 1.0; 0.0 0.0 0.0]
-    return NetworkData(8, length(init_nodes), init_nodes, term_nodes, free_flow_time,
-                      capacity, b, power, travel_demand, 3)
+    return NetworkData(
+        8,
+        length(init_nodes),
+        init_nodes,
+        term_nodes,
+        free_flow_time,
+        capacity,
+        b,
+        power,
+        travel_demand,
+        3,
+    )
 end
 
 # ## Direct modelling via MathOptInterface
@@ -106,15 +116,15 @@ function build_moi_model(net_data, removed_edges, use_big_m=true)
     end
     edge_list = [(net_data.init_nodes[i], net_data.term_nodes[i]) for i in 1:num_edges]
     edge_dict = Dict(edge_list[i] => i for i in eachindex(edge_list))
-    incoming = Dict{Int, Vector{Int}}()
-    outgoing = Dict{Int, Vector{Int}}()
-    
+    incoming = Dict{Int,Vector{Int}}()
+    outgoing = Dict{Int,Vector{Int}}()
+
     for (idx, (src, dst)) in enumerate(edge_list)
         if !haskey(outgoing, src)
             outgoing[src] = Int[]
         end
         push!(outgoing[src], idx)
-        
+
         if !haskey(incoming, dst)
             incoming[dst] = Int[]
         end
@@ -125,12 +135,12 @@ function build_moi_model(net_data, removed_edges, use_big_m=true)
             terms = MOI.ScalarAffineTerm{Float64}[]
             if haskey(outgoing, node)
                 for edge_idx in outgoing[node]
-                    push!(terms, MOI.ScalarAffineTerm(1.0, x[(dest-1)*num_edges + edge_idx]))
+                    push!(terms, MOI.ScalarAffineTerm(1.0, x[(dest-1)*num_edges+edge_idx]))
                 end
             end
             if haskey(incoming, node)
                 for edge_idx in incoming[node]
-                    push!(terms, MOI.ScalarAffineTerm(-1.0, x[(dest-1)*num_edges + edge_idx]))
+                    push!(terms, MOI.ScalarAffineTerm(-1.0, x[(dest-1)*num_edges+edge_idx]))
                 end
             end
             if node == dest
@@ -140,19 +150,15 @@ function build_moi_model(net_data, removed_edges, use_big_m=true)
             else
                 rhs = 0.0
             end
-            MOI.add_constraint(optimizer, 
-                             MOI.ScalarAffineFunction(terms, 0.0),
-                             MOI.EqualTo(rhs))
+            MOI.add_constraint(optimizer, MOI.ScalarAffineFunction(terms, 0.0), MOI.EqualTo(rhs))
         end
     end
     for edge_idx in 1:num_edges
         terms = [MOI.ScalarAffineTerm(1.0, x_agg[edge_idx])]
         for dest in 1:num_zones
-            push!(terms, MOI.ScalarAffineTerm(-1.0, x[(dest-1)*num_edges + edge_idx]))
+            push!(terms, MOI.ScalarAffineTerm(-1.0, x[(dest-1)*num_edges+edge_idx]))
         end
-        MOI.add_constraint(optimizer,
-                         MOI.ScalarAffineFunction(terms, 0.0),
-                         MOI.EqualTo(0.0))
+        MOI.add_constraint(optimizer, MOI.ScalarAffineFunction(terms, 0.0), MOI.EqualTo(0.0))
     end
     max_flow = 1.5 * sum(net_data.travel_demand)
     for (y_idx, edge) in enumerate(removed_edges)
@@ -162,21 +168,26 @@ function build_moi_model(net_data, removed_edges, use_big_m=true)
             if use_big_m
                 terms = [
                     MOI.ScalarAffineTerm(1.0, x[var_idx]),
-                    MOI.ScalarAffineTerm(-max_flow, y[y_idx])
+                    MOI.ScalarAffineTerm(-max_flow, y[y_idx]),
                 ]
-                MOI.add_constraint(optimizer,
-                                 MOI.ScalarAffineFunction(terms, 0.0),
-                                 MOI.LessThan(0.0))
+                MOI.add_constraint(
+                    optimizer,
+                    MOI.ScalarAffineFunction(terms, 0.0),
+                    MOI.LessThan(0.0),
+                )
             else
                 indicator_func = MOI.VectorAffineFunction(
                     [
                         MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(1.0, y[y_idx])),
-                        MOI.VectorAffineTerm(2, MOI.ScalarAffineTerm(1.0, x[var_idx]))
+                        MOI.VectorAffineTerm(2, MOI.ScalarAffineTerm(1.0, x[var_idx])),
                     ],
-                    [0.0, 0.0]
+                    [0.0, 0.0],
+                )
+                MOI.add_constraint(
+                    optimizer,
+                    indicator_func,
+                    MOI.Indicator{MOI.ACTIVATE_ON_ZERO}(MOI.EqualTo(0.0)),
                 )
-                MOI.add_constraint(optimizer, indicator_func,
-                                 MOI.Indicator{MOI.ACTIVATE_ON_ZERO}(MOI.EqualTo(0.0)))
             end
         end
     end
@@ -213,7 +224,7 @@ function build_objective_and_gradient(net_data, removed_edges, cost_per_edge)
         end
         design_start = num_zones * num_edges + num_edges + 1
         for i in 1:num_removed
-            total += cost_per_edge[i] * x[design_start + i - 1]
+            total += cost_per_edge[i] * x[design_start+i-1]
         end
         return total
     end
@@ -229,16 +240,16 @@ function build_objective_and_gradient(net_data, removed_edges, cost_per_edge)
             b = net_data.b[i]
             cap = net_data.capacity[i]
             p = net_data.power[i]
-            storage[agg_start + i - 1] = t0 * (1 + b * flow^p / cap^p)
+            storage[agg_start+i-1] = t0 * (1 + b * flow^p / cap^p)
         end
         for dest in 1:num_zones
             for edge in 1:num_edges
-                storage[(dest - 1) * num_edges + edge] = storage[agg_start + edge - 1]
+                storage[(dest-1)*num_edges+edge] = storage[agg_start+edge-1]
             end
         end
         design_start = num_zones * num_edges + num_edges + 1
         for i in 1:num_removed
-            storage[design_start + i - 1] = cost_per_edge[i]
+            storage[design_start+i-1] = cost_per_edge[i]
         end
         return storage
     end
@@ -288,8 +299,8 @@ x_moi, _, result_moi = Boscia.solve(f_moi, grad_moi!, lmo_moi, settings=settings
 struct ShortestPathLMO <: FrankWolfe.LinearMinimizationOracle
     graph::Graphs.SimpleDiGraph{Int}
     net_data::NetworkData
-    link_dic::SparseMatrixCSC{Int, Int}
-    edge_list::Vector{Tuple{Int, Int}}
+    link_dic::SparseMatrixCSC{Int,Int}
+    edge_list::Vector{Tuple{Int,Int}}
 end
 
 # Add demand to flow vector following shortest path
@@ -298,14 +309,14 @@ function add_demand_to_path!(x, demand, state, origin, destination, link_dic, ed
     parent = -1
     edge_count = length(edge_list)
     agg_start = edge_count * num_zones
-    
+
     while parent != origin && origin != destination && current != 0
         parent = state.parents[current]
         if parent != 0
             link_idx = link_dic[parent, current]
             if link_idx != 0
-                x[(destination - 1) * edge_count + link_idx] += demand
-                x[agg_start + link_idx] += demand
+                x[(destination-1)*edge_count+link_idx] += demand
+                x[agg_start+link_idx] += demand
             end
         end
         current = parent
@@ -316,28 +327,40 @@ end
 function all_or_nothing_assignment(travel_time_vector, net_data, graph, link_dic, edge_list)
     num_zones = net_data.num_zones
     edge_count = net_data.num_edges
-    travel_time = travel_time_vector[num_zones * edge_count + 1 : (num_zones + 1) * edge_count]
+    travel_time = travel_time_vector[(num_zones*edge_count+1):((num_zones+1)*edge_count)]
     x = zeros(length(travel_time_vector))
-    
+
     for origin in 1:num_zones
         state = Graphs.dijkstra_shortest_paths(graph, origin)
-        
+
         for destination in 1:num_zones
             demand = net_data.travel_demand[origin, destination]
             if demand > 0
-                add_demand_to_path!(x, demand, state, origin, destination, 
-                                  link_dic, edge_list, num_zones)
+                add_demand_to_path!(
+                    x,
+                    demand,
+                    state,
+                    origin,
+                    destination,
+                    link_dic,
+                    edge_list,
+                    num_zones,
+                )
             end
         end
     end
-    
+
     return x
 end
 
-function Boscia.bounded_compute_extreme_point(lmo::ShortestPathLMO, direction, 
-                                               lower_bounds, upper_bounds, int_vars)
-    x = all_or_nothing_assignment(direction, lmo.net_data, lmo.graph, 
-                                  lmo.link_dic, lmo.edge_list)
+function Boscia.bounded_compute_extreme_point(
+    lmo::ShortestPathLMO,
+    direction,
+    lower_bounds,
+    upper_bounds,
+    int_vars,
+)
+    x = all_or_nothing_assignment(direction, lmo.net_data, lmo.graph, lmo.link_dic, lmo.edge_list)
     for (i, var_idx) in enumerate(int_vars)
         if direction[var_idx] < 0
             x[var_idx] = upper_bounds[i]
@@ -369,14 +392,19 @@ end
 # The gradient function computes derivatives of the objective with respect to:
 # - Aggregate flows: d/d(flow) of BPR function + penalty gradient w.r.t. flows
 # - Design variables: cost_per_edge[i] + penalty gradient w.r.t. design variables
-function build_objective_and_gradient_with_penalty(net_data, removed_edges, cost_per_edge, 
-                                                    penalty_weight=1e6, penalty_exponent=2.0)
+function build_objective_and_gradient_with_penalty(
+    net_data,
+    removed_edges,
+    cost_per_edge,
+    penalty_weight=1e6,
+    penalty_exponent=2.0,
+)
     num_zones = net_data.num_zones
     num_edges = net_data.num_edges
     num_removed = length(removed_edges)
     edge_list = [(net_data.init_nodes[i], net_data.term_nodes[i]) for i in 1:num_edges]
-    removed_edge_indices = [findfirst(e -> e == removed_edge, edge_list) 
-                            for removed_edge in removed_edges]
+    removed_edge_indices =
+        [findfirst(e -> e == removed_edge, edge_list) for removed_edge in removed_edges]
     max_flow = 1.5 * sum(net_data.travel_demand)
     function f(x)
         x = max.(x, 0.0)
@@ -394,11 +422,11 @@ function build_objective_and_gradient_with_penalty(net_data, removed_edges, cost
         end
         design_start = num_zones * num_edges + num_edges + 1
         for i in 1:num_removed
-            total += cost_per_edge[i] * x[design_start + i - 1]
+            total += cost_per_edge[i] * x[design_start+i-1]
         end
         for (y_idx, edge_idx) in enumerate(removed_edge_indices)
             if edge_idx !== nothing
-                y_val = x[design_start + y_idx - 1]
+                y_val = x[design_start+y_idx-1]
                 for dest in 1:num_zones
                     flow_idx = (dest - 1) * num_edges + edge_idx
                     flow_val = x[flow_idx]
@@ -421,28 +449,29 @@ function build_objective_and_gradient_with_penalty(net_data, removed_edges, cost
             b = net_data.b[i]
             cap = net_data.capacity[i]
             p = net_data.power[i]
-            storage[agg_start + i - 1] = t0 * (1 + b * flow^p / cap^p)
+            storage[agg_start+i-1] = t0 * (1 + b * flow^p / cap^p)
         end
         for dest in 1:num_zones
             for edge in 1:num_edges
-                storage[(dest - 1) * num_edges + edge] = storage[agg_start + edge - 1]
+                storage[(dest-1)*num_edges+edge] = storage[agg_start+edge-1]
             end
         end
         design_start = num_zones * num_edges + num_edges + 1
         for i in 1:num_removed
-            storage[design_start + i - 1] = cost_per_edge[i]
+            storage[design_start+i-1] = cost_per_edge[i]
         end
         for (y_idx, edge_idx) in enumerate(removed_edge_indices)
             if edge_idx !== nothing
-                y_val = x[design_start + y_idx - 1]
+                y_val = x[design_start+y_idx-1]
                 for dest in 1:num_zones
                     flow_idx = (dest - 1) * num_edges + edge_idx
                     flow_val = x[flow_idx]
                     violation = max(0.0, flow_val - max_flow * y_val)
                     if violation > 1e-10
-                        grad_coeff = penalty_weight * penalty_exponent * violation^(penalty_exponent - 1)
+                        grad_coeff =
+                            penalty_weight * penalty_exponent * violation^(penalty_exponent - 1)
                         storage[flow_idx] += grad_coeff
-                        storage[design_start + y_idx - 1] += grad_coeff * (-max_flow)
+                        storage[design_start+y_idx-1] += grad_coeff * (-max_flow)
                     end
                 end
             end
@@ -465,8 +494,7 @@ for i in 1:net_data.num_edges
     push!(edge_list_custom, (net_data.init_nodes[i], net_data.term_nodes[i]))
 end
 
-link_dic = sparse(net_data.init_nodes, net_data.term_nodes, 
-                 collect(1:net_data.num_edges))
+link_dic = sparse(net_data.init_nodes, net_data.term_nodes, collect(1:net_data.num_edges))
 
 custom_lmo = ShortestPathLMO(graph, net_data, link_dic, edge_list_custom)
 
@@ -476,19 +504,25 @@ num_edges = net_data.num_edges
 num_removed = length(removed_edges)
 total_vars = num_zones * num_edges + num_edges + num_removed
 
-int_vars = collect((num_zones * num_edges + num_edges + 1):total_vars) # last num_removed variables
+int_vars = collect((num_zones*num_edges+num_edges+1):total_vars) # last num_removed variables
 lower_bounds = zeros(Float64, num_removed)  # Binary: lower bound = 0
 upper_bounds = ones(Float64, num_removed)   # Binary: upper bound = 1
 
 # To have Boscia handle the bounds, we need to wrap our LMO in an instance of `ManagedLMO`.
 bounded_lmo = Boscia.ManagedLMO(custom_lmo, lower_bounds, upper_bounds, int_vars, total_vars)
 
-f_custom, grad_custom! = build_objective_and_gradient_with_penalty(net_data, removed_edges, cost_per_edge, 
-                                                      penalty_weight, penalty_exponent)
+f_custom, grad_custom! = build_objective_and_gradient_with_penalty(
+    net_data,
+    removed_edges,
+    cost_per_edge,
+    penalty_weight,
+    penalty_exponent,
+)
 
 settings_custom = Boscia.create_default_settings()
 settings_custom.branch_and_bound[:verbose] = true
 
-x_custom, _, result_custom = Boscia.solve(f_custom, grad_custom!, bounded_lmo, settings=settings_custom)
+x_custom, _, result_custom =
+    Boscia.solve(f_custom, grad_custom!, bounded_lmo, settings=settings_custom)
 
 @show x_custom
diff --git a/examples/docs-02-graph-isomorphism.jl b/examples/docs-02-graph-isomorphism.jl
index fd1c3c4c3..40d5d7e53 100644
--- a/examples/docs-02-graph-isomorphism.jl
+++ b/examples/docs-02-graph-isomorphism.jl
@@ -256,6 +256,6 @@ println("Certificate verified: graphs are isomorphic (A ≈ X' * B * X)")
 # ```
 
 B = randomNonIsomorphic(A)
-x, _, result = Boscia.solve(f, grad!, blmo, settings = settings)
+x, _, result = Boscia.solve(f, grad!, blmo, settings=settings)
 @assert result[:dual_bound] > 0.0
 println("Graphs are not isomorphic (lower bound > 0)")
diff --git a/examples/docs-03-optimal-design.jl b/examples/docs-03-optimal-design.jl
index ddfea5262..20092d790 100644
--- a/examples/docs-03-optimal-design.jl
+++ b/examples/docs-03-optimal-design.jl
@@ -54,7 +54,7 @@ function f_a(x)
     X = Symmetric(X)
     U = cholesky(X)
     X_inv = U \ I
-    return LinearAlgebra.tr(X_inv) 
+    return LinearAlgebra.tr(X_inv)
 end
 
 function grad_a!(storage, x)
@@ -62,7 +62,7 @@ function grad_a!(storage, x)
     X = Symmetric(X * X)
     F = cholesky(X)
     for i in 1:length(x)
-        storage[i] = LinearAlgebra.tr(-(F \ A[i, :]) * transpose(A[i, :])) 
+        storage[i] = LinearAlgebra.tr(-(F \ A[i, :]) * transpose(A[i, :]))
     end
     return storage
 end
@@ -123,31 +123,31 @@ end
 # If $x$ does not yet satisfy the knapsack constraint, we increase the values of $X$, first by sampling 
 # from the linearly independent rows and then by adding 1 to the smallest value of $x$ while respecting the upper bounds $u$.
 function linearly_independent_rows(A; u=fill(1, size(A, 1)))
-S = []
-m, n = size(A)
-for i in 1:m
-    if iszero(u[i])
-        continue
-    end
-    S_i = vcat(S, i)
-    if rank(A[S_i, :]) == length(S_i)
-        S = S_i
-    end
-    if length(S) == n # we only n linearly independent points
-        return S
+    S = []
+    m, n = size(A)
+    for i in 1:m
+        if iszero(u[i])
+            continue
+        end
+        S_i = vcat(S, i)
+        if rank(A[S_i, :]) == length(S_i)
+            S = S_i
+        end
+        if length(S) == n # we only n linearly independent points
+            return S
+        end
     end
-end
-return S # then x= zeros(m) and x[S] = 1
+    return S # then x= zeros(m) and x[S] = 1
 end
 function add_to_min(x, u)
-perm = sortperm(x)
-for i in perm
-    if x[i] < u[i]
-        x[i] += 1
-        break
+    perm = sortperm(x)
+    for i in perm
+        if x[i] < u[i]
+            x[i] += 1
+            break
+        end
     end
-end
-return x
+    return x
 end
 function domain_point(local_bounds)
     lb = fill(0.0, m)
diff --git a/examples/lasso.jl b/examples/lasso.jl
index 5a369bcfc..d253ae4a4 100644
--- a/examples/lasso.jl
+++ b/examples/lasso.jl
@@ -39,7 +39,7 @@ const A_g = rand(rng, Float64, n, p)
 k_int = convert(Int64, k)
 
 for i in 1:k_int
-    for _ in 1:group_size-1
+    for _ in 1:(group_size-1)
         β_sol[rand(rng, ((i-1)*group_size+1):(i*group_size))] = 0
     end
 end
@@ -114,7 +114,7 @@ push!(groups, ((k_int-1)*group_size+1):p)
 
     function f(x)
         return sum((y_g - A_g * x[1:p]) .^ 2) +
-               lambda_0_g * sum(x[p+1:2p]) +
+               lambda_0_g * sum(x[(p+1):2p]) +
                lambda_2_g * FrankWolfe.norm(x[1:p])^2
     end
     function grad!(storage, x)
@@ -132,7 +132,7 @@ push!(groups, ((k_int-1)*group_size+1):p)
     x, _, result = Boscia.solve(f, grad!, blmo, settings=settings)
 
     # println("Solution: $(x[1:p])")
-    z = x[p+1:2p]
+    z = x[(p+1):2p]
     @test sum(z) <= k
     for i in 1:k_int
         @test sum(z[groups[i]]) >= 1
diff --git a/examples/oed_utils.jl b/examples/oed_utils.jl
index d21c9f849..44ea71344 100644
--- a/examples/oed_utils.jl
+++ b/examples/oed_utils.jl
@@ -350,4 +350,3 @@ function build_domain_oracle2(A, n)
         return minimum(eigvals(X)) > sqrt(eps())
     end
 end
-
diff --git a/examples/plot_utilities.jl b/examples/plot_utilities.jl
index c246a0766..8c1a5827d 100644
--- a/examples/plot_utilities.jl
+++ b/examples/plot_utilities.jl
@@ -57,19 +57,19 @@ plot_bounds_progress(result, "output.png",
 function plot_bounds_progress(
     result::Dict,
     filename::String;
-    title_prefix::String = "",
-    font_family::String = "serif",
-    font_size::Int = 11,
-    use_latex::Bool = true,
-    latex_preamble::String = "\\usepackage{charter}\\usepackage[charter]{mathdesign}",
-    lower_color::String = "C0",
-    upper_color::String = "C1",
-    linewidth::Real = 2,
-    figsize::Tuple{Real,Real} = (12, 4),
-    dpi::Int = 300,
-    show_grid::Bool = true,
-    grid_alpha::Real = 0.3,
-    legend_loc::String = "best",
+    title_prefix::String="",
+    font_family::String="serif",
+    font_size::Int=11,
+    use_latex::Bool=true,
+    latex_preamble::String="\\usepackage{charter}\\usepackage[charter]{mathdesign}",
+    lower_color::String="C0",
+    upper_color::String="C1",
+    linewidth::Real=2,
+    figsize::Tuple{Real,Real}=(12, 4),
+    dpi::Int=300,
+    show_grid::Bool=true,
+    grid_alpha::Real=0.3,
+    legend_loc::String="best",
 )
     # Set up fonts
     PyPlot.rc("text", usetex=use_latex)
@@ -145,17 +145,19 @@ setup_plot_font(use_latex=false, font_family="sans-serif", font_size=12)
 ```
 """
 function setup_plot_font(;
-    font_family::String = "serif",
-    font_size::Int = 11,
-    use_latex::Bool = true,
-    latex_preamble::String = "\\usepackage{charter}\\usepackage[charter]{mathdesign}",
+    font_family::String="serif",
+    font_size::Int=11,
+    use_latex::Bool=true,
+    latex_preamble::String="\\usepackage{charter}\\usepackage[charter]{mathdesign}",
 )
     PyPlot.rc("text", usetex=use_latex)
     if use_latex
         PyPlot.rc("text.latex", preamble=latex_preamble)
     end
     PyPlot.rc("font", family=font_family, size=font_size)
-    println("Configured plotting with font_family=$font_family, font_size=$font_size, use_latex=$use_latex")
+    return println(
+        "Configured plotting with font_family=$font_family, font_size=$font_size, use_latex=$use_latex",
+    )
 end
 
 """
@@ -185,4 +187,3 @@ function find_fonts(pattern::String)
     fonts = list_available_fonts()
     return filter(f -> occursin(Regex(pattern, "i"), f), fonts)
 end
-
diff --git a/examples/poisson_reg.jl b/examples/poisson_reg.jl
index 6309b011c..4b483e877 100644
--- a/examples/poisson_reg.jl
+++ b/examples/poisson_reg.jl
@@ -101,7 +101,7 @@ Ns = 0.10
         w = @view(θ[1:p])
         b = θ[end]
         storage[1:p] .= 2α .* w
-        storage[p+1:2p] .= 0
+        storage[(p+1):2p] .= 0
         storage[end] = 0
         for i in 1:n
             xi = @view(Xs[:, i])
@@ -119,5 +119,5 @@ Ns = 0.10
     #@show x
     @show result[:raw_solution]
     @test f(x) <= f(result[:raw_solution]) + 1e-6
-    @test sum(x[p+1:2p]) <= k
+    @test sum(x[(p+1):2p]) <= k
 end
diff --git a/examples/quadratic_over_birkhoff.jl b/examples/quadratic_over_birkhoff.jl
index a7cf97f63..1b16b4ac4 100644
--- a/examples/quadratic_over_birkhoff.jl
+++ b/examples/quadratic_over_birkhoff.jl
@@ -70,9 +70,9 @@ end
     f, grad! = build_objective(n)
 
     x = zeros(n, n)
-    int_vars = collect(1:n^2)
+    int_vars = collect(1:(n^2))
     @testset "Birkhoff BLMO (BPCG)" begin
-        lmo = CLO.BirkhoffLMO(n, collect(1:n^2))
+        lmo = CLO.BirkhoffLMO(n, collect(1:(n^2)))
 
         settings = Boscia.create_default_settings()
         settings.branch_and_bound[:verbose] = true
@@ -83,7 +83,7 @@ end
 
     x_dicg = zeros(n, n)
     @testset "Birkhoff BLMO (DICG)" begin
-        lmo = CLO.BirkhoffLMO(n, collect(1:n^2))
+        lmo = CLO.BirkhoffLMO(n, collect(1:(n^2)))
 
         settings = Boscia.create_default_settings()
         settings.branch_and_bound[:verbose] = true
diff --git a/examples/sparse_reg.jl b/examples/sparse_reg.jl
index 5e42f88fe..12581c0ee 100644
--- a/examples/sparse_reg.jl
+++ b/examples/sparse_reg.jl
@@ -49,7 +49,7 @@ const M = 2 * var(A)
     MOI.set(o, MOI.Silent(), true)
     MOI.empty!(o)
     x = MOI.add_variables(o, 2p)
-    for i in p+1:2p
+    for i in (p+1):2p
         MOI.add_constraint(o, x[i], MOI.GreaterThan(0.0))
         MOI.add_constraint(o, x[i], MOI.LessThan(1.0))
         MOI.add_constraint(o, x[i], MOI.ZeroOne()) # or MOI.Integer()
@@ -68,19 +68,19 @@ const M = 2 * var(A)
     end
     MOI.add_constraint(
         o,
-        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(ones(p), x[p+1:2p]), 0.0),
+        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(ones(p), x[(p+1):2p]), 0.0),
         MOI.LessThan(k),
     )
     blmo = Boscia.MathOptBLMO(o)
 
     function f(x)
         xv = @view(x[1:p])
-        return norm(y - A * xv)^2 + lambda_0 * sum(x[p+1:2p]) + lambda_2 * norm(xv)^2
+        return norm(y - A * xv)^2 + lambda_0 * sum(x[(p+1):2p]) + lambda_2 * norm(xv)^2
     end
 
     function grad!(storage, x)
         storage[1:p] .= 2 * (transpose(A) * A * x[1:p] - transpose(A) * y + lambda_2 * x[1:p])
-        storage[p+1:2p] .= lambda_0
+        storage[(p+1):2p] .= lambda_0
         return storage
     end
 
diff --git a/test/LMO_test.jl b/test/LMO_test.jl
index 89ecb21a6..f972bcdc7 100644
--- a/test/LMO_test.jl
+++ b/test/LMO_test.jl
@@ -140,7 +140,7 @@ end
 end
 
 n = 20
-x_sol = rand(rng, 1:floor(Int, n/4), n)
+x_sol = rand(rng, 1:floor(Int, n / 4), n)
 N = sum(x_sol)
 dir = vcat(fill(1, floor(Int, n / 2)), fill(-1, floor(Int, n / 2)), fill(0, mod(n, 2)))
 diffi = x_sol + 0.3 * dir
@@ -178,7 +178,7 @@ diffi = x_sol + 0.3 * dir
 end
 
 n = 20
-x_sol = rand(rng, 1:floor(Int, n/4), n)
+x_sol = rand(rng, 1:floor(Int, n / 4), n)
 diffi = x_sol + 0.3 * rand(rng, [-1, 1], n)
 
 @testset "Unit Simplex LMO" begin
@@ -199,7 +199,7 @@ diffi = x_sol + 0.3 * rand(rng, [-1, 1], n)
 end
 
 n = 20
-x_sol = rand(1:floor(Int, n/4), n)
+x_sol = rand(1:floor(Int, n / 4), n)
 diffi = x_sol + 0.3 * rand([-1, 1], n)
 
 @testset "Reverse Knapsack LMO" begin
diff --git a/test/heuristics.jl b/test/heuristics.jl
index 73ef06087..705ba72eb 100644
--- a/test/heuristics.jl
+++ b/test/heuristics.jl
@@ -19,7 +19,7 @@ seed = rand(UInt64)
 rng = StableRNG(seed)
 
 n = 20
-x_sol = rand(rng, 1:floor(Int, n/4), n)
+x_sol = rand(rng, 1:floor(Int, n / 4), n)
 N = sum(x_sol)
 dir = vcat(fill(1, floor(Int, n / 2)), fill(-1, floor(Int, n / 2)), fill(0, mod(n, 2)))
 diffi = x_sol + 0.3 * dir
@@ -52,7 +52,7 @@ diffi = x_sol + 0.3 * dir
 end
 
 n = 20
-x_sol = rand(rng, 1:floor(Int, n/4), n)
+x_sol = rand(rng, 1:floor(Int, n / 4), n)
 diffi = x_sol + 0.3 * rand(rng, [-1, 1], n)
 
 @testset "Hyperplane Aware Rounding - Unit Simplex" begin

From 74ba511451ec9aeaa8050ee886bd3a201c43416d Mon Sep 17 00:00:00 2001
From: LJS42 <1457928939@qq.com>
Date: Wed, 10 Dec 2025 12:06:05 +0100
Subject: [PATCH 12/20] dep

---
 Project.toml | 2 --
 1 file changed, 2 deletions(-)

diff --git a/Project.toml b/Project.toml
index 4f9af0ed7..d77196ac5 100644
--- a/Project.toml
+++ b/Project.toml
@@ -7,7 +7,6 @@ authors = ["ZIB IOL"]
 Bonobo = "f7b14807-3d4d-461a-888a-05dd4bca8bc3"
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
 FrankWolfe = "f55ce6ea-fdc5-4628-88c5-0087fe54bd30"
-JuliaFormatter = "98e50ef6-434e-11e9-1051-2b60c6c9e899"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee"
 MathOptSetDistances = "3b969827-a86c-476c-9527-bb6f1a8fbad5"
@@ -35,7 +34,6 @@ DoubleFloats = "1"
 FrankWolfe = "0.6"
 Graphs = "1.13"
 HiGHS = "1"
-JuliaFormatter = "1.0.62"
 LinearAlgebra = "1.6"
 MathOptInterface = "1"
 MathOptSetDistances = "0.2"

From 1ee885820142d6dc985e861e6cb0f23346945c1e Mon Sep 17 00:00:00 2001
From: LJS42 <1457928939@qq.com>
Date: Wed, 17 Dec 2025 16:29:41 +0100
Subject: [PATCH 13/20] fix

---
 src/polytope_blmos.jl | 19 ++++++----
 test/LMO_test.jl      | 82 ++++++++++++-------------------------------
 test/heuristics.jl    |  4 +--
 3 files changed, 37 insertions(+), 68 deletions(-)

diff --git a/src/polytope_blmos.jl b/src/polytope_blmos.jl
index a3616ae1e..b91298553 100644
--- a/src/polytope_blmos.jl
+++ b/src/polytope_blmos.jl
@@ -723,11 +723,17 @@ end
 BLMO denotes the L2normBall, It is unit ball which means R = 1
 """
 
-struct L2normBallBLMO <: FrankWolfe.LinearMinimizationOracle end
-
-function bounded_compute_extreme_point(blmo::L2normBallBLMO, d, lb, ub, int_vars; kwargs...)
+function bounded_compute_extreme_point(
+    lmo::FrankWolfe.LpNormBallLMO{T,2},
+    d,
+    lb,
+    ub,
+    int_vars;
+    kwargs...,
+) where {T}
     max = 0
     max_idx = 0
+    #the extreme point would be either a point with 1 in a integer entry or in the hyperplane in all integer entry equal to 0
     for idx in int_vars
         i = findfirst(x -> x == idx, int_vars)
         if lb[i] > 0
@@ -765,7 +771,7 @@ function bounded_compute_extreme_point(blmo::L2normBallBLMO, d, lb, ub, int_vars
     return v
 end
 
-function is_simple_linear_feasible(lmo::L2normBallBLMO, v)
+function is_simple_linear_feasible(lmo::FrankWolfe.LpNormBallLMO{T,2}, v) where {T}
     if norm(v) > 1 + 1e-6
         @debug "norm(v) : $(norm(v)) > 1"
         return false
@@ -773,14 +779,15 @@ function is_simple_linear_feasible(lmo::L2normBallBLMO, v)
     return true
 end
 
-function check_feasibility(lmo::L2normBallBLMO, lb, ub, int_vars, n)
+function check_feasibility(lmo::FrankWolfe.LpNormBallLMO{T,2}, lb, ub, int_vars, n) where {T}
     if any(lb .> 1) || any(ub .< -1)
         return INFEASIBLE
     end
     count = 0
+    #If there is at least two entry have to larger or equal to 1 than that is definitely infeasible
     for idx in 1:length(int_vars)
         if !(lb[idx] <= 0 <= ub[idx])
-            count = count + 1
+            count += count
         end
         if count > 1
             return INFEASIBLE
diff --git a/test/LMO_test.jl b/test/LMO_test.jl
index f972bcdc7..2bbe7dc63 100644
--- a/test/LMO_test.jl
+++ b/test/LMO_test.jl
@@ -140,7 +140,7 @@ end
 end
 
 n = 20
-x_sol = rand(rng, 1:floor(Int, n / 4), n)
+x_sol = rand(rng, 1:floor(Int, n/4), n)
 N = sum(x_sol)
 dir = vcat(fill(1, floor(Int, n / 2)), fill(-1, floor(Int, n / 2)), fill(0, mod(n, 2)))
 diffi = x_sol + 0.3 * dir
@@ -178,7 +178,7 @@ diffi = x_sol + 0.3 * dir
 end
 
 n = 20
-x_sol = rand(rng, 1:floor(Int, n / 4), n)
+x_sol = rand(rng, 1:floor(Int, n/4), n)
 diffi = x_sol + 0.3 * rand(rng, [-1, 1], n)
 
 @testset "Unit Simplex LMO" begin
@@ -199,7 +199,7 @@ diffi = x_sol + 0.3 * rand(rng, [-1, 1], n)
 end
 
 n = 20
-x_sol = rand(1:floor(Int, n / 4), n)
+x_sol = rand(1:floor(Int, n/4), n)
 diffi = x_sol + 0.3 * rand([-1, 1], n)
 
 @testset "Reverse Knapsack LMO" begin
@@ -243,7 +243,7 @@ end
     end
 
     # Create L2normBallBLMO
-    blmo = Boscia.L2normBallBLMO()
+    blmo = FrankWolfe.LpNormBallLMO{Float64,2}(1.0)
 
     # Bounds: each variable in [-1, 1] due to L2 ball constraint
     lower_bounds = fill(-2.0, num_int)
@@ -259,72 +259,34 @@ end
     @test isapprox(f(x), f(result[:raw_solution]), atol=1e-6, rtol=1e-3)
 end
 
-@testset "L2normBall BLMO integer 1" begin
+@testset "L2normBall BLMO integer" begin
     n = 20
-
-    # Generate a solution inside the L2 ball (||x|| <= 1)
-    # Round some coordinates to integers for testing
     num_int = 5
-    int_indices = sort(rand(rng, 1:n, num_int))
-    x_sol = zeros(n)
-    x_sol[int_indices[2]] = 1
-
-    function f(x)
-        return 0.5 * sum((x[i] - x_sol[i])^2 for i in eachindex(x))
-    end
-
-    function grad!(storage, x)
-        @. storage = x - x_sol
-    end
 
-    # Create L2normBallBLMO
-    blmo = Boscia.L2normBallBLMO()
+    for sign in (1, -1)
+        int_indices = sort(rand(rng, 1:n, num_int))
 
-    # Bounds: each variable in [-1, 1] due to L2 ball constraint
-    lower_bounds = fill(-2.0, num_int)
-    upper_bounds = fill(3.0, num_int)
+        x_sol = zeros(n)
+        x_sol[int_indices[2]] = sign
 
-    # Some variables are integer
-    int_vars = collect(int_indices)
+        function f(x)
+            return 0.5 * sum((x[i] - x_sol[i])^2 for i in eachindex(x))
+        end
 
-    x, _, result = Boscia.solve(f, grad!, blmo, lower_bounds, upper_bounds, int_vars, n)
+        function grad!(storage, x)
+            @. storage = x - x_sol
+        end
 
-    # Check objective value
-    @test sum(isapprox.(x, x_sol, atol=1e-6, rtol=1e-2)) == n  # Should improve or stay similar
-    @test isapprox(f(x), f(result[:raw_solution]), atol=1e-6, rtol=1e-3)
-end
+        blmo = FrankWolfe.LpNormBallLMO{Float64,2}(1.0)
 
-@testset "L2normBall BLMO integer 2" begin
-    n = 20
+        lower_bounds = fill(-2.0, num_int)
+        upper_bounds = fill(3.0, num_int)
 
-    # Generate a solution inside the L2 ball (||x|| <= 1)
-    # Round some coordinates to integers for testing
-    num_int = 5
-    int_indices = sort(rand(rng, 1:n, num_int))
-    x_sol = zeros(n)
-    x_sol[int_indices[2]] = -1
+        int_vars = collect(int_indices)
 
-    function f(x)
-        return 0.5 * sum((x[i] - x_sol[i])^2 for i in eachindex(x))
-    end
+        x, _, result = Boscia.solve(f, grad!, blmo, lower_bounds, upper_bounds, int_vars, n)
 
-    function grad!(storage, x)
-        @. storage = x - x_sol
+        @test sum(isapprox.(x, x_sol; atol=1e-6, rtol=1e-2)) == n
+        @test isapprox(f(x), f(result[:raw_solution]); atol=1e-6, rtol=1e-3)
     end
-
-    # Create L2normBallBLMO
-    blmo = Boscia.L2normBallBLMO()
-
-    # Bounds: each variable in [-1, 1] due to L2 ball constraint
-    lower_bounds = fill(-2.0, num_int)
-    upper_bounds = fill(3.0, num_int)
-
-    # Some variables are integer
-    int_vars = collect(int_indices)
-
-    x, _, result = Boscia.solve(f, grad!, blmo, lower_bounds, upper_bounds, int_vars, n)
-
-    # Check objective value
-    @test sum(isapprox.(x, x_sol, atol=1e-6, rtol=1e-2)) == n  # Should improve or stay similar
-    @test isapprox(f(x), f(result[:raw_solution]), atol=1e-6, rtol=1e-3)
 end
diff --git a/test/heuristics.jl b/test/heuristics.jl
index 705ba72eb..73ef06087 100644
--- a/test/heuristics.jl
+++ b/test/heuristics.jl
@@ -19,7 +19,7 @@ seed = rand(UInt64)
 rng = StableRNG(seed)
 
 n = 20
-x_sol = rand(rng, 1:floor(Int, n / 4), n)
+x_sol = rand(rng, 1:floor(Int, n/4), n)
 N = sum(x_sol)
 dir = vcat(fill(1, floor(Int, n / 2)), fill(-1, floor(Int, n / 2)), fill(0, mod(n, 2)))
 diffi = x_sol + 0.3 * dir
@@ -52,7 +52,7 @@ diffi = x_sol + 0.3 * dir
 end
 
 n = 20
-x_sol = rand(rng, 1:floor(Int, n / 4), n)
+x_sol = rand(rng, 1:floor(Int, n/4), n)
 diffi = x_sol + 0.3 * rand(rng, [-1, 1], n)
 
 @testset "Hyperplane Aware Rounding - Unit Simplex" begin

From 206c9b442ccc2497d0b1b3de5ad3749300a0e843 Mon Sep 17 00:00:00 2001
From: LJS42 <1457928939@qq.com>
Date: Wed, 17 Dec 2025 16:54:32 +0100
Subject: [PATCH 14/20] format

---
 test/LMO_test.jl   | 6 +++---
 test/heuristics.jl | 4 ++--
 2 files changed, 5 insertions(+), 5 deletions(-)

diff --git a/test/LMO_test.jl b/test/LMO_test.jl
index 2bbe7dc63..0b3851171 100644
--- a/test/LMO_test.jl
+++ b/test/LMO_test.jl
@@ -140,7 +140,7 @@ end
 end
 
 n = 20
-x_sol = rand(rng, 1:floor(Int, n/4), n)
+x_sol = rand(rng, 1:floor(Int, n / 4), n)
 N = sum(x_sol)
 dir = vcat(fill(1, floor(Int, n / 2)), fill(-1, floor(Int, n / 2)), fill(0, mod(n, 2)))
 diffi = x_sol + 0.3 * dir
@@ -178,7 +178,7 @@ diffi = x_sol + 0.3 * dir
 end
 
 n = 20
-x_sol = rand(rng, 1:floor(Int, n/4), n)
+x_sol = rand(rng, 1:floor(Int, n / 4), n)
 diffi = x_sol + 0.3 * rand(rng, [-1, 1], n)
 
 @testset "Unit Simplex LMO" begin
@@ -199,7 +199,7 @@ diffi = x_sol + 0.3 * rand(rng, [-1, 1], n)
 end
 
 n = 20
-x_sol = rand(1:floor(Int, n/4), n)
+x_sol = rand(1:floor(Int, n / 4), n)
 diffi = x_sol + 0.3 * rand([-1, 1], n)
 
 @testset "Reverse Knapsack LMO" begin
diff --git a/test/heuristics.jl b/test/heuristics.jl
index 73ef06087..705ba72eb 100644
--- a/test/heuristics.jl
+++ b/test/heuristics.jl
@@ -19,7 +19,7 @@ seed = rand(UInt64)
 rng = StableRNG(seed)
 
 n = 20
-x_sol = rand(rng, 1:floor(Int, n/4), n)
+x_sol = rand(rng, 1:floor(Int, n / 4), n)
 N = sum(x_sol)
 dir = vcat(fill(1, floor(Int, n / 2)), fill(-1, floor(Int, n / 2)), fill(0, mod(n, 2)))
 diffi = x_sol + 0.3 * dir
@@ -52,7 +52,7 @@ diffi = x_sol + 0.3 * dir
 end
 
 n = 20
-x_sol = rand(rng, 1:floor(Int, n/4), n)
+x_sol = rand(rng, 1:floor(Int, n / 4), n)
 diffi = x_sol + 0.3 * rand(rng, [-1, 1], n)
 
 @testset "Hyperplane Aware Rounding - Unit Simplex" begin

From 1509784f0687dba06066a1aa9335b294c134a871 Mon Sep 17 00:00:00 2001
From: LJS42 <1457928939@qq.com>
Date: Wed, 14 Jan 2026 10:24:53 +0100
Subject: [PATCH 15/20] doc

---
 src/polytope_blmos.jl |  5 +++--
 test/LMO_test.jl      | 15 ++++-----------
 2 files changed, 7 insertions(+), 13 deletions(-)

diff --git a/src/polytope_blmos.jl b/src/polytope_blmos.jl
index b91298553..d8a45b09b 100644
--- a/src/polytope_blmos.jl
+++ b/src/polytope_blmos.jl
@@ -718,10 +718,11 @@ function rounding_hyperplane_heuristic(
 end
 
 """
-    2normBallBLMO()
+    FrankWolfe.LpNormBallLMO{T,2}
 
-BLMO denotes the L2normBall, It is unit ball which means R = 1
+BLMO denotes the only unit L2normBall, It is unit ball which means R = 1
 """
+const _L2norm_BALL_DOC = nothing
 
 function bounded_compute_extreme_point(
     lmo::FrankWolfe.LpNormBallLMO{T,2},
diff --git a/test/LMO_test.jl b/test/LMO_test.jl
index 0b3851171..fb73c871c 100644
--- a/test/LMO_test.jl
+++ b/test/LMO_test.jl
@@ -227,13 +227,6 @@ end
     x_sol = randn(rng, n)
     x_sol = x_sol ./ (norm(x_sol) * 1.5)
 
-    # Round some coordinates to integers for testing
-    num_int = 5
-    int_indices = sort(rand(rng, 1:n, num_int))
-    for idx in int_indices
-        x_sol[idx] = 0
-    end
-
     function f(x)
         return 0.5 * sum((x[i] - x_sol[i])^2 for i in eachindex(x))
     end
@@ -246,8 +239,8 @@ end
     blmo = FrankWolfe.LpNormBallLMO{Float64,2}(1.0)
 
     # Bounds: each variable in [-1, 1] due to L2 ball constraint
-    lower_bounds = fill(-2.0, num_int)
-    upper_bounds = fill(2.0, num_int)
+    lower_bounds = fill(-1.0, num_int)
+    upper_bounds = fill(1.0, num_int)
 
     # Some variables are integer
     int_vars = collect(int_indices)
@@ -279,8 +272,8 @@ end
 
         blmo = FrankWolfe.LpNormBallLMO{Float64,2}(1.0)
 
-        lower_bounds = fill(-2.0, num_int)
-        upper_bounds = fill(3.0, num_int)
+        lower_bounds = fill(-1.0, num_int)
+        upper_bounds = fill(1.0, num_int)
 
         int_vars = collect(int_indices)
 

From 947f088063487d2e0d3cea9517717997e712f03d Mon Sep 17 00:00:00 2001
From: LJS42 <1457928939@qq.com>
Date: Wed, 14 Jan 2026 10:53:19 +0100
Subject: [PATCH 16/20] document

---
 src/polytope_blmos.jl | 11 ++++++++---
 1 file changed, 8 insertions(+), 3 deletions(-)

diff --git a/src/polytope_blmos.jl b/src/polytope_blmos.jl
index d8a45b09b..9bdca495a 100644
--- a/src/polytope_blmos.jl
+++ b/src/polytope_blmos.jl
@@ -718,11 +718,16 @@ function rounding_hyperplane_heuristic(
 end
 
 """
-    FrankWolfe.LpNormBallLMO{T,2}
+# FrankWolfe.LpNormBallLMO{T,2}
 
-BLMO denotes the only unit L2normBall, It is unit ball which means R = 1
+Unit L2-norm ball linear minimization oracle (BLMO).
+
+This type is defined in the `FrankWolfe` package.  
+It solves the linear minimization problem over the Euclidean unit ball.
+
+You can use it with the helpers in this package for integer-constrained Frank-Wolfe optimization.
 """
-const _L2norm_BALL_DOC = nothing
+const L2BallLMO = FrankWolfe.LpNormBallLMO{Float64,2} 
 
 function bounded_compute_extreme_point(
     lmo::FrankWolfe.LpNormBallLMO{T,2},

From d3f9cf2fea095dd9f1e24d4a8dcb925ce630de51 Mon Sep 17 00:00:00 2001
From: LJS42 <1457928939@qq.com>
Date: Wed, 14 Jan 2026 10:55:28 +0100
Subject: [PATCH 17/20] int

---
 test/LMO_test.jl | 7 +++++++
 1 file changed, 7 insertions(+)

diff --git a/test/LMO_test.jl b/test/LMO_test.jl
index fb73c871c..61aa00668 100644
--- a/test/LMO_test.jl
+++ b/test/LMO_test.jl
@@ -227,6 +227,13 @@ end
     x_sol = randn(rng, n)
     x_sol = x_sol ./ (norm(x_sol) * 1.5)
 
+    # Round some coordinates to integers for testing
+    num_int = 3
+    int_indices = sort(rand(rng, 1:n, num_int))
+    for idx in int_indices
+        x_sol[idx] = 0
+    end
+
     function f(x)
         return 0.5 * sum((x[i] - x_sol[i])^2 for i in eachindex(x))
     end

From 1700eb3999cef8b27580573c52554f31b9efc3e0 Mon Sep 17 00:00:00 2001
From: LJS42 <1457928939@qq.com>
Date: Wed, 14 Jan 2026 10:57:04 +0100
Subject: [PATCH 18/20] format

---
 src/polytope_blmos.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/polytope_blmos.jl b/src/polytope_blmos.jl
index 9bdca495a..9ff6f908b 100644
--- a/src/polytope_blmos.jl
+++ b/src/polytope_blmos.jl
@@ -727,7 +727,7 @@ It solves the linear minimization problem over the Euclidean unit ball.
 
 You can use it with the helpers in this package for integer-constrained Frank-Wolfe optimization.
 """
-const L2BallLMO = FrankWolfe.LpNormBallLMO{Float64,2} 
+const L2BallLMO = FrankWolfe.LpNormBallLMO{Float64,2}
 
 function bounded_compute_extreme_point(
     lmo::FrankWolfe.LpNormBallLMO{T,2},

From 50db5993292a207fca91c8c96ae34da9d1a9f10c Mon Sep 17 00:00:00 2001
From: LJS42 <1457928939@qq.com>
Date: Wed, 21 Jan 2026 14:29:53 +0100
Subject: [PATCH 19/20] document and rhs

---
 src/polytope_blmos.jl | 54 +++++++++++++++++++++++++++++++++++++++++--
 test/LMO_test.jl      | 35 ++++++++++++++++++++++++++++
 2 files changed, 87 insertions(+), 2 deletions(-)

diff --git a/src/polytope_blmos.jl b/src/polytope_blmos.jl
index 9ff6f908b..7dfce28aa 100644
--- a/src/polytope_blmos.jl
+++ b/src/polytope_blmos.jl
@@ -726,9 +726,19 @@ This type is defined in the `FrankWolfe` package.
 It solves the linear minimization problem over the Euclidean unit ball.
 
 You can use it with the helpers in this package for integer-constrained Frank-Wolfe optimization.
-"""
-const L2BallLMO = FrankWolfe.LpNormBallLMO{Float64,2}
 
+# bounded_compute_extreme_point
+
+Compute an extreme point of an L2-norm unit ball under integer constraints.
+
+The function returns a vector `v` that maximizes `dot(v, d)` over the L2 unit ball,
+while respecting integer variable bounds (`lb`, `ub`).  
+
+- For integer variables, if bounds force ±1, return that vertex immediately.
+- Otherwise, choose the integer index with largest |d[i]| and assign ±1 if it improves the objective.
+- Continuous variables follow the negative gradient direction and are normalized to unit L2-norm.
+- Integer entries not assigned ±1 are set to zero in the continuous hyperplane.
+"""
 function bounded_compute_extreme_point(
     lmo::FrankWolfe.LpNormBallLMO{T,2},
     d,
@@ -786,6 +796,11 @@ function is_simple_linear_feasible(lmo::FrankWolfe.LpNormBallLMO{T,2}, v) where
 end
 
 function check_feasibility(lmo::FrankWolfe.LpNormBallLMO{T,2}, lb, ub, int_vars, n) where {T}
+    if lmo.right_hand_side != 1
+        println("right_hand_side > 1: cannot handle this case now")
+        return INFEASIBLE 
+    end
+
     if any(lb .> 1) || any(ub .< -1)
         return INFEASIBLE
     end
@@ -801,3 +816,38 @@ function check_feasibility(lmo::FrankWolfe.LpNormBallLMO{T,2}, lb, ub, int_vars,
     end
     return OPTIMAL
 end
+
+@testset "L2normBall BLMO rhs > 1 triggers infeasible" begin
+    n = 10
+    num_int = 3
+
+    rng = Random.default_rng()
+    int_indices = sort(rand(rng, 1:n, num_int))
+
+    x_sol = zeros(n)
+    x_sol[int_indices[1]] = 1.0
+
+    function f(x)
+        return 0.5 * sum((x[i] - x_sol[i])^2 for i in eachindex(x))
+    end
+
+    function grad!(storage, x)
+        @. storage = x - x_sol
+    end
+
+    blmo = FrankWolfe.LpNormBallLMO{Float64,2}(1.5)  
+
+    lower_bounds = fill(-1.0, num_int)
+    upper_bounds = fill(1.0, num_int)
+    int_vars = collect(int_indices)
+
+    result = try
+        x, _, sol = Boscia.solve(f, grad!, blmo, lower_bounds, upper_bounds, int_vars, n)
+        false  
+    catch e
+        println("Caught error as expected: ", e)
+        true   
+    end
+
+    @test result == true
+end
diff --git a/test/LMO_test.jl b/test/LMO_test.jl
index 61aa00668..88b933d9c 100644
--- a/test/LMO_test.jl
+++ b/test/LMO_test.jl
@@ -290,3 +290,38 @@ end
         @test isapprox(f(x), f(result[:raw_solution]); atol=1e-6, rtol=1e-3)
     end
 end
+
+@testset "L2normBall BLMO rhs > 1 triggers infeasible" begin
+    n = 10
+    num_int = 3
+
+    rng = Random.default_rng()
+    int_indices = sort(rand(rng, 1:n, num_int))
+
+    x_sol = zeros(n)
+    x_sol[int_indices[1]] = 1.0
+
+    function f(x)
+        return 0.5 * sum((x[i] - x_sol[i])^2 for i in eachindex(x))
+    end
+
+    function grad!(storage, x)
+        @. storage = x - x_sol
+    end
+
+    blmo = FrankWolfe.LpNormBallLMO{Float64,2}(1.5)  
+
+    lower_bounds = fill(-1.0, num_int)
+    upper_bounds = fill(1.0, num_int)
+    int_vars = collect(int_indices)
+
+    result = try
+        x, _, sol = Boscia.solve(f, grad!, blmo, lower_bounds, upper_bounds, int_vars, n)
+        false 
+    catch e
+        println("Caught error as expected: ", e)
+        true   
+    end
+
+    @test result == true
+end

From ce497c6e464fc00037ccd62024cb910b197bacb8 Mon Sep 17 00:00:00 2001
From: LJS42 <1457928939@qq.com>
Date: Wed, 21 Jan 2026 14:31:04 +0100
Subject: [PATCH 20/20] format

---
 src/polytope_blmos.jl | 37 +------------------------------------
 test/LMO_test.jl      |  6 +++---
 2 files changed, 4 insertions(+), 39 deletions(-)

diff --git a/src/polytope_blmos.jl b/src/polytope_blmos.jl
index 7dfce28aa..9c5ba2c56 100644
--- a/src/polytope_blmos.jl
+++ b/src/polytope_blmos.jl
@@ -798,7 +798,7 @@ end
 function check_feasibility(lmo::FrankWolfe.LpNormBallLMO{T,2}, lb, ub, int_vars, n) where {T}
     if lmo.right_hand_side != 1
         println("right_hand_side > 1: cannot handle this case now")
-        return INFEASIBLE 
+        return INFEASIBLE
     end
 
     if any(lb .> 1) || any(ub .< -1)
@@ -816,38 +816,3 @@ function check_feasibility(lmo::FrankWolfe.LpNormBallLMO{T,2}, lb, ub, int_vars,
     end
     return OPTIMAL
 end
-
-@testset "L2normBall BLMO rhs > 1 triggers infeasible" begin
-    n = 10
-    num_int = 3
-
-    rng = Random.default_rng()
-    int_indices = sort(rand(rng, 1:n, num_int))
-
-    x_sol = zeros(n)
-    x_sol[int_indices[1]] = 1.0
-
-    function f(x)
-        return 0.5 * sum((x[i] - x_sol[i])^2 for i in eachindex(x))
-    end
-
-    function grad!(storage, x)
-        @. storage = x - x_sol
-    end
-
-    blmo = FrankWolfe.LpNormBallLMO{Float64,2}(1.5)  
-
-    lower_bounds = fill(-1.0, num_int)
-    upper_bounds = fill(1.0, num_int)
-    int_vars = collect(int_indices)
-
-    result = try
-        x, _, sol = Boscia.solve(f, grad!, blmo, lower_bounds, upper_bounds, int_vars, n)
-        false  
-    catch e
-        println("Caught error as expected: ", e)
-        true   
-    end
-
-    @test result == true
-end
diff --git a/test/LMO_test.jl b/test/LMO_test.jl
index 88b933d9c..c5e4cdcad 100644
--- a/test/LMO_test.jl
+++ b/test/LMO_test.jl
@@ -309,7 +309,7 @@ end
         @. storage = x - x_sol
     end
 
-    blmo = FrankWolfe.LpNormBallLMO{Float64,2}(1.5)  
+    blmo = FrankWolfe.LpNormBallLMO{Float64,2}(1.5)
 
     lower_bounds = fill(-1.0, num_int)
     upper_bounds = fill(1.0, num_int)
@@ -317,10 +317,10 @@ end
 
     result = try
         x, _, sol = Boscia.solve(f, grad!, blmo, lower_bounds, upper_bounds, int_vars, n)
-        false 
+        false
     catch e
         println("Caught error as expected: ", e)
-        true   
+        true
     end
 
     @test result == true