diff --git a/examples/birkhoff.jl b/examples/birkhoff.jl
index ecea58e1c..41ee6fac8 100644
--- a/examples/birkhoff.jl
+++ b/examples/birkhoff.jl
@@ -27,8 +27,10 @@ import HiGHS
 # The variables are ordered (Y, X, theta) in the MOI model
 # the objective only uses the last n^2 variables
 # Small dimensions since the size of the problem grows quickly (2 k n^2 + k variables)
-n = 3
-k = 2
+
+Random.seed!(5)
+n = 4
+k = 3
 
 # generate random doubly stochastic matrix
 const Xstar = rand(n, n)
@@ -94,7 +96,7 @@ function build_birkhoff_lmo()
 end
 
 lmo = build_birkhoff_lmo()
-x, _, _ = Boscia.solve(f, grad!, lmo, verbose=true)
+x, _, _ = Boscia.solve(f, grad!, lmo, verbose=true, lazy=false, variant=Boscia.DICG())
 
 
 # TODO the below needs to be fixed
diff --git a/examples/co_localization.jl b/examples/co_localization.jl
new file mode 100644
index 000000000..5179a5b9a
--- /dev/null
+++ b/examples/co_localization.jl
@@ -0,0 +1,272 @@
+using Boscia
+using FrankWolfe
+using Test
+using Random
+using FiniteDifferences
+using SCIP
+using LinearAlgebra
+using DataFrames
+using CSV
+using MathOptInterface: MathOptInterface
+const MOI = MathOptInterface
+using HiGHS: HiGHS
+
+
+
+# For bug hunting:
+seed = rand(UInt64)
+@show seed
+Random.seed!(seed)
+
+include("co_localization_BLMO.jl")
+
+
+
+
+
+function build_objective()
+
+	param_mu = 0.4
+	param_lambda = 0.1
+
+	linVector = CSV.read("./data/qp_linVector.csv", DataFrame; header = false) |> Matrix
+	lapMatrix = CSV.read("./data/qp_lapMatrix.csv", DataFrame; header = false) |> Matrix
+	rrMatrix = CSV.read("./data/qp_rrMatrix.csv", DataFrame; header = false) |> Matrix
+	A = lapMatrix + param_mu * rrMatrix
+	@show size(A)
+	b = param_lambda * linVector
+
+	function f(x)
+
+		return (1 / 2) * (transpose(x) * A * x) + sum(transpose(b) * x)
+	end
+
+	function grad!(storage, x)
+		@. storage = $(A * x) + b
+	end
+
+	return f, grad!
+end
+
+
+
+
+
+
+function run_args_boscia(mode, f, grad!, lmo, verbose, time_limit; variant = Boscia.BPCG(), fw_epsilon = 1e-2, dual_gap = 1e-6, rel_dual_gap = 1e-2, min_node_fw_epsilon = 1e-6)
+	if contains(mode, "lazy")
+		if contains(mode, "ws")
+			if contains(mode, "full")
+				x, _, result =
+					Boscia.solve(
+						f,
+						grad!,
+						lmo,
+						verbose = verbose,
+						time_limit = time_limit,
+						lazy = true,
+						variant = variant,
+						use_DICG_warm_start = true,
+						use_full_ws = true,
+						fw_epsilon = fw_epsilon,
+						dual_gap = dual_gap,
+						rel_dual_gap = rel_dual_gap,
+						min_node_fw_epsilon = min_node_fw_epsilon,
+					)
+			else
+				x, _, result =
+					Boscia.solve(
+						f,
+						grad!,
+						lmo,
+						verbose = verbose,
+						time_limit = time_limit,
+						lazy = true,
+						variant = variant,
+						use_DICG_warm_start = true,
+						use_full_ws = false,
+						fw_epsilon = fw_epsilon,
+						dual_gap = dual_gap,
+						rel_dual_gap = rel_dual_gap,
+						min_node_fw_epsilon = min_node_fw_epsilon,
+					)
+			end
+		else
+			x, _, result = Boscia.solve(
+				f,
+				grad!,
+				lmo,
+				verbose = verbose,
+				time_limit = time_limit,
+				lazy = true,
+				variant = variant,
+				use_DICG_warm_start = false,
+				fw_epsilon = fw_epsilon,
+				dual_gap = dual_gap,
+				rel_dual_gap = rel_dual_gap,
+				min_node_fw_epsilon = min_node_fw_epsilon,
+			)
+		end
+	else
+		if contains(mode, "ws")
+			if contains(mode, "full")
+				x, _, result =
+					Boscia.solve(
+						f,
+						grad!,
+						lmo,
+						verbose = verbose,
+						time_limit = time_limit,
+						lazy = false,
+						variant = variant,
+						use_DICG_warm_start = true,
+						use_full_ws = true,
+						fw_epsilon = fw_epsilon,
+						dual_gap = dual_gap,
+						rel_dual_gap = rel_dual_gap,
+						min_node_fw_epsilon = min_node_fw_epsilon,
+					)
+			else
+				x, _, result = Boscia.solve(
+					f,
+					grad!,
+					lmo,
+					verbose = verbose,
+					time_limit = time_limit,
+					lazy = false,
+					variant = variant,
+					use_DICG_warm_start = true,
+					fw_epsilon = fw_epsilon,
+					dual_gap = dual_gap,
+					rel_dual_gap = rel_dual_gap,
+					min_node_fw_epsilon = min_node_fw_epsilon,
+				)
+			end
+		else
+			x, _, result = Boscia.solve(
+				f,
+				grad!,
+				lmo,
+				verbose = verbose,
+				time_limit = time_limit,
+				lazy = false,
+				variant = variant,
+				use_DICG_warm_start = false,
+				fw_epsilon = fw_epsilon,
+				dual_gap = dual_gap,
+				rel_dual_gap = rel_dual_gap,
+				min_node_fw_epsilon = min_node_fw_epsilon,
+			)
+		end
+	end
+
+	return x, result
+end
+
+function cl_boscia(
+	seed,
+	boxes_per_img,
+	n_imgs;
+	mode = "custom",
+	int_vars = nothing,
+	lower_bounds = nothing,
+	upper_bounds = nothing,
+	dual_gap = 1e-6,
+	rel_dual_gap = 1e-2,
+	verbose = true,
+	time_limit = 3600,
+	write = true,
+	fw_epsilon = 1e-2,
+	min_node_fw_epsilon = 1e-6,
+)
+	@show seed
+	Random.seed!(seed)
+	dim = 630
+
+	f, grad! = build_objective()
+
+	sblmo = FlowPolytopeBLMO(10, 63)
+	lb = zeros(dim)
+	ub = ones(dim)
+	int_vars = collect(1:dim)
+	lmo = Boscia.ManagedBoundedLMO(sblmo, lb, ub, int_vars, dim)
+
+	custom_modes = ["custom", "custom_lazy", "custom_dicg", "custom_dicg_lazy", "custom_dicg_ws", "custom_dicg_lazy_ws", "custom_dicg_lazy_full_ws", "custom_dicg_full_ws"]
+	mip_modes = ["mip", "mip_lazy", "mip_dicg", "mip_dicg_lazy", "mip_dicg_ws", "mip_dicg_lazy_ws", "mip_dicg_lazy_full_ws", "mip_dicg_full_ws"]
+
+	if mode in custom_modes
+	elseif mode in mip_modes
+
+	else
+		error("Mode not known")
+	end
+
+	if contains(mode, "dicg")
+		x, result = run_args_boscia(mode, f, grad!, lmo, verbose, time_limit; variant = Boscia.DICG(), fw_epsilon = fw_epsilon, dual_gap = dual_gap, rel_dual_gap = rel_dual_gap, min_node_fw_epsilon = min_node_fw_epsilon)
+	else
+		x, result = run_args_boscia(mode, f, grad!, lmo, verbose, time_limit; variant = Boscia.BPCG(), fw_epsilon = fw_epsilon, dual_gap = dual_gap, rel_dual_gap = rel_dual_gap, min_node_fw_epsilon = min_node_fw_epsilon)
+	end
+
+	total_time_in_sec = result[:total_time_in_sec]
+	status = result[:status]
+	if occursin("Optimal", result[:status])
+		status = "OPTIMAL"
+	elseif occursin("Time", result[:status])
+		status = "TIME_LIMIT"
+	end
+
+
+	if write
+		lb_list = result[:list_lb]
+		ub_list = result[:list_ub]
+		time_list = result[:list_time]
+		list_lmo_calls = result[:list_lmo_calls_acc]
+		list_active_set_size_cb = result[:list_active_set_size]
+		list_discarded_set_size_cb = result[:list_discarded_set_size]
+		list_local_tightening = result[:local_tightenings]
+		list_global_tightening = result[:global_tightenings]
+
+		dir_path = joinpath(@__DIR__, "./csv")
+
+		mkpath(dir_path)
+
+		df_full = DataFrame(
+			time = time_list / 1000,
+			lowerBound = lb_list,
+			upperBound = ub_list,
+			termination = status,
+			LMOcalls = list_lmo_calls,
+			localTighteings = list_local_tightening,
+			globalTightenings = list_global_tightening,
+			list_active_set_size_cb = list_active_set_size_cb,
+			list_discarded_set_size_cb = list_discarded_set_size_cb,
+		)
+		file_name_full = joinpath(
+			dir_path,
+			"full_run_boscia_" * mode * "_.csv",
+		)
+		CSV.write(file_name_full, df_full, append = false)
+
+
+
+		@show result[:primal_objective]
+		df = DataFrame(
+			time = total_time_in_sec,
+			solution = result[:primal_objective],
+			dual_gap = result[:dual_gap],
+			rel_dual_gap = result[:rel_dual_gap],
+			termination = status,
+			ncalls = result[:lmo_calls],
+		)
+
+
+
+		file_name = joinpath(
+			dir_path,
+			"boscia_" * mode * "_" * ".csv",
+		)
+		CSV.write(file_name, df, append = false, writeheader = true)
+	end
+end
+
+cl_boscia(4, 20, 33; mode = "custom_lazy", time_limit = 3600, dual_gap = 1e-6, rel_dual_gap = 1e-1, fw_epsilon = 1e-3, min_node_fw_epsilon = 1e-6)
diff --git a/examples/co_localization_BLMO.jl b/examples/co_localization_BLMO.jl
new file mode 100644
index 000000000..9d34375ba
--- /dev/null
+++ b/examples/co_localization_BLMO.jl
@@ -0,0 +1,209 @@
+## FlowPolytope BLMO
+using Boscia
+using FrankWolfe
+using LinearAlgebra
+using SparseArrays
+using Random
+
+Random.seed!(4)
+
+struct FlowPolytopeBLMO <: Boscia.SimpleBoundableLMO
+	m::Int
+	n::Int
+	atol::Float64
+	rtol::Float64
+end
+
+FlowPolytopeBLMO(m, n) =
+	FlowPolytopeBLMO(m, n, 1e-6, 1e-3)
+
+"""
+Computes the extreme point given an direction d, the current lower and upper bounds on the integer variables, and the set of integer variables.
+"""
+function Boscia.bounded_compute_extreme_point(sblmo::FlowPolytopeBLMO, d, lb, ub, int_vars; kwargs...)
+	m = sblmo.m
+	n = sblmo.n
+	d = reshape(d, m, n)
+	v = spzeros(m * n)
+
+	for j in collect(1:n)
+		min = Inf
+		min_i = -1
+		for i in collect(1:m)
+			idx = Int((j - 1) * m + i)
+			int_vars_idx = findfirst(x -> x == idx, int_vars)
+			if int_vars_idx !== nothing
+				if lb[int_vars_idx] == 1.0
+					min_i = i
+					break
+				end
+				if ub[int_vars_idx] == 0.0
+					continue
+				end
+			end
+
+			if d[i, j] < min
+				min = d[i, j]
+				min_i = i
+			end
+		end
+
+		v[Int((j - 1) * m + min_i)] = 1.0
+
+	end
+	v = Vector(v)
+	return v
+end
+
+# sblmo = FlowPolytopeBLMO(3, 4)
+# d = rand(3, 4)
+# lb = zeros(12)
+# ub = ones(12)
+# lb[5] = 1.0
+# int_vars = collect(1:12)
+# v = Boscia.bounded_compute_extreme_point(sblmo, d, lb, ub, int_vars)
+# @show d
+# @show v
+
+
+
+"""
+Computes the inface extreme point given an direction d, x, the current lower and upper bounds on the integer variables, and the set of integer variables.
+"""
+function Boscia.bounded_compute_inface_extreme_point(
+	sblmo::FlowPolytopeBLMO,
+	d,
+	x,
+	lb,
+	ub,
+	int_vars;
+	kwargs...,
+)
+	m = sblmo.m
+	n = sblmo.n
+	d = reshape(d, m, n)
+	x = reshape(x, m, n)
+	v = spzeros(m * n)
+
+	for j in collect(1:n)
+		min = Inf
+		min_i = -1
+		for i in collect(1:m)
+			if x[i, j] ≈ 0.0
+				continue
+			end
+
+			if x[i, j] ≈ 1.0
+				min_i = i
+				break
+			end
+			idx = Int((j - 1) * m + i)
+			int_vars_idx = findfirst(x -> x == idx, int_vars)
+			if int_vars_idx !== nothing
+				if lb[int_vars_idx] == 1.0
+					min_i = i
+					break
+				end
+				if ub[int_vars_idx] == 0.0
+					continue
+				end
+			end
+
+			if d[i, j] < min
+				min = d[i, j]
+				min_i = i
+			end
+		end
+		v[Int((j - 1) * m + min_i)] = 1.0
+
+	end
+	v = Vector(v)
+	return v
+
+
+end
+
+"""
+LMO-like operation which computes a vertex minimizing in `direction` on the face defined by the current fixings.
+Fixings are maintained by the oracle (or deduced from `x` itself).
+"""
+function Boscia.bounded_dicg_maximum_step(
+	sblmo::FlowPolytopeBLMO,
+	direction,
+	x,
+	lb,
+	ub,
+	int_vars;
+	kwargs...,
+)
+	T = promote_type(eltype(x), eltype(direction))
+	gamma_max = one(T)
+	for idx in eachindex(x)
+		if direction[idx] != 0.0
+			# iterate already on the boundary
+			if (direction[idx] < 0 && x[idx] ≈ 1) || (direction[idx] > 0 && x[idx] ≈ 0)
+				return zero(gamma_max)
+			end
+			# clipping with the zero boundary
+			if direction[idx] > 0
+				gamma_max = min(gamma_max, x[idx] / direction[idx])
+			else
+				@assert direction[idx] < 0
+				gamma_max = min(gamma_max, -(1 - x[idx]) / direction[idx])
+			end
+		end
+	end
+	return gamma_max
+
+end
+
+"""
+"""
+function Boscia.is_decomposition_invariant_oracle_simple(sblmo::FlowPolytopeBLMO)
+	return true
+end
+
+"""
+"""
+function Boscia.dicg_split_vertices_set_simple(sblmo::FlowPolytopeBLMO, x, vidx)
+	x0_left = copy(x)
+	x0_right = copy(x)
+	return x0_left, x0_right
+end
+
+"""
+The sum of each row and column has to be equal to 1.
+"""
+function Boscia.is_simple_linear_feasible(sblmo::FlowPolytopeBLMO, v::AbstractVector)
+	m = sblmo.m
+	n = sblmo.n
+	v = reshape(v, m, n)
+	for j in 1:n
+
+
+		if !isapprox(sum(v[:, j]), 1.0)
+			@debug "Column sum not 1: $(sum(v[((j-1)*m+1):(j*m)]))"
+			return false
+		end
+
+		# # append by column ? row sum : column sum
+		# if !isapprox(sum(v[j:n:n^2]), 1.0, atol = 1e-6, rtol = 1e-3)
+		# 	@debug "Row sum not 1: $(sum(v[i:n:n^2]))"
+		# 	return false
+		# end
+	end
+	return true
+end
+
+
+function Boscia.build_dicg_start_point(lmo)
+	# We pick a random point.
+	n = lmo.blmo.n
+	d = ones(n)
+	x0 = FrankWolfe.compute_extreme_point(lmo, d)
+	return x0
+end
+
+function domain_orcal(x)
+	return true
+end
diff --git a/examples/experiments/examples b/examples/experiments/examples
new file mode 100644
index 000000000..30736e52f
--- /dev/null
+++ b/examples/experiments/examples
@@ -0,0 +1 @@
+read
diff --git a/examples/experiments/src_birkhoff.7z b/examples/experiments/src_birkhoff.7z
new file mode 100644
index 000000000..dcb36847f
Binary files /dev/null and b/examples/experiments/src_birkhoff.7z differ
diff --git a/examples/experiments/src_birkhoff.zip b/examples/experiments/src_birkhoff.zip
new file mode 100644
index 000000000..e95381fcc
Binary files /dev/null and b/examples/experiments/src_birkhoff.zip differ
diff --git a/src/Boscia.jl b/src/Boscia.jl
index 6bb6a11c5..511ebacce 100644
--- a/src/Boscia.jl
+++ b/src/Boscia.jl
@@ -3,6 +3,14 @@ module Boscia
 using FrankWolfe
 import FrankWolfe: compute_extreme_point
 export compute_extreme_point
