diff --git a/Project.toml b/Project.toml
index 39ee062aa..162f2b66a 100644
--- a/Project.toml
+++ b/Project.toml
@@ -46,6 +46,7 @@ Statistics = "1.6"
 Suppressor = "0.2"
 Test = "1"
 julia = "1.6"
+JuMP = "1"
 
 [extras]
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
@@ -60,6 +61,7 @@ SCIP = "82193955-e24f-5292-bf16-6f2c5261a85f"
 StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
 Suppressor = "fd094767-a336-5f1f-9728-57cf17d0bbfb"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
 
 [targets]
-test = ["Test", "HiGHS", "Distributions", "SCIP", "DoubleFloats", "StableRNGs", "Aqua", "Suppressor", "Graphs", "Bonobo", "CSV", "CombinatorialLinearOracles"]
+test = ["Test", "HiGHS", "Distributions", "SCIP", "DoubleFloats", "StableRNGs", "Aqua", "Suppressor", "Graphs", "Bonobo", "CSV", "CombinatorialLinearOracles", "JuMP"]
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 0a3766c8f..93db85397 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 7914cf096..3e9425f86 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 67a25e209..0b4fb6ecf 100644
--- a/src/polytope_blmos.jl
+++ b/src/polytope_blmos.jl
@@ -1025,3 +1025,248 @@ end
 _compute_k_hypersimplex(lmo::FrankWolfe.HyperSimplexLMO, direction) = lmo.K
 _compute_k_hypersimplex(lmo::FrankWolfe.UnitHyperSimplexLMO, direction) =
     min(lmo.K, length(direction), sum(<(0), direction))
