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/src/polytope_blmos.jl b/src/polytope_blmos.jl
index 1656ddb21..9c5ba2c56 100644
--- a/src/polytope_blmos.jl
+++ b/src/polytope_blmos.jl
@@ -716,3 +716,103 @@ function rounding_hyperplane_heuristic(
     end
     return [z], false
 end
+
+"""
+# FrankWolfe.LpNormBallLMO{T,2}
+
+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.
+
+# 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,
+    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
+            v = zeros(length(d))
+            v[idx] = 1
+            return v
+        elseif ub[i] < 0
+            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
+    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))
+        if d[max_idx] < 0
+            v[max_idx] = 1
+        elseif d[max_idx] > 0
+            v[max_idx] = -1
+        end
+    end
+    return v
+end
+
+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
+    end
+    return true
+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
+    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
+        end
+        if count > 1
+            return INFEASIBLE
+        end
+    end
+    return OPTIMAL
+end
diff --git a/test/LMO_test.jl b/test/LMO_test.jl
index f44f4415c..c5e4cdcad 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,110 @@ 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 continuous" 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 = 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
+
+    function grad!(storage, x)
+        @. storage = x - x_sol
+    end
+
+    # Create L2normBallBLMO
+    blmo = FrankWolfe.LpNormBallLMO{Float64,2}(1.0)
+
+    # Bounds: each variable in [-1, 1] due to L2 ball constraint
+    lower_bounds = fill(-1.0, num_int)
+    upper_bounds = fill(1.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
+
+@testset "L2normBall BLMO integer" begin
+    n = 20
+    num_int = 5
+
+    for sign in (1, -1)
+        int_indices = sort(rand(rng, 1:n, num_int))
+
+        x_sol = zeros(n)
+        x_sol[int_indices[2]] = sign
+
+        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.0)
+
+        lower_bounds = fill(-1.0, num_int)
+        upper_bounds = fill(1.0, num_int)
+
+        int_vars = collect(int_indices)
+
+        x, _, result = Boscia.solve(f, grad!, blmo, lower_bounds, upper_bounds, int_vars, 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
+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/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 5c4af8879..4681970e2 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