+import FrankWolfe: is_decomposition_invariant_oracle
+export is_decomposition_invariant_oracle
+
+import FrankWolfe: compute_inface_extreme_point
+export compute_inface_extreme_point
+
+import FrankWolfe: dicg_maximum_step
+export dicg_maximum_step
 using Random
 using LinearAlgebra
 import Bonobo
diff --git a/src/MOI_bounded_oracle.jl b/src/MOI_bounded_oracle.jl
index 89b3f20d2..153d502f7 100644
--- a/src/MOI_bounded_oracle.jl
+++ b/src/MOI_bounded_oracle.jl
@@ -1,12 +1,25 @@
 """
     BoundedLinearMinimizationOracle for solvers supporting MathOptInterface.
 """
+
+# Store extra information of solving inface extrem points.
+# The keys of MOI_attribute should correspond to specific MOI_attribute names.
+mutable struct Inface_point_solve_data
+    MOI_attribute::Dict
+    function Inface_point_solve_data()
+        MOI_attribute = Dict()
+        return new(MOI_attribute)
+    end
+end
+
 struct MathOptBLMO{OT<:MOI.AbstractOptimizer} <: BoundedLinearMinimizationOracle
     o::OT
     use_modify::Bool
+    inface_point_solve_data::Inface_point_solve_data
     function MathOptBLMO(o, use_modify=true)
         MOI.set(o, MOI.ObjectiveSense(), MOI.MIN_SENSE)
-        return new{typeof(o)}(o, use_modify)
+        inface_point_solve_data = Inface_point_solve_data()
+        return new{typeof(o)}(o, use_modify, inface_point_solve_data)
     end
 end
 