+
+
+
+"""
+    bounded_compute_extreme_point(lmo::KSparseLMO, direction, lb, ub, int_vars)
+
+    KSparse: C = B_1(τK) ∩ B_∞(τ)
+
+Compute an extreme point of a K-sparse LMO along a given direction.  
+
+Compute an extreme point of a K-sparse LMO using a greedy strategy:  
+integer variables are fixed first, then the remaining entries are filled greedily according to the largest absolute values in `direction`, respecting the K-sparse and bounds constraints.
+"""
+function bounded_compute_extreme_point(
+    lmo::FrankWolfe.KSparseLMO{T},
+    direction,
+    lb,
+    ub,
+    int_vars;
+    atol=1e-6,
+    rtol=1e-4,
+    v=nothing,
+    kwargs...,
+) where {T}
+    K = lmo.K
+    K_indices = sortperm(direction, by=abs, rev=true)
+    v = spzeros(Float64, length(direction))
+
+    for i in int_vars
+        idx = findfirst(==(i), int_vars)
+        if lb[idx] > 0
+            v[i] = ceil(lb[idx])
+        elseif ub[idx] < 0
+            v[i] = floor(ub[idx])
+        end
+    end
+
+    for i in K_indices
+
+        rem = K * lmo.right_hand_side - sum(abs, v)
+
+        if isapprox(rem, 0.0; atol=1e-6, rtol=1e-8)
+            break
+        end
+
+        if i in int_vars
+            idx = findfirst(==(i), int_vars)
+            lower_eff = clamp(ceil(max(-lmo.right_hand_side, -rem)), ceil(lb[idx]), floor(ub[idx]))
+            upper_eff = clamp(floor(min(lmo.right_hand_side, rem)), ceil(lb[idx]), floor(ub[idx]))
+        else
+            lower_eff = max(-rem, -lmo.right_hand_side)
+            upper_eff = min(rem, lmo.right_hand_side)
+        end
+        if direction[i] > 0
+            v[i] = lower_eff
+        elseif direction[i] < 0
+            v[i] = upper_eff
+        else
+            v[i] = clamp(0, lower_eff, upper_eff)
+        end
+    end
+    return v
+end
+
+"""
+    bounded_compute_inface_extreme_point(lmo::KSparseLMO, direction, x, lb, ub, int_vars)
+
+    KSparse: C = B_1(τK) ∩ B_∞(τ)
+
+This function computes an extreme point of a K-sparse LMO using a greedy strategy:  
+integer variables are fixed first, then the remaining entries are filled greedily according to the largest absolute values in `direction`, respecting the K-sparse and bounds constraints.
+This routine is similar to `compute_bounded_extreme_point`, but here we additionally fix the variable entries to the current minimal face of the iterate.
+"""
+function bounded_compute_inface_extreme_point(
+    sblmo::FrankWolfe.KSparseLMO{T},
+    direction,
+    x,
+    lb,
+    ub,
+    int_vars;
+    atol=1e-6,
+    rtol=1e-4,
+    kwargs...,
+) where {T}
+    n = length(direction)
+    rhs = sblmo.right_hand_side
+    K = sblmo.K
+
+    v = copy(x)
+    fixed_vars = Int[]
+
+    # Identify fixed coordinates (already on the face boundary)
+    for i in 1:n
+        if i in int_vars
+            idx = findfirst(==(i), int_vars)
+            if isapprox(x[i], ceil(lb[idx]); atol=atol, rtol=rtol) ||
+               isapprox(x[i], floor(ub[idx]); atol=atol, rtol=rtol)
+                push!(fixed_vars, i)
+            end
+        else
+            if isapprox(x[i], rhs; atol=atol, rtol=rtol) ||
+               isapprox(x[i], -rhs; atol=atol, rtol=rtol)
+                push!(fixed_vars, i)
+            end
+        end
+    end
+
+    free_idx = setdiff(1:n, fixed_vars)
+    if isempty(free_idx)
+        return v  # already at in-face extreme point
+    end
+
+    # Construct new in-face extreme point
+    v[free_idx] .= 0
+    free_idx = sort(free_idx, by=i -> abs(direction[i]), rev=true)
+
+    for i in intersect(int_vars, free_idx)
+        idx = findfirst(==(i), int_vars)
+        if lb[idx] > 0
+            v[i] = ceil(lb[idx])
+        elseif ub[idx] < 0
+            v[i] = floor(ub[idx])
+        end
+    end
+
+    for i in free_idx
+
+        rem = K * sblmo.right_hand_side - sum(abs, v)
+
+        if isapprox(rem, 0.0; atol=1e-6, rtol=1e-8)
+            break
+        end
+
+        if i in int_vars
+            idx = findfirst(==(i), int_vars)
+            lower_eff =
+                clamp(ceil(max(-sblmo.right_hand_side, -rem)), ceil(lb[idx]), floor(ub[idx]))
+            upper_eff = clamp(floor(min(sblmo.right_hand_side, rem)), ceil(lb[idx]), floor(ub[idx]))
+        else
+            lower_eff = max(-rem, -sblmo.right_hand_side)
+            upper_eff = min(rem, sblmo.right_hand_side)
+        end
+        if direction[i] > 0
+            v[i] = lower_eff
+        elseif direction[i] < 0
+            v[i] = upper_eff
+        else
+            v[i] = clamp(0, lower_eff, upper_eff)
+        end
+    end
+    return v
+end
+
+"""
+    bounded_dicg_maximum_step(sblmo::KSparseLMO, direction, x, lb, ub, int_vars; tol=1e-6)
+
+Compute the maximum feasible step size `γ` along `direction` from point `x`  
+for the K-sparse LMO, respecting integer and continuous bounds.  
+The step is limited by hitting either the bounds or the K-sparse LMO constraints.
+"""
+function bounded_dicg_maximum_step(
+    sblmo::FrankWolfe.KSparseLMO{T},
+    direction,
+    x,
+    lb,
+    ub,
+    int_vars;
+    tol=1e-6,
+    kwargs...,
+) where {T}
+    gamma_max = one(eltype(direction))
+    for idx in eachindex(x)
+        di = direction[idx]
+        if di > tol
+            if idx in int_vars
+                int_idx = findfirst(==(idx), int_vars)
+                gamma_max = min(gamma_max, (x[idx] - lb[int_idx]) / di)
+            else
+                gamma_max = min(gamma_max, (x[idx] - (-sblmo.rhs)) / di)
+            end
+        elseif di < -tol
+            if idx in int_vars
+                int_idx = findfirst(==(idx), int_vars)
+                gamma_max = min(gamma_max, (ub[int_idx] - x[idx]) / -di)
+            else
+                gamma_max = min(gamma_max, (sblmo.rhs - x[idx]) / -di)
+            end
+        end
+
+        if isapprox(gamma_max, 0.0; atol=tol)
+            return 0.0
+        end
+    end
+
+    return max(gamma_max, 0.0)
+end
+
+"""
+    is_simple_linear_feasible(sblmo::KSparseLMO, x, int_vars; tol=1e-8)
+
+If L1 norm is larger than K * τ or Infinite norm is larger than τ, than the point is not feasible.
+"""
+function is_simple_linear_feasible(
+    sblmo::FrankWolfe.KSparseLMO{T},
+    v;
+    int_vars=Int[],
+    tol=1e-8,
+) where {T}
+    τ = sblmo.right_hand_side
+    K = sblmo.K
+
+    l1_norm = sum(abs, v)
+    if l1_norm > K * τ + tol
+        @debug "v violates sparsity constraint: ‖v‖₀ = $l1_norm > K*τ = $(K*τ)"
+        return false
+    end
+
+    if any(abs.(v) .> τ + tol)
+        @debug "v violates bound constraint: some |vᵢ| > τ = $τ (v=$v)"
+        return false
+    end
+
+    return true
+end
+
+"""
+    check_feasibility(sblmo::KSparseLMO, lb, ub, int_vars, n; tol=1e-8)
+
+    If the absolute value of lb or ub is already larger than τ, or the sum of each dimensions is larger than K * τ then there is no feasible point in this region.
+"""
+function check_feasibility(sblmo::FrankWolfe.KSparseLMO{T}, lb, ub, int_vars, n; tol=1e-8) where {T}
+    K = sblmo.K
+    τ = sblmo.right_hand_side
+    n_int = length(lb)
+
+    for i in 1:n_int
+        v = clamp(τ, lb[i], ub[i])
+        if abs.(v) .> τ + tol
+            return INFEASIBLE
+        end
+    end
+    return OPTIMAL
+end
+
+is_decomposition_invariant_oracle_simple(::FrankWolfe.KSparseLMO{T}) where {T} = true
diff --git a/test/LMO_test.jl b/test/LMO_test.jl
index 955e553bd..33a295be3 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
@@ -288,3 +289,43 @@ 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
+
+n = 20
+n_int = 10
+sparsity = 0.3
+x_sol = [rand() < sparsity ? 0 : rand(1:floor(Int, n / 4)) for _ in 1:n]
+
+@testset "KSparse LMO" begin
+    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
+
+    K = count(!iszero, x_sol)
+    τ = norm(x_sol, Inf)
+
+    sblmo = FrankWolfe.KSparseLMO{Float64}(K, τ)
+
+    x, _, result =
+        Boscia.solve(f, grad!, sblmo, fill(0.0, n_int), fill(100.0, n_int), collect(1:n_int), n)
+
+    settings = Boscia.create_default_settings()
+    settings.frank_wolfe[:variant] = Boscia.DecompositionInvariantConditionalGradient()
+    x_dicg, _, result_dicg = Boscia.solve(
+        f,
+        grad!,
+        sblmo,
+        fill(0.0, n),
+        fill(100.0, n),
+        collect(1:n),
+        n,
+        settings=settings,
+    )
+
+    @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)
+    @test sum(isapprox.(x_dicg, x_sol, atol=1e-6, rtol=1e-2)) == n
+    @test isapprox(f(x_dicg), f(result[:raw_solution]), atol=1e-6, rtol=1e-3)
+end
diff --git a/test/custom_lmo_testing.jl b/test/custom_lmo_testing.jl
new file mode 100644
index 000000000..03235ad90
--- /dev/null
+++ b/test/custom_lmo_testing.jl
@@ -0,0 +1,110 @@
+using MathOptInterface
+using Test
+using HiGHS
+using JuMP
+using Boscia
+using FrankWolfe
+
+
+@testset "LMO vs MOI on the same feasible region - K" begin
+    n = 10
+    K = 5
+    rhs = 4.0
+
+    direction = randn(n)
+    lb = fill(-7.0, 3)
+    ub = fill(-2.0, 3)
+    int_vars = [2, 5, 7]
+
+    # --- LMO ---
+    lmo = FrankWolfe.KSparseLMO(K, rhs)
+
+    v_lmo = Boscia.bounded_compute_extreme_point(lmo, direction, lb, ub, int_vars)
+
+    # --- MOI ---
+    model = Model(HiGHS.Optimizer)
+
+    # variables v with bounds
+    @variable(model, v[i=1:n])
+
+    # integer variables
+    for i in int_vars
+        set_integer(v[i])
+    end
+
+    @variable(model, abs_v[i=1:n] >= 0)
+
+    for i in 1:n
+        if i in int_vars
+            idx = findfirst(==(i), int_vars)
+            @constraint(model, lb[idx] <= v[i] <= ub[idx])
+        end
+        @constraint(model, abs_v[i] <= rhs)
+        @constraint(model, abs_v[i] >= v[i])
+        @constraint(model, abs_v[i] >= -v[i])
+    end
+
+    @constraint(model, sum(abs_v[i] for i in 1:n) <= K * rhs)
+
+    # objective
+    @objective(model, Min, sum(direction[i] * v[i] for i in 1:n))
+
+    optimize!(model)
+
+    v_moi = value.(v)
+
+    # compare objective values (not pointwise equality!)
+    @test isapprox(v_lmo, v_moi; atol=1e-6)
+end
+
+@testset "LMO on inface vs MOI on the same feasible region - K" begin
+    n = 10
+    K = 5
+    rhs = 3.0
+
+    direction = randn(n)
+    lb = fill(1.0, 3)
+    ub = fill(7.0, 3)
+    int_vars = [2, 5, 7]
+    x = zeros(n)
+
+    # --- LMO ---
+    lmo = FrankWolfe.KSparseLMO(K, rhs)
+
+    v_lmo = Boscia.bounded_compute_inface_extreme_point(lmo, direction, x, lb, ub, int_vars)
+
+    # --- MOI ---
+    model = Model(HiGHS.Optimizer)
+
+    # variables v with bounds
+    @variable(model, v[i=1:n])
+
+    # integer variables
+    for i in int_vars
+        set_integer(v[i])
+    end
+
+    @variable(model, abs_v[i=1:n] >= 0)
+
+    for i in 1:n
+        if i in int_vars
+            idx = findfirst(==(i), int_vars)
+            @constraint(model, lb[idx] <= v[i] <= ub[idx])
+        end
+        @constraint(model, abs_v[i] <= rhs)
+        @constraint(model, abs_v[i] >= v[i])
+        @constraint(model, abs_v[i] >= -v[i])
+    end
+
+    @constraint(model, sum(abs_v[i] for i in 1:n) <= K * rhs)
+
+    # objective
+    @objective(model, Min, sum(direction[i] * v[i] for i in 1:n))
+
+    optimize!(model)
+
+    v_moi = value.(v)
+
+    # compare objective values (not pointwise equality!)
+    @test isapprox(v_lmo, v_moi; atol=1e-6)
+end
diff --git a/test/heuristics.jl b/test/heuristics.jl
index 0d9774ceb..52363bbcd 100644
--- a/test/heuristics.jl
+++ b/test/heuristics.jl
@@ -8,6 +8,7 @@ import Bonobo
 using HiGHS
 using Printf
 using Dates
+using LinearAlgebra
 const MOI = MathOptInterface
 const MOIU = MOI.Utilities
 using StableRNGs
diff --git a/test/indicator_test.jl b/test/indicator_test.jl
index a9e8d3bea..dac0e6bbb 100644
--- a/test/indicator_test.jl
+++ b/test/indicator_test.jl
@@ -57,7 +57,7 @@ rng = StableRNG(seed)
     @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/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