@@ -342,9 +355,16 @@ end
 Get solve time, number of nodes and number of simplex iterations.
 """
 function get_BLMO_solve_data(blmo::MathOptBLMO)
-    opt_times = MOI.get(blmo.o, MOI.SolveTimeSec())
-    numberofnodes = MOI.get(blmo.o, MOI.NodeCount())
-    simplex_iterations = MOI.get(blmo.o, MOI.SimplexIterations())
+    if !isempty(blmo.inface_point_solve_data.MOI_attribute)
+        opt_times = blmo.inface_point_solve_data.MOI_attribute[MOI.SolveTimeSec()]
+        numberofnodes = blmo.inface_point_solve_data.MOI_attribute[MOI.NodeCount()]
+        simplex_iterations = blmo.inface_point_solve_data.MOI_attribute[MOI.SimplexIterations()]
+        empty!(blmo.inface_point_solve_data.MOI_attribute)
+    else
+        opt_times = MOI.get(blmo.o, MOI.SolveTimeSec())
+        numberofnodes = MOI.get(blmo.o, MOI.NodeCount())
+        simplex_iterations = MOI.get(blmo.o, MOI.SimplexIterations())
+    end
     return opt_times, numberofnodes, simplex_iterations
 end
 
@@ -599,6 +619,50 @@ function Bonobo.get_branching_variable(
     return max_idx
 end
 
+
+function is_decomposition_invariant_oracle(blmo::MathOptBLMO)
+    true
+end
+
+function compute_inface_extreme_point(blmo::MathOptBLMO, direction, x; kwargs...)
+    MOI_attribute = Dict()
+    MOI_attribute[MOI.SolveTimeSec()] = 0.0
+    MOI_attribute[MOI.NodeCount()] = 0.0
+    MOI_attribute[MOI.SimplexIterations()] = 0.0
+    blmo.inface_point_solve_data.MOI_attribute = MOI_attribute
+    lmo = convert(FrankWolfe.MathOptLMO, blmo)
+    a = FrankWolfe.compute_inface_extreme_point(lmo,
+        direction,
+        x;
+        solve_data=blmo.inface_point_solve_data.MOI_attribute,
+        kwargs,
+    )
+    @assert blmo isa MathOptBLMO
+    return a
+end
+
+function dicg_maximum_step(blmo::MathOptBLMO, direction, x; kwargs...)
+    lmo = convert(FrankWolfe.MathOptLMO, blmo)
+    return FrankWolfe.dicg_maximum_step(
+        lmo,
+        direction,
+        x;
+        kwargs...,
+    )
+end
+
+# Provide specific active_set split method for MathOptBLMO in DICG.
+function dicg_split_vertices_set!(blmo::MathOptBLMO, active_set::FrankWolfe.ActiveSet{T,R}, tree, vidx::Int;kwargs...)where {T,R}
+    x = FrankWolfe.get_active_set_iterate(active_set)
+    x0_left = copy(x)
+    x0_right = copy(x)
+    x0_left[vidx] = floor(x[vidx])
+    x0_right[vidx] = ceil(x[vidx])
+    as_left = FrankWolfe.ActiveSet([(1.0, x0_left)])
+    as_right = FrankWolfe.ActiveSet([(1.0, x0_right)])
+    return as_left, as_right
+end
+
 """
 Solve function in case of MathOptLMO. 
 Converts the lmo into a MathOptBLMO and calls the solve function below.
diff --git a/src/blmo_interface.jl b/src/blmo_interface.jl
index be6c4a6cd..e50d174bf 100644
--- a/src/blmo_interface.jl
+++ b/src/blmo_interface.jl
@@ -185,6 +185,13 @@ function get_variables_pointers(blmo::BoundedLinearMinimizationOracle, tree)
     return collect(1:N)
 end
 
+"""
+Given split variable index, return corresponding active sets for DICG as a warm start.
+"""
+function dicg_split_vertices_set!(blmo::BoundedLinearMinimizationOracle, active_set::FrankWolfe.ActiveSet{T,R}, tree, vidx::Int;kwargs...)where {T,R}
+    error("Not implemented yet")
+end
+
 
 ## Logs
 """
diff --git a/src/callbacks.jl b/src/callbacks.jl
index a7aee6d32..e898df074 100644
--- a/src/callbacks.jl
+++ b/src/callbacks.jl
@@ -7,89 +7,137 @@ Time limit is checked.
 If the vertex is providing a better incumbent, it is added as solution.
 """
 function build_FW_callback(
-    tree,
-    min_number_lower,
-    check_rounding_value::Bool,
-    fw_iterations,
-    min_fw_iterations,
-    time_ref,
-    time_limit,
+	tree,
+	min_number_lower,
+	check_rounding_value::Bool,
+	fw_iterations,
+	min_fw_iterations,
+	time_ref,
+	time_limit;
+	use_DICG = false,
 )
-    vars = get_variables_pointers(tree.root.problem.tlmo.blmo, tree)
-    # variable to only fetch heuristics when the counter increases
-    ncalls = -1
-    return function fw_callback(state, active_set, kwargs...)
-        @assert isapprox(sum(active_set.weights), 1.0)
-        @assert sum(active_set.weights .< 0) == 0
-        # TODO deal with vertices becoming infeasible with conflicts
-        @debug begin
-            if !is_linear_feasible(tree.root.problem.tlmo, state.v)
-                @info "$(state.v)"
-                check_infeasible_vertex(tree.root.problem.tlmo.blmo, tree)
-                @assert is_linear_feasible(tree.root.problem.tlmo, state.v)
-            end  
-            if state.step_type != FrankWolfe.ST_SIMPLEXDESCENT && !is_integer_feasible(tree, state.v)
-                @info "Vertex not integer feasible! Here are the integer variables: $(state.v[tree.root.problem.integer_variables])"
-                @assert is_integer_feasible(tree, state.v)
-            end    
-        end
-        push!(fw_iterations, state.t)
+	vars = get_variables_pointers(tree.root.problem.tlmo.blmo, tree)
+	# variable to only fetch heuristics when the counter increases
+	ncalls = -1
+	if !use_DICG
+		return function (state, active_set, kwargs...)
+			return process_FW_callback_logic(
+				tree, state, vars, fw_iterations, ncalls, min_fw_iterations,
+				min_number_lower, time_ref, time_limit, use_DICG; active_set = active_set, kwargs,
+			)
+		end
+	else
+		return function (state, pre_computed_set, kwargs...)
+			return process_FW_callback_logic(
+				tree, state, vars, fw_iterations, ncalls, min_fw_iterations,
+				min_number_lower, time_ref, time_limit, use_DICG; pre_computed_set = pre_computed_set, kwargs,
+			)
+		end
+	end
+end
 
-        if state.lmo !== nothing  # can happen with using Blended Conditional Gradient
-            if ncalls != state.lmo.ncalls
-                ncalls = state.lmo.ncalls
-                (best_v, best_val) = find_best_solution(
-                    tree.root.problem.f,
-                    tree.root.problem.tlmo.blmo,
-                    vars,
-                    tree.root.options[:domain_oracle],
-                )
-                if best_val < tree.incumbent
-                    node = tree.nodes[tree.root.current_node_id[]]
-                    add_new_solution!(tree, node, best_val, best_v, :Solver)
-                    Bonobo.bound!(tree, node.id)
-                end
-            end
-        end
+function process_FW_callback_logic(
+	tree,
+	state,
+	vars,
+	fw_iterations,
+	ncalls,
+	min_fw_iterations,
+	min_number_lower,
+	time_ref,
+	time_limit,
+	use_DICG;
+	active_set = nothing,
+	pre_computed_set = nothing,
+	kwargs...,
+)
 
-        if (state.primal - state.dual_gap > tree.incumbent + 1e-2) &&
-           tree.num_nodes != 1 &&
-           state.t > min_fw_iterations
-            return false
-        end
+	if !use_DICG
+		@assert isapprox(sum(active_set.weights), 1.0)
+		@assert sum(active_set.weights .< 0) == 0
+	end
 
-        if tree.root.options[:domain_oracle](state.v) && state.step_type != FrankWolfe.ST_SIMPLEXDESCENT
-            val = tree.root.problem.f(state.v)
-            if val < tree.incumbent
-                #TODO: update solution without adding node
-                node = tree.nodes[tree.root.current_node_id[]]
-                add_new_solution!(tree, node, val, copy(state.v), :vertex)
-                Bonobo.bound!(tree, node.id)
-            end
-        end
+	# TODO deal with vertices becoming infeasible with conflicts
+	@debug begin
+		if !is_linear_feasible(tree.root.problem.tlmo, state.v)
+			@info "$(state.v)"
+			check_infeasible_vertex(tree.root.problem.tlmo.blmo, tree)
+			@assert is_linear_feasible(tree.root.problem.tlmo, state.v)
+		end
+		if state.step_type != FrankWolfe.ST_SIMPLEXDESCENT && !is_integer_feasible(tree, state.v)
+			@info "Vertex not integer feasible! Here are the integer variables: $(state.v[tree.root.problem.integer_variables])"
+			@assert is_integer_feasible(tree, state.v)
+		end
+	end
+	push!(fw_iterations, state.t)
 
-        node = tree.nodes[tree.root.current_node_id[]]
-        if length(node.active_set) > 1 &&
-           !isempty(tree.nodes) &&
-           min_number_lower <= length(values(tree.nodes))
-            counter = 0
-            for n in values(tree.nodes)
-                if n.lb < val
-                    counter += 1
-                end
-                if counter > min_number_lower
-                    return false
-                end
-            end
-        end
+	if state.lmo !== nothing  # can happen with using Blended Conditional Gradient
+		if ncalls != state.lmo.ncalls
+			ncalls = state.lmo.ncalls
+			(best_v, best_val) = find_best_solution(
+				tree.root.problem.f,
+				tree.root.problem.tlmo.blmo,
+				vars,
+				tree.root.options[:domain_oracle],
+			)
+			if best_val < tree.incumbent
+				node = tree.nodes[tree.root.current_node_id[]]
+				add_new_solution!(tree, node, best_val, best_v, :Solver)
+				Bonobo.bound!(tree, node.id)
+			end
+		end
+	end
+
+	if (state.primal - state.dual_gap > tree.incumbent + 1e-2) &&
+	   tree.num_nodes != 1 &&
+	   state.t > min_fw_iterations
+		return false
+	end
 
-        # check for time limit
-        if isfinite(time_limit) && Dates.now() >= time_ref + Dates.Second(time_limit)
-            return false
+	if tree.root.options[:domain_oracle](state.v) && state.step_type != FrankWolfe.ST_SIMPLEXDESCENT
+		val = tree.root.problem.f(state.v)
+		if val < tree.incumbent
+			#TODO: update solution without adding node
+			node = tree.nodes[tree.root.current_node_id[]]
+			add_new_solution!(tree, node, val, copy(state.v), :vertex)
+			Bonobo.bound!(tree, node.id)
+		end
+	end
+
+	node = tree.nodes[tree.root.current_node_id[]]
+	if length(node.active_set) > 1 &&
+	   !isempty(tree.nodes) &&
+	   min_number_lower <= length(values(tree.nodes))
+		counter = 0
+		for n in values(tree.nodes)
+			if n.lb < val
+				counter += 1
+			end
+			if counter > min_number_lower
+				return false
+			end
+		end
+	end
+
+	# check for time limit
+	if isfinite(time_limit) && Dates.now() >= time_ref + Dates.Second(time_limit)
+		return false
+	end
+
+        if pre_computed_set !== nothing
+                if state.step_type !== FrankWolfe.last || state.step_type !== FrankWolfe.pp
+                    idx = findfirst(x -> x == state.v, pre_computed_set)
+                    if idx == nothing
+                        if length(pre_computed_set) > (length(state.v) + 1)
+                            deleteat!(pre_computed_set, 1)
+                        end
+                        push!(pre_computed_set, state.v)
+                    end
+                end
         end
+        
+	return true
 
-        return true
-    end
 end
 
 """
diff --git a/src/custom_bonobo.jl b/src/custom_bonobo.jl
index b1966ec54..abadd136c 100644
--- a/src/custom_bonobo.jl
+++ b/src/custom_bonobo.jl
@@ -90,6 +90,7 @@ function Bonobo.optimize!(
             0,
             0,
             0.0,
+            [],
         )
         callback(tree, dummy_node, node_infeasible=true)
     end
diff --git a/src/frank_wolfe_variants.jl b/src/frank_wolfe_variants.jl
index cd3e3f21c..16fc502fa 100644
--- a/src/frank_wolfe_variants.jl
+++ b/src/frank_wolfe_variants.jl
@@ -4,7 +4,7 @@ Frank-Wolfe variant used to compute the problems at node level.
 A `FrankWolfeVariant` must implement
 ```
 solve_frank_wolfe(fw::FrankWolfeVariant, f, grad!, lmo, active_set, line_search, epsilon, max_iteration,
-    added_dropped_vertices, use_extra_vertex_storage, callback, lazy, timeout, verbose, workspace))
+	added_dropped_vertices, use_extra_vertex_storage, callback, lazy, timeout, verbose, workspace))
 ```
 It may also implement `build_frank_wolfe_workspace(x)` which creates a
 workspace structure that is passed as last argument to `solve_frank_wolfe`.
@@ -15,8 +15,8 @@ abstract type FrankWolfeVariant end
 Base.print(io::IO, fw::FrankWolfeVariant) = print(io, split(string(typeof(fw)), ".")[end])
 
 """
-    solve_frank_wolfe(fw::FrankWolfeVariant, f, grad!, lmo, active_set, line_search, epsilon, max_iteration,
-    added_dropped_vertices, use_extra_vertex_storage, callback, lazy, timeout, verbose, workspace)
+	solve_frank_wolfe(fw::FrankWolfeVariant, f, grad!, lmo, active_set, line_search, epsilon, max_iteration,
+	added_dropped_vertices, use_extra_vertex_storage, callback, lazy, timeout, verbose, workspace)
 
 Returns the optimal solution x to the node problem, its primal and dual gap and the active set. 
 """
@@ -26,7 +26,7 @@ build_frank_wolfe_workspace(::FrankWolfeVariant, x) = nothing
 
 
 """
-    Away-Frank-Wolfe
+	Away-Frank-Wolfe
 
 In every iteration, it computes the worst performing vertex, called away vertex, in the active set with regard to the gradient.
 If enough local progress can be made, weight is shifted from the away vertex to all other vertices. 
@@ -53,6 +53,7 @@ function solve_frank_wolfe(
     timeout=Inf,
     verbose=false,
     workspace=nothing,
+    pre_computed_set=nothing,
 )
     x, _, primal, dual_gap, _, active_set = FrankWolfe.away_frank_wolfe(
         f,
@@ -79,7 +80,7 @@ Base.print(io::IO, ::AwayFrankWolfe) = print(io, "Away-Frank-Wolfe")
 
 
 """
-    Blended Conditional Gradient
+	Blended Conditional Gradient
 """
 struct Blended <: FrankWolfeVariant end
 
@@ -101,6 +102,7 @@ function solve_frank_wolfe(
     timeout=Inf,
     verbose=false,
     workspace=nothing,
+    pre_computed_set=nothing,
 )
     x, _, primal, dual_gap, _, active_set = blended_conditional_gradient(
         f,
@@ -125,7 +127,7 @@ end
 Base.print(io::IO, ::Blended) = print(io, "Blended Conditional Gradient")
 
 """
-    Blended Pairwise Conditional Gradient
+	Blended Pairwise Conditional Gradient
 """
 struct BPCG <: FrankWolfeVariant end
 
@@ -147,6 +149,7 @@ function solve_frank_wolfe(
     timeout=Inf,
     verbose=false,
     workspace=nothing,
+    pre_computed_set=nothing,
 )
     x, _, primal, dual_gap, _, active_set = FrankWolfe.blended_pairwise_conditional_gradient(
         f,
@@ -165,18 +168,84 @@ function solve_frank_wolfe(
         timeout=timeout,
         verbose=verbose,
     )
-
     return x, primal, dual_gap, active_set
 end
 
 Base.print(io::IO, ::BPCG) = print(io, "Blended Pairwise Conditional Gradient")
 
 """
-    Vanilla-Frank-Wolfe
+   DICG-Frank-Wolfe
+
+The Decomposition-invariant Frank-Wolfe. 
+
+"""
+struct DICG <: FrankWolfeVariant end
+
+function solve_frank_wolfe(
+    frank_wolfe_variant::DICG,
+    f,
+    grad!,
+    lmo,
+    active_set;
+    line_search::FrankWolfe.LineSearchMethod=FrankWolfe.Adaptive(),
+    epsilon=1e-7,
+    max_iteration=10000,
+    add_dropped_vertices=false,
+    use_extra_vertex_storage=false,
+    extra_vertex_storage=nothing,
+    callback=nothing,
+    lazy=false,
+    lazy_tolerance=2.0,
+    timeout=Inf,
+    verbose=false,
+    workspace=nothing,
+    pre_computed_set=nothing,
+)
+    # We keep track of computed extreme points by creating logging callback.
+    function make_callback(pre_computed_set)
+	return function DICG_callback(state)
+		if !callback(state, pre_computed_set)
+			return false
+		end
+		return true
+	end
+    end
+    
+    x0 = dicg_start_point_initialize(lmo, active_set, pre_computed_set)
+	
+    if x0 == nothing
+	return NaN, Inf, Inf, pre_computed_set
+    else
+	@assert is_linear_feasible(lmo, x0)
+    end
 
-The standard variant of Frank-Wolfe. In each iteration, the vertex v minimizing ∇f * (x-v) is computed. 
+    DICG_callback = make_callback(pre_computed_set)
+
+    x, _, primal, dual_gap, _ = FrankWolfe.decomposition_invariant_conditional_gradient(
+        f,
+        grad!,
+        lmo,
+        x0;
+        line_search=line_search,
+        epsilon=epsilon,
+        max_iteration=max_iteration,
+        verbose=verbose,
+        timeout=timeout,
+        lazy=lazy,
+        linesearch_workspace=workspace,
+        lazy_tolerance=lazy_tolerance,
+        callback=DICG_callback,
+        extra_vertex_storage=pre_computed_set,
+    )
+
+    return x, primal, dual_gap, pre_computed_set
+end
+
+Base.print(io::IO, ::DICG) = print(io, "Decompostion-Invariant-Frank-Wolfe")
+
+"""
+	Vanilla-Frank-Wolfe
 
-Lazification cannot be used in this setting.
 """
 struct VanillaFrankWolfe <: FrankWolfeVariant end
 
@@ -198,27 +267,28 @@ function solve_frank_wolfe(
     timeout=Inf,
     verbose=false,
     workspace=nothing,
+    pre_computed_set=nothing,
 )
-    # If the flag away_steps is set to false, away_frank_wolfe performs Vanilla.
     # Observe that the lazy flag is only observed if away_steps is set to true, so it can neglected. 
     x, _, primal, dual_gap, _, active_set = FrankWolfe.away_frank_wolfe(
         f,
         grad!,
         lmo,
         active_set,
-        away_steps=false,
         epsilon=epsilon,
         max_iteration=max_iteration,
         line_search=line_search,
         callback=callback,
+        lazy=lazy,
+        lazy_tolerance=lazy_tolerance,
         timeout=timeout,
         add_dropped_vertices=add_dropped_vertices,
         use_extra_vertex_storage=use_extra_vertex_storage,
         extra_vertex_storage=extra_vertex_storage,
         verbose=verbose,
+        away_steps=false,
     )
-
-    return x, primal, dual_gap, active_set
+	return x, primal, dual_gap, active_set
 end
 
 Base.print(io::IO, ::VanillaFrankWolfe) = print(io, "Vanilla-Frank-Wolfe")
diff --git a/src/interface.jl b/src/interface.jl
index 92e161ea9..0ca9a041f 100644
--- a/src/interface.jl
+++ b/src/interface.jl
@@ -109,8 +109,9 @@ function solve(
     use_shadow_set=true,
     custom_heuristics=[Heuristic()],
     rounding_prob=1.0,
-    clean_solutions=false, 
+    clean_solutions=false,
     max_clean_iter=10,
+    use_DICG_warm_start=false,
     kwargs...,
 )
     if verbose
@@ -129,6 +130,8 @@ function solve(
         println("\t Additional kwargs: ", join(keys(kwargs), ","))
     end
 
+
+    println("Boscia-version")
     n, _ = get_list_of_variables(blmo)
 
     integer_variables = Vector{Int}()
@@ -174,6 +177,8 @@ function solve(
     end
     vertex_storage = FrankWolfe.DeletedVertexStorage(typeof(v)[], 1)
 
+    pre_computed_set = use_DICG_warm_start ? [v] : nothing
+
     m = SimpleOptimizationProblem(f, grad!, n, integer_variables, time_lmo, global_bounds)
     nodeEx = FrankWolfeNode(
         NodeInfo(1, f(v), f(v)),
@@ -187,6 +192,7 @@ function solve(
         0,
         0,
         0.0,
+        [v],
     )
 
     # Create standard heuristics
@@ -230,7 +236,7 @@ function solve(
                 :usePostsolve => use_postsolve,
                 :variant => variant,
                 :use_shadow_set => use_shadow_set,
-                :heuristics => heuristics, 
+                :heuristics => heuristics,
                 :heu_ncalls => 0,
                 :max_clean_iter => max_clean_iter,
                 :clean_solutions => clean_solutions,
@@ -254,6 +260,7 @@ function solve(
             local_tightenings=0,
             local_potential_tightenings=0,
             dual_gap=-Inf,
+            pre_computed_set=pre_computed_set,
         ),
     )
 
@@ -313,7 +320,16 @@ function solve(
         num_int,
     )
 
-    fw_callback = build_FW_callback(tree, min_number_lower, true, fw_iterations, min_fw_iterations, time_ref, tree.root.options[:time_limit])
+    fw_callback = build_FW_callback(
+        tree, 
+        min_number_lower, 
+        true, 
+        fw_iterations, 
+        min_fw_iterations, 
+        time_ref, 
+        tree.root.options[:time_limit], 
+        use_DICG=tree.root.options[:variant] == DICG(),
+    )
 
     tree.root.options[:callback] = fw_callback
     tree.root.current_node_id[] = Bonobo.get_next_node(tree, tree.options.traverse_strategy).id
diff --git a/src/managed_blmo.jl b/src/managed_blmo.jl
index 3a580c716..95439a556 100644
--- a/src/managed_blmo.jl
+++ b/src/managed_blmo.jl
@@ -69,6 +69,53 @@ function compute_extreme_point(blmo::ManagedBoundedLMO, d; kwargs...)
     return v
 end
 
+function is_decomposition_invariant_oracle(blmo::ManagedBoundedLMO)
+    return is_decomposition_invariant_oracle_simple(blmo.simple_lmo)
+end
+
+# Provide FrankWolfe.compute_inface_extreme_point
+function compute_inface_extreme_point(blmo::ManagedBoundedLMO, direction, x; kwargs...)
+    time_ref = Dates.now()
+    a = bounded_compute_inface_extreme_point(
+                blmo.simple_lmo,
+                direction,
+                x,
+                blmo.lower_bounds,
+                blmo.upper_bounds,
+                blmo.int_vars,
+                )
+    
+    blmo.solving_time = float(Dates.value(Dates.now() - time_ref))
+    return a
+end
+
+# Check if the given point a is on the minimal face of x
+function is_inface_feasible(blmo::ManagedBoundedLMO, a, x)
+	return is_simple_inface_feasible(blmo.simple_lmo, a, x, blmo.lower_bounds, blmo.upper_bounds, blmo.int_vars)
+end
+
+#Provide FrankWolfe.dicg_maximum_step
+function dicg_maximum_step(blmo::ManagedBoundedLMO, x, direction; kwargs...)
+    return bounded_dicg_maximum_step(
+                blmo.simple_lmo,
+                x,
+                direction,
+                blmo.lower_bounds, 
+                blmo.upper_bounds, 
+                blmo.int_vars,
+                )
+end
+
+# Provide specific active_set split method for simple_lmo in DICG.
+function dicg_split_vertices_set!(blmo::ManagedBoundedLMO, active_set::FrankWolfe.ActiveSet{T,R}, tree, vidx::Int;kwargs...)where {T,R}
+    x = FrankWolfe.get_active_set_iterate(active_set)
+    x0_left, x0_right = dicg_split_vertices_set_simple(blmo.simple_lmo, x, vidx)
+    as_left = FrankWolfe.ActiveSet([(1.0, x0_left)])
+    as_right = FrankWolfe.ActiveSet([(1.0, x0_right)])
+    return as_left, as_right
+end
+
+
 # Read global bounds from the problem.
 function build_global_bounds(blmo::ManagedBoundedLMO, integer_variables)
     global_bounds = IntegerBounds()
diff --git a/src/node.jl b/src/node.jl
index 1de2146fb..bf5a6f660 100644
--- a/src/node.jl
+++ b/src/node.jl
@@ -14,14 +14,14 @@ This needs to be added by every `AbstractNode` as `std::NodeInfo`
 
 This variant is more flexibel than Bonobo.BnBNodeInfo.
 """
-mutable struct NodeInfo{T<:Real} 
-    id :: Int
-    lb :: T
-    ub :: T 
+mutable struct NodeInfo{T<:Real}
+    id::Int
+    lb::T
+    ub::T
 end
 
-function Base.convert(::Type{NodeInfo{T}}, std::Bonobo.BnBNodeInfo) where T<:Real
-    return NodeInfo(std.id, T(std.lb), T(std.ub)) 
+function Base.convert(::Type{NodeInfo{T}}, std::Bonobo.BnBNodeInfo) where {T<:Real}
+    return NodeInfo(std.id, T(std.lb), T(std.ub))
 end
 
 """
@@ -40,6 +40,7 @@ A node in the branch-and-bound tree storing information for a Frank-Wolfe subpro
     All other integer bounds are stored in the root.
 'level' stores the level in the tree
 'fw_dual_gap_limit' set the tolerance for the dual gap in the FW algorithms
+'pre_computed_set' stores specifically the extreme points computed in DICG for warm-start.
 """
 mutable struct FrankWolfeNode{
     AT<:FrankWolfe.ActiveSet,
@@ -58,8 +59,8 @@ mutable struct FrankWolfeNode{
     local_tightenings::Int
     local_potential_tightenings::Int
     dual_gap::Float64
+    pre_computed_set::Any
 end
-
 """
 Create the information of the new branching nodes 
 based on their parent and the index of the branching variable
@@ -68,7 +69,6 @@ function Bonobo.get_branching_nodes_info(tree::Bonobo.BnBTree, node::FrankWolfeN
     if !is_valid_split(tree, vidx)
         error("Splitting on the same index as parent! Abort!")
     end
-
     # get iterate, primal and lower bound
     x = Bonobo.get_relaxed_values(tree, node)
     primal = tree.root.problem.f(x)
@@ -78,35 +78,59 @@ function Bonobo.get_branching_nodes_info(tree::Bonobo.BnBTree, node::FrankWolfeN
     # In case of strong convexity, check if a child can be pruned
     prune_left, prune_right = prune_children(tree, node, lower_bound_base, x, vidx)
 
-    # Split active set
-    active_set_left, active_set_right =
-        split_vertices_set!(node.active_set, tree, vidx, node.local_bounds)
+    #different ways to split active set
+    if tree.root.options[:variant] != DICG()
+
+        # Keep the same pre_computed_set
+        pre_computed_set_left, pre_computed_set_right =
+            node.pre_computed_set, node.pre_computed_set
+
+        # Split active set
+        active_set_left, active_set_right =
+            split_vertices_set!(node.active_set, tree, vidx, node.local_bounds)
+    else
+
+        if node.pre_computed_set !== nothing
+            # Split pre_computed_set
+            pre_computed_set_left, pre_computed_set_right =
+                split_pre_computed_set!(x, node.pre_computed_set, tree, vidx, node.local_bounds)
+        else
+            pre_computed_set_left, pre_computed_set_right = node.pre_computed_set, node.pre_computed_set
+        end
+        # Only support SBLMO polytope.
+        # User should implement specific dicg_split_vertices_set!() for different polytopes.
+        active_set_left, active_set_right =
+            dicg_split_vertices_set!(tree.root.problem.tlmo.blmo, node.active_set, tree, vidx)
+    end
+
     discarded_set_left, discarded_set_right =
         split_vertices_set!(node.discarded_vertices, tree, vidx, x, node.local_bounds)
 
-    # Sanity check
-    @assert isapprox(sum(active_set_left.weights), 1.0) "sum weights left: $(sum(active_set_left.weights))"
-    @assert sum(active_set_left.weights .< 0) == 0
-    for v in active_set_left.atoms
-        if !(v[vidx] <= floor(x[vidx]) + tree.options.atol)
-            error("active_set_left\n$(v)\n$vidx, $(x[vidx]), $(v[vidx])")
+    if tree.root.options[:variant] != DICG()
+        # Sanity check
+        @assert isapprox(sum(active_set_left.weights), 1.0) "sum weights left: $(sum(active_set_left.weights))"
+        @assert sum(active_set_left.weights .< 0) == 0
+        for v in active_set_left.atoms
+            if !(v[vidx] <= floor(x[vidx]) + tree.options.atol)
+                error("active_set_left\n$(v)\n$vidx, $(x[vidx]), $(v[vidx])")
+            end
         end
-    end
-    @assert isapprox(sum(active_set_right.weights), 1.0) "sum weights right: $(sum(active_set_right.weights))"
-    @assert sum(active_set_right.weights .< 0) == 0
-    for v in active_set_right.atoms
-        if !(v[vidx] >= ceil(x[vidx]) - tree.options.atol)
-            error("active_set_right\n$(v)\n$vidx, $(x[vidx]), $(v[vidx])")
+        @assert isapprox(sum(active_set_right.weights), 1.0) "sum weights right: $(sum(active_set_right.weights))"
+        @assert sum(active_set_right.weights .< 0) == 0
+        for v in active_set_right.atoms
+            if !(v[vidx] >= ceil(x[vidx]) - tree.options.atol)
+                error("active_set_right\n$(v)\n$vidx, $(x[vidx]), $(v[vidx])")
+            end
         end
-    end
-    for v in discarded_set_left.storage
-        if !(v[vidx] <= floor(x[vidx]) + tree.options.atol)
-            error("storage left\n$(v)\n$vidx, $(x[vidx]), $(v[vidx])")
+        for v in discarded_set_left.storage
+            if !(v[vidx] <= floor(x[vidx]) + tree.options.atol)
+                error("storage left\n$(v)\n$vidx, $(x[vidx]), $(v[vidx])")
+            end
         end
-    end
-    for v in discarded_set_right.storage
-        if !(v[vidx] >= ceil(x[vidx]) - tree.options.atol)
-            error("storage right\n$(v)\n$vidx, $(x[vidx]), $(v[vidx])")
+        for v in discarded_set_right.storage
+            if !(v[vidx] >= ceil(x[vidx]) - tree.options.atol)
+                error("storage right\n$(v)\n$vidx, $(x[vidx]), $(v[vidx])")
+            end
         end
     end
 
@@ -151,6 +175,7 @@ function Bonobo.get_branching_nodes_info(tree::Bonobo.BnBTree, node::FrankWolfeN
         local_tightenings=0,
         local_potential_tightenings=0,
         dual_gap=NaN,
+        pre_computed_set=pre_computed_set_left,
     )
     node_info_right = (
         active_set=active_set_right,
@@ -163,6 +188,7 @@ function Bonobo.get_branching_nodes_info(tree::Bonobo.BnBTree, node::FrankWolfeN
         local_tightenings=0,
         local_potential_tightenings=0,
         dual_gap=NaN,
+        pre_computed_set=pre_computed_set_right,
     )
 
     domain_right = !isempty(active_set_right)
@@ -229,20 +255,22 @@ function Bonobo.evaluate_node!(tree::Bonobo.BnBTree, node::FrankWolfeNode)
         error("Feasible region unbounded! Please check your constraints!")
         return NaN, NaN
     end
-
-    # Check feasibility of the iterate
-    active_set = node.active_set
-    x = FrankWolfe.compute_active_set_iterate!(node.active_set)
-    @assert is_linear_feasible(tree.root.problem.tlmo, x)
-    for (_, v) in node.active_set
-        @assert is_linear_feasible(tree.root.problem.tlmo, v)
+    
+    if tree.root.options[:variant] != DICG()
+        # Check feasibility of the iterate
+        active_set = node.active_set
+        x = FrankWolfe.compute_active_set_iterate!(node.active_set)
+        @assert is_linear_feasible(tree.root.problem.tlmo, x)
+        for (_, v) in node.active_set
+            @assert is_linear_feasible(tree.root.problem.tlmo, v)
+        end
     end
 
     # time tracking FW
     time_ref = Dates.now()
     domain_oracle = tree.root.options[:domain_oracle]
 
-    x, primal, dual_gap, active_set = solve_frank_wolfe(
+    x, primal, dual_gap, atoms_set = solve_frank_wolfe(
         tree.root.options[:variant],
         tree.root.problem.f,
         tree.root.problem.g,
@@ -258,14 +286,25 @@ function Bonobo.evaluate_node!(tree::Bonobo.BnBTree, node::FrankWolfeNode)
         extra_vertex_storage=node.discarded_vertices,
         callback=tree.root.options[:callback],
         verbose=tree.root.options[:fwVerbose],
+        pre_computed_set=node.pre_computed_set,
     )
 
+    if typeof(atoms_set).name.wrapper == FrankWolfe.ActiveSet
+        # update active set of the node
+        node.active_set = atoms_set
+    else
+	# update set of computed atoms and active set
+        if isa(x, Vector)
+            node.pre_computed_set = atoms_set
+            node.active_set = FrankWolfe.ActiveSet([(1.0, x)])
+        else
+            return NaN, NaN
+        end
+    end
+    
     node.fw_time = Dates.now() - time_ref
     node.dual_gap = dual_gap
 
-    # update active set of the node
-    node.active_set = active_set
-
     # tightening bounds at node level
     dual_tightening(tree, node, x, dual_gap)
 
diff --git a/src/polytope_blmos.jl b/src/polytope_blmos.jl
index 5eae736e4..b0212d0f4 100644
--- a/src/polytope_blmos.jl
+++ b/src/polytope_blmos.jl
@@ -27,7 +27,7 @@ end
 
 function is_simple_linear_feasible(sblmo::CubeSimpleBLMO, v)
     for i in setdiff(eachindex(v), sblmo.int_vars)
-        if !(sblmo.lower_bounds[i] ≤ v[i] + 1e-6 || !(v[i] - 1e-6 ≤ blmo.upper_bounds[i]))
+        if !(sblmo.lower_bounds[i] ≤ v[i] + 1e-6) || !(v[i] - 1e-6 ≤ sblmo.upper_bounds[i])
             @debug(
                 "Vertex entry: $(v[i]) Lower bound: $(blmo.bounds[i, :greaterthan]) Upper bound: $(blmo.bounds[i, :lessthan]))"
             )
@@ -37,6 +37,81 @@ function is_simple_linear_feasible(sblmo::CubeSimpleBLMO, v)
     return true
 end
 
+function is_simple_inface_feasible(sblmo::CubeSimpleBLMO, a, x, lb, ub, int_vars; kwargs...)
+	return is_simple_inface_feasible_subroutine(sblmo, a, x, lb, ub, int_vars; kwargs)
+end
+
+function is_decomposition_invariant_oracle_simple(sblmo::CubeSimpleBLMO)
+    return true
+end
+
+"""
+If the entry in x is at the boundary, choose the corresponding bound.
+Otherwise, if the entry in direction is positve, choose the lower bound. Else, choose the upper bound.
+"""
+function bounded_compute_inface_extreme_point(sblmo::CubeSimpleBLMO, direction, x, lb, ub, int_vars; kwargs...)
+    n = length(x)
+    v = copy(x)
+    non_fixed_int = findall(lb .!= ub)
+    non_fixed_int_idx = int_vars[non_fixed_int]
+
+    idx = collect(1:n)
+    non_int_idx = setdiff(idx, int_vars)
+    non_fixed_idx = vcat(non_fixed_int_idx, non_int_idx)
+
+     # For non_fixed coordinates, zero-sum means that they are all fixed to origin.
+    sx = sum(x[non_fixed_idx])
+    if sx <= 0
+        return v
+    end
+    # Fix the point to the same face.
+    # Zero will be return only if d_i is greater than zero.
+    for idx in non_fixed_idx
+        if x[idx] > 0 
+            if x[idx] ≈ 1
+                v[idx] = 1
+            else
+                v[idx] = direction[idx] >  0 ? 0 : 1
+            end
+        end
+    end
+    return v       
+end
+
+"""
+Compute the maximum step size for each entry and return the minium of all the possible step sizes.
+"""
+function bounded_dicg_maximum_step(sblmo::CubeSimpleBLMO, direction, x, lb, ub, int_vars; kwargs...)
+    gamma_max = one(eltype(direction))
+    for idx in eachindex(x)
+        di = direction[idx]
+        if idx in int_vars
+            i = findfirst(x -> x == idx, int_vars)
+            if di < 0
+                gamma_max = min(gamma_max, (ub[i] - x[idx]) / -di)
+            elseif di > 0
+                gamma_max = min(gamma_max, (x[idx] - lb[i]) / di)
+            end
+        else
+            if di < 0
+                gamma_max = min(gamma_max, (sblmo.upper_bounds[idx] - x[idx]) / -di)
+            elseif di > 0
+                gamma_max = min(gamma_max, (x[idx] - sblmo.lower_bounds[idx]) / di)
+            end
+        end
+
+    end
+    return gamma_max
+end
+
+function dicg_split_vertices_set_simple(sblmo::CubeSimpleBLMO, x, vidx)
+    x0_left = copy(x)
+    x0_right = copy(x)
+    x0_left[vidx] = floor(x[vidx])
+    x0_right[vidx] = ceil(x[vidx])
+    return x0_left, x0_right
+end
+
 """
     ProbablitySimplexSimpleBLMO(N)
 
@@ -46,6 +121,14 @@ struct ProbabilitySimplexSimpleBLMO <: SimpleBoundableLMO
     N::Float64
 end
 
+function is_decomposition_invariant_oracle_simple(sblmo::ProbabilitySimplexSimpleBLMO)
+    return true  
+end
+
+function is_simple_inface_feasible(sblmo::ProbabilitySimplexSimpleBLMO, a, x, lb, ub, int_vars; kwargs...)
+	return is_simple_inface_feasible_subroutine(sblmo, a, x, lb, ub, int_vars; kwargs)
+end
+
 """
 Assign the largest possible values to the entries corresponding to the smallest entries of d.
 """
@@ -68,6 +151,64 @@ function bounded_compute_extreme_point(sblmo::ProbabilitySimplexSimpleBLMO, d, l
     return v
 end
 
+function bounded_compute_inface_extreme_point(sblmo::ProbabilitySimplexSimpleBLMO, direction, x, lb, ub, int_vars; kwargs...)
+    a = zeros(length(x))
+    if sblmo.N in lb
+        idx = findfirst(x->x==sblmo.N, lb)
+        a[idx] = sblmo.N
+        return a
+    end
+    min_val = Inf
+    min_idx = -1
+    for idx in eachindex(direction)
+        val = direction[idx]
+        if val < min_val && x[idx] > 0
+            min_val = val
+            min_idx = idx
+        end
+    end
+    a[min_idx] = 1.0
+    return a
+end
+
+function bounded_dicg_maximum_step(
+    sblmo::ProbabilitySimplexSimpleBLMO,
+    direction,
+    x,
+    lb,
+    ub,
+    int_vars;
+    tol = 1e-6,
+    kwargs...,
+)
+    # the direction should never violate the simplex constraint because it would correspond to a gamma_max > 1
+    gamma_max = one(eltype(direction))
+    for idx in eachindex(x)
+        di = direction[idx]
+        if di > tol
+            gamma_max = min(gamma_max, (x[idx] - lb[idx]) / di)
+        elseif di < -tol
+            gamma_max = min(gamma_max, (ub[idx] - x[idx]) / -di)
+        end
+
+        if gamma_max == 0.0
+            return 0.0
+        end
+    end
+    return gamma_max
+end
+
+function dicg_split_vertices_set_simple(sblmo::ProbabilitySimplexSimpleBLMO, x, vidx)
+    n = length(x)
+    x0_left = copy(x)
+    sum_val = sum(x) - x[vidx]
+    x0_left .+= (n-1) / sum_val
+    x0_left[vidx] = floor(x[vidx])
+    x0_right = zeros(length(x))
+    x0_right[vidx] = 1.0
+    return x0_left, x0_right
+end
+
 function is_simple_linear_feasible(sblmo::ProbabilitySimplexSimpleBLMO, v)
     if sum(v .≥ 0) < length(v)
         @debug "v has negative entries: $(v)"
@@ -156,6 +297,17 @@ struct UnitSimplexSimpleBLMO <: SimpleBoundableLMO
     N::Float64
 end
 
+function is_decomposition_invariant_oracle_simple(sblmo::UnitSimplexSimpleBLMO)
+    return true  
+end
+
+function is_simple_inface_feasible(sblmo::UnitSimplexSimpleBLMO, a, x, lb, ub, int_vars; kwargs...)
+    if isapprox(sum(x), N; atol = atol, rtol = rtol) && !isapprox(sum(a), N; atol = atol, rtol = rtol)
+        return false
+    end
+    return is_simple_inface_feasible_subroutine(sblmo, a, x, lb, ub, int_vars; kwargs)
+end
+
 """
 For all positive entries of d, assign the corresponding lower bound.
 For non-positive entries, assign largest possible value in increasing order.
@@ -182,6 +334,80 @@ function bounded_compute_extreme_point(sblmo::UnitSimplexSimpleBLMO, d, lb, ub,
     return v
 end
 
+function bounded_compute_inface_extreme_point(sblmo::UnitSimplexSimpleBLMO, direction, x, lb, ub, int_vars; kwargs...)
+    if sblmo.N in lb
+        idx = findfirst(x->x==sblmo.N, lb)
+        a = zeros(length(x))
+        a[idx] = sblmo.N
+        return a
+    end
+        
+     # For non_fixed dimensions, zero-vector x means fixing to all coordinate faces, return zero-vector
+    sx = sum(x)
+    if sx <= 0
+        return zeros(length(x))
+    end
+    
+    min_val = Inf
+    min_idx = -1
+
+    for idx in eachindex(direction)
+        val = direction[idx]
+        if val < min_val && x[idx] > 0
+            min_val = val
+            min_idx = idx
+        end
+    end
+        
+    if sx ≉ sblmo.N && min_val > 0
+        return zeros(length(x))
+    end
+    a = zeros(length(x))
+    a[min_idx] = sblmo.N
+    return a
+end
+
+
+function Boscia.bounded_dicg_maximum_step(
+    sblmo::UnitSimplexSimpleBLMO,
+    direction,
+    x,
+    lb,
+    ub,
+    int_vars;
+    tol = 1e-6,
+    kwargs...,
+)
+    # the direction should never violate the simplex constraint because it would correspond to a gamma_max > 1.
+    gamma_max = one(eltype(direction))
+    for idx in eachindex(x)
+        di = direction[idx]
+        if di > tol
+            gamma_max = min(gamma_max, (x[idx] - lb[idx]) / di)
+        elseif di < -tol
+            gamma_max = min(gamma_max, (ub[idx] - x[idx]) / -di)
+        end
+		
+        if gamma_max == 0.0
+            return 0.0
+        end
+    end
+
+	# the sum of entries should be smaller than N.
+	if sum(direction) < 0.0
+		gamma_max = min(gamma_max, (sum(x)-sblmo.N) / sum(direction))
+	end
+    return gamma_max
+end
+
+function dicg_split_vertices_set_simple(sblmo::UnitSimplexSimpleBLMO, x, vidx)
+    x0_left = copy(x)
+    x0_left[vidx] = floor(x[vidx])
+    x0_right = zeros(length(x))
+    x0_right[vidx] = 1.0
+    return x0_left, x0_right
+end
+
 function is_simple_linear_feasible(sblmo::UnitSimplexSimpleBLMO, v)
     if sum(v .≥ 0) < length(v)
         @debug "v has negative entries: $(v)"
@@ -224,3 +450,23 @@ function rounding_hyperplane_heuristic(tree::Bonobo.BnBTree, tlmo::TimeTrackingL
     end
     return [z], true
 end
+
+function is_simple_inface_feasible_subroutine(sblmo::SimpleBoundableLMO, a, x, lb, ub, int_vars; atol = 1e-6, rtol = 1e-5, kwargs...)
+    for i in eachindex(x)
+        if i in int_vars
+            idx = findfirst(x -> x == i, int_vars)
+            if isapprox(x[idx], lb[idx]; atol = atol, rtol = rtol) && !isapprox(a[i], lb[idx]; atol = atol, rtol = rtol)
+                return false
+            elseif isapprox(x[idx], ub[idx]; atol = atol, rtol = rtol) && !isapprox(a[i], ub[idx]; atol = atol, rtol = rtol)
+                return false
+            end
+        else
+            if isapprox(x[i], sblmo.lower_bounds[i]; atol = atol, rtol = rtol) && !isapprox(a[i], sblmo.lower_bounds[i]; atol = atol, rtol = rtol)
+                return false
+            elseif isapprox(x[i], sblmo.upper_bounds[i]; atol = atol, rtol = rtol) && !isapprox(a[i], sblmo.upper_bounds[i]; atol = atol, rtol = rtol)
+                return false
+            end
+        end
+    end
+    return true
+end
diff --git a/src/time_tracking_lmo.jl b/src/time_tracking_lmo.jl
index d9b50141c..cd8281387 100644
--- a/src/time_tracking_lmo.jl
+++ b/src/time_tracking_lmo.jl
@@ -20,6 +20,40 @@ TimeTrackingLMO(blmo::BoundedLinearMinimizationOracle) =
 TimeTrackingLMO(blmo::BoundedLinearMinimizationOracle, int_vars) =
     TimeTrackingLMO(blmo, Float64[], Int[], Int[], 0, int_vars)
 
+
+is_decomposition_invariant_oracle(tlmo::TimeTrackingLMO) = is_decomposition_invariant_oracle(tlmo.blmo)
+
+function is_inface_feasible(tlmo::TimeTrackingLMO, a, x)
+	return is_inface_feasible(tlmo.blmo, a, x)
+end
+
+function compute_inface_extreme_point(tlmo::TimeTrackingLMO, direction, x; lazy=false, kwargs...)
+    tlmo.ncalls += 1
+    free_model(tlmo.blmo)
+    a = compute_inface_extreme_point(tlmo.blmo, direction, x)
+    
+    if !is_linear_feasible(tlmo, a)
+        @debug "Vertex not linear feasible $(a)"
+        @assert is_linear_feasible(tlmo, a)
+    end
+
+    opt_times, numberofnodes, simplex_iterations = get_BLMO_solve_data(tlmo.blmo)
+
+    push!(tlmo.optimizing_times, opt_times)
+    push!(tlmo.optimizing_nodes, numberofnodes)
+    push!(tlmo.simplex_iterations, simplex_iterations)
+
+    free_model(tlmo.blmo)
+    
+    return a
+end
+
+function dicg_maximum_step(tlmo::TimeTrackingLMO, direction, x)
+    gamma_max = dicg_maximum_step(tlmo.blmo, direction, x)
+    return gamma_max
+end
+
+
 # if we want to reset the info between nodes in Bonobo
 function reset!(tlmo::TimeTrackingLMO)
     empty!(tlmo.optimizing_times)
diff --git a/src/utilities.jl b/src/utilities.jl
index 0bb0041bb..046b306cf 100644
--- a/src/utilities.jl
+++ b/src/utilities.jl
@@ -119,6 +119,59 @@ function split_vertices_set!(
     return (active_set, right_as)
 end
 
+function split_pre_computed_set!(
+    x, 
+    pre_computed_set::Vector, 
+    tree, vidx::Int,
+    local_bounds::IntegerBounds;
+    atol = 1e-5, 
+    rtol = 1e-5, 
+    kwargs...
+)
+    pre_computed_set_left = []
+    pre_computed_set_right = []
+    for atom in pre_computed_set
+        if !is_bound_feasible(local_bounds, atom)
+            continue
+        end
+        if atom[vidx] >= ceil(x[vidx]) || isapprox(atom[vidx], ceil(x[vidx]), atol = atol, rtol = rtol)
+            push!(pre_computed_set_right, atom)
+        elseif atom[vidx] <= floor(x[vidx]) || isapprox(atom[vidx], floor(x[vidx]), atol = atol, rtol = rtol)
+            push!(pre_computed_set_left, atom)
+        end
+    end
+    return pre_computed_set_left, pre_computed_set_right
+end
+
+function build_dicg_start_point(lmo::TimeTrackingLMO)
+	blmo = lmo.blmo
+	n, _ = get_list_of_variables(blmo)
+	d = ones(n)
+	x0 = FrankWolfe.compute_extreme_point(lmo, d)
+	return x0
+end
+
+function dicg_start_point_initialize(lmo::TimeTrackingLMO, active_set::FrankWolfe.ActiveSet{T, R}, pre_computed_set) where {T, R}
+	if lmo.ncalls == 0
+		return FrankWolfe.get_active_set_iterate(active_set)
+	end
+	if pre_computed_set === nothing
+		x0 = build_dicg_start_point(lmo)
+	else
+		if !isempty(pre_computed_set)
+			# We pick a point by averaging the pre_computed_atoms as warm-start.  
+			num_pre_computed_set = length(pre_computed_set)
+			x0 = sum(pre_computed_set) / num_pre_computed_set
+		else
+			# We pick a random point.
+			d = FrankWolfe.get_active_set_iterate(active_set)
+			d = ones(length(d))
+			x0 = FrankWolfe.compute_extreme_point(lmo, d)
+		end
+	end
+end
+
+
 """
 Split a discarded vertices set between left and right children.
 """
diff --git a/test/simple_test_ODWB.jl b/test/simple_test_ODWB.jl
new file mode 100644
index 000000000..c94c2ffc8
--- /dev/null
+++ b/test/simple_test_ODWB.jl
@@ -0,0 +1,151 @@
+############## A optimal design ######################################################################
+# Problem described here: https://link.springer.com/article/10.1007/s11222-014-9476-y
+# "A first-order algorithm for the A-optimal experimental design Problem: a mathematical programming approach"
+
+# min 1/(trace(∑x_i v_iv_i^T))
+# s.t. \sum x_i = s
+#       lb ≤ x ≤ ub
+#       x ∈ Z^m
+
+# v_i ∈ R^n
+# n - number of parameters
+# m - number of possible experiments
+# A = [v_1^T,.., v_m^T], so the rows of A correspond to the different experiments
+
+
+################################ D optimal design ########################################################################
+# Problem described here: https://arxiv.org/pdf/2302.07386.pdf
+# "Branch-and-Bound for D-Optimality with fast local search and bound tightening"
+
+# min log(1/(det(∑x_i v_iv_i^T)))
+# s.t. \sum x_i = s
+#       lb ≤ x ≤ ub
+#       x ∈ Z^m
+
+# v_i ∈ R^n
+# n - number of parameters
+# m - number of possible experiments
+# A = [v_1^T,.., v_m^T], so the rows of A correspond to the different experiments
+
+################################ D-fusion design ########################################################################
+# Problem described here: https://arxiv.org/pdf/2302.07386.pdf
+# "Branch-and-Bound for D-Optimality with fast local search and bound tightening"
+
+# min log(1/(det(∑x_i v_iv_i^T)))
+# s.t. \sum x_i = s
+#       lb ≤ x ≤ ub
+#       x ∈ Z^m
+
+# v_i ∈ R^n
+# n - number of parameters
+# m - number of possible experiments
+# A = [v_1^T,.., v_m^T], so the rows of A correspond to the different experiments
+
+using Boscia
+using FrankWolfe
+using Bonobo
+using Random
+using SCIP
+using JuMP
+using Hypatia
+import Hypatia.Cones: vec_length, vec_copyto!, svec_length, svec_side
+import Hypatia.Cones: smat_to_svec!, svec_to_smat!
+const Cones = Hypatia.Cones
+using Pajarito
+using PajaritoExtras # https://github.com/chriscoey/PajaritoExtras.jl
+using HiGHS
+using LinearAlgebra
+using Statistics
+using Distributions
+import MathOptInterface
+using Printf
+using Dates
+using Test
+using ProfileView
+using DataFrames
+using CSV
+const MOI = MathOptInterface
+const MOIU = MOI.Utilities
+
+import MathOptSetDistances
+const MOD = MathOptSetDistances
+
+include("utilities.jl")
+include("opt_design_boscia.jl")
+include("scip_oa.jl")
+include("opt_design_scip.jl")
+include("spectral_functions_JuMP.jl")
+include("opt_design_pajarito.jl")
+include("opt_design_custom_BB.jl")
+
+seed = rand(UInt64)
+#seed = 0xad2cba3722a98b62
+@show seed
+Random.seed!(seed)
+
+dimensions = [20,30]
+facs = [10,4]
+time_limit = 300
+verbose = false
+
+m = 20
+k = 4
+n = Int(floor(m/k))
+
+A, _, N, ub, _ = build_data(seed, m, n, false, false; scaling_C=false)
+ub = ones(m)
+o = SCIP.Optimizer()
+MOI.set(o, MOI.Silent(), true)
+lmo, x = build_lmo(o, m, N, ub)
+blmo = build_blmo(m, N, ub)
+branching_strategy = Bonobo.MOST_INFEASIBLE()
+heu = Boscia.Heuristic()
+
+result = 0.0
+domain_oracle = build_domain_oracle(A, n)
+f, grad! = build_a_criterion(A, false, μ=1e-4, build_safe=true, long_run=false)
+_, active_set, S = build_start_point2(A, m, n, N, ub)
+z = greedy_incumbent(A, m, n, N, ub)
+#x, _, result = Boscia.solve(f, grad!, lmo; verbose=false, time_limit=10, active_set=active_set, domain_oracle=domain_oracle, start_solution=z, dual_tightening=false, global_dual_tightening=false, lazy_tolerance=2.0, branching_strategy=branching_strategy, use_shadow_set=false, custom_heuristics=[heu]) 
+
+function Boscia.is_decomposition_invariant_oracle_simple(sblmo::Boscia.ProbabilitySimplexSimpleBLMO)
+    return true  
+end
+
+function Boscia.bounded_compute_inface_extreme_point(sblmo::Boscia.ProbabilitySimplexSimpleBLMO, d, x, lb, ub, int_vars; kwargs...)
+    a = zeros(length(d))
+    indices = collect(1:length(d))
+    perm = sortperm(d)
+
+    # The lower bounds always have to be met. 
+    a[int_vars] = lb
+
+    for i in indices[perm]
+        if x[i] !== 0.0
+            if i in int_vars
+                idx = findfirst(x -> x == i, int_vars)
+                a[i] += min(ub[idx] - lb[idx], sblmo.N - sum(a))
+            else
+                a[i] += sblmo.N - sum(a)
+            end
+        end
+    end
+    return a
+end
+
+function Boscia.bounded_dicg_maximum_step(sblmo::Boscia.ProbabilitySimplexSimpleBLMO, direction, x, lb, ub, int_vars; kwargs...)
+    gamma_max = one(eltype(direction))
+    @inbounds for idx in eachindex(x)
+        di = direction[idx]
+        if di < 0
+            gamma_max = min(gamma_max, ub[idx]-x[idx])
+        elseif di > 0
+            gamma_max = min(gamma_max, x[idx]-lb[idx])
+        end
+    end
+    return gamma_max
+end
+
+
+ProfileView.@profview x, _, result = Boscia.solve(f, grad!, blmo; verbose=true, time_limit=time_limit, active_set=active_set, domain_oracle=domain_oracle, start_solution=z, dual_tightening=false, global_dual_tightening=false, lazy_tolerance=2.0, branching_strategy=branching_strategy, use_shadow_set=false, custom_heuristics=[heu])  
+#x, _, result = Boscia.solve(f, grad!, blmo; verbose=true, time_limit=time_limit, active_set=active_set, domain_oracle=domain_oracle, start_solution=z, dual_tightening=false, global_dual_tightening=false, lazy_tolerance=2.0, branching_strategy=branching_strategy, use_shadow_set=false, custom_heuristics=[heu], variant=Boscia.DICG(),lazy=false)
\ No newline at end of file