From 3fc7a605ef02e2fd7c6945927001df112dff8e5f Mon Sep 17 00:00:00 2001
From: Wenjie Xiao <116566662+WenjieXiao-2022@users.noreply.github.com>
Date: Mon, 1 Sep 2025 16:03:06 +0200
Subject: [PATCH 01/21] add Birkhoff_Polytope

---
 src/Birkhoff_polytope.jl | 445 +++++++++++++++++++++++++++++++++++++++
 1 file changed, 445 insertions(+)
 create mode 100644 src/Birkhoff_polytope.jl

diff --git a/src/Birkhoff_polytope.jl b/src/Birkhoff_polytope.jl
new file mode 100644
index 0000000..ebbf32e
--- /dev/null
+++ b/src/Birkhoff_polytope.jl
@@ -0,0 +1,445 @@
+"""
+	BirkhoffBLMO
+
+A bounded LMO for the Birkhoff polytope. This oracle computes an extreme point subject to  
+node-specific bounds on the integer variables.
+"""
+struct BirkhoffBLMO <: Boscia.BoundedLinearMinimizationOracle
+	append_by_column::Bool
+	dim::Int
+	lower_bounds::Vector{Float64}
+	upper_bounds::Vector{Float64}
+	int_vars::Vector{Int}
+	fixed_to_one_rows::Vector{Int}
+	fixed_to_one_cols::Vector{Int}
+	index_map_rows::Vector{Int}
+	index_map_cols::Vector{Int}
+	atol::Float64
+	rtol::Float64
+end
+
+BirkhoffBLMO(dim, lower_bounds, upper_bounds, int_vars; append_by_column = true) =
+	BirkhoffBLMO(append_by_column, dim, lower_bounds, upper_bounds, int_vars, Int[], Int[], collect(1:dim), collect(1:dim), 1e-6, 1e-3)
+
+## Necessary
+
+"""
+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.compute_extreme_point(blmo::BirkhoffBLMO, d; kwargs...)
+    n = blmo.dim
+
+    if size(d, 2) == 1
+        d = blmo.append_by_column ? reshape(d, (n, n)) : transpose(reshape(d, (n, n)))
+    end
+
+    fixed_to_one_rows = blmo.fixed_to_one_rows
+    fixed_to_one_cols = blmo.fixed_to_one_cols
+    index_map_rows = blmo.index_map_rows
+    index_map_cols = blmo.index_map_cols
+	ub = blmo.upper_bounds
+
+    nreduced = length(index_map_rows)
+    type = typeof(d[1, 1])
+    d2 = ones(Union{type,Missing}, nreduced, nreduced)
+    for j in 1:nreduced
+		col_orig = index_map_cols[j]
+        for i in 1:nreduced
+            row_orig = index_map_rows[i]
+            if blmo.append_by_column
+                orig_linear_idx = (col_orig - 1) * n + row_orig
+            else
+                orig_linear_idx = (row_orig - 1) * n + col_orig
+            end
+            # interdict arc when fixed to zero
+            if ub[orig_linear_idx] <= eps()
+                if blmo.append_by_column
+                    d2[i, j] = missing
+                else
+                    d2[j, i] = missing
+                end
+            else
+                if blmo.append_by_column
+                    d2[i, j] = d[row_orig, col_orig]
+                else
+                    d2[j, i] = d[col_orig, row_orig]
+                end
+            end
+        end
+    end
+
+    m = SparseArrays.spzeros(n, n)
+    for (i, j) in zip(fixed_to_one_rows, fixed_to_one_cols)
+        m[i, j] = 1
+    end
+
+    res_mat = Hungarian.munkres(d2)
+    (rows, cols, vals) = SparseArrays.findnz(res_mat)
+    @inbounds for i in eachindex(cols)
+        m[index_map_rows[rows[i]], index_map_cols[cols[i]]] = (vals[i] == 2)
+    end
+
+    m = if blmo.append_by_column
+        # Convert sparse matrix to sparse vector by columns
+        I, J, V = SparseArrays.findnz(m)
+        linear_indices = (J .- 1) .* n .+ I
+        SparseArrays.sparsevec(linear_indices, V, n^2)
+    else
+        # Convert sparse matrix to sparse vector by rows (transpose first)
+        mt = SparseArrays.sparse(LinearAlgebra.transpose(m))
+        I, J, V = SparseArrays.findnz(mt)
+        linear_indices = (J .- 1) .* n .+ I
+        SparseArrays.sparsevec(linear_indices, V, n^2)
+    end
+    return m
+end
+
+"""
+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.compute_inface_extreme_point(
+	blmo::BirkhoffBLMO,
+	direction,
+	x;
+	kwargs...,
+)
+	n = blmo.dim
+
+	if size(direction, 2) == 1
+		direction =
+			blmo.append_by_column ? reshape(direction, (n, n)) :
+			transpose(reshape(direction, (n, n)))
+	end
+
+	if size(x, 2) == 1
+		x = blmo.append_by_column ? reshape(x, (n, n)) : transpose(reshape(x, (n, n)))
+	end
+
+
+	fixed_to_one_rows = copy(blmo.fixed_to_one_rows)
+	fixed_to_one_cols = copy(blmo.fixed_to_one_cols)
+	index_map_rows = copy(blmo.index_map_rows)
+	index_map_cols = copy(blmo.index_map_cols)
+	ub = blmo.upper_bounds
+
+	nreduced = length(blmo.index_map_rows)
+
+	delete_index_map_rows = Int[]
+	delete_index_map_cols = Int[]
+	for j in 1:nreduced
+		for i in 1:nreduced
+			row_orig = index_map_rows[i]
+			col_orig = index_map_cols[j]
+			if x[row_orig, col_orig] >= 1-eps()
+				push!(fixed_to_one_rows, row_orig)
+				push!(fixed_to_one_cols, col_orig)
+
+				push!(delete_index_map_rows, i)
+				push!(delete_index_map_cols, j)
+			end
+		end
+	end
+
+	unique!(delete_index_map_rows)
+	unique!(delete_index_map_cols)
+	sort!(delete_index_map_rows)
+	sort!(delete_index_map_cols)
+	deleteat!(index_map_rows, delete_index_map_rows)
+	deleteat!(index_map_cols, delete_index_map_cols)
+
+	nreduced = length(index_map_rows)
+	type = typeof(direction[1, 1])
+	d2 = ones(Union{type, Missing}, nreduced, nreduced)
+	for j in 1:nreduced
+		for i in 1:nreduced
+			row_orig = index_map_rows[i]
+			col_orig = index_map_cols[j]
+			if blmo.append_by_column
+				orig_linear_idx = (col_orig-1)*n+row_orig
+			else
+				orig_linear_idx = (row_orig-1)*n+col_orig
+			end
+			# interdict arc when fixed to zero
+			if ub[orig_linear_idx] <= eps() || x[row_orig, col_orig] <= eps()
+				if blmo.append_by_column
+					d2[i, j] = missing
+				else
+					d2[j, i] = missing
+				end
+			else
+				if blmo.append_by_column
+					d2[i, j] = direction[row_orig, col_orig]
+				else
+					d2[j, i] = direction[col_orig, row_orig]
+				end
+			end
+		end
+	end
+
+	m = SparseArrays.spzeros(n, n)
+	for (i, j) in zip(fixed_to_one_rows, fixed_to_one_cols)
+		m[i, j] = 1
+	end
+
+	res_mat = Hungarian.munkres(d2)
+	(rows, cols, vals) = SparseArrays.findnz(res_mat)
+	@inbounds for i in eachindex(cols)
+		m[index_map_rows[rows[i]], index_map_cols[cols[i]]] = (vals[i] == 2)
+	end
+
+	m = if blmo.append_by_column
+		# Convert sparse matrix to sparse vector by columns
+		I, J, V = SparseArrays.findnz(m)
+		linear_indices = (J .- 1) .* n .+ I
+		SparseArrays.sparsevec(linear_indices, V, n^2)
+	else
+		# Convert sparse matrix to sparse vector by rows (transpose first)
+		mt = SparseArrays.sparse(LinearAlgebra.transpose(m))
+		I, J, V = SparseArrays.findnz(mt)
+		linear_indices = (J .- 1) .* n .+ I
+		SparseArrays.sparsevec(linear_indices, V, n^2)
+	end
+
+	return m
+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.dicg_maximum_step(blmo::BirkhoffBLMO, direction, x; kwargs...)
+	n = blmo.dim
+
+	direction =
+		blmo.append_by_column ? reshape(direction, (n, n)) : transpose(reshape(direction, (n, n)))
+	x = blmo.append_by_column ? reshape(x, (n, n)) : transpose(reshape(x, (n, n)))
+	return FrankWolfe.dicg_maximum_step(FrankWolfe.BirkhoffPolytopeLMO(), direction, x)
+end
+
+function Boscia.is_decomposition_invariant_oracle(blmo::BirkhoffBLMO)
+	return true
+end
+
+"""
+The sum of each row and column has to be equal to 1.
+"""
+function Boscia.is_linear_feasible(blmo::BirkhoffBLMO, v::AbstractVector)
+	n = blmo.dim
+	for i in 1:n
+		# append by column ? column sum : row sum 
+		if !isapprox(sum(v[((i-1)*n+1):(i*n)]), 1.0, atol = 1e-6, rtol = 1e-3)
+			@debug "Column sum not 1: $(sum(v[((i-1)*n+1):(i*n)]))"
+			return false
+		end
+		# append by column ? row sum : column sum
+		if !isapprox(sum(v[i: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
+
+# Read global bounds from the problem.
+function Boscia.build_global_bounds(blmo::BirkhoffBLMO, integer_variables)
+	global_bounds = Boscia.IntegerBounds()
+	for (idx, int_var) in enumerate(blmo.int_vars)
+		push!(global_bounds, (int_var, blmo.lower_bounds[idx]), :greaterthan)
+		push!(global_bounds, (int_var, blmo.upper_bounds[idx]), :lessthan)
+	end
+	return global_bounds
+end
+
+# Get list of variables indices. 
+# If the problem has n variables, they are expected to contiguous and ordered from 1 to n.
+function Boscia.get_list_of_variables(blmo::BirkhoffBLMO)
+	n = blmo.dim^2
+	return n, collect(1:n)
+end
+
+# Get list of integer variables
+function Boscia.get_integer_variables(blmo::BirkhoffBLMO)
+	return blmo.int_vars
+end
+
+# Get the index of the integer variable the bound is working on.
+function Boscia.get_int_var(blmo::BirkhoffBLMO, cidx)
+	return blmo.int_vars[cidx]
+end
+
+# Get the list of lower bounds.
+function Boscia.get_lower_bound_list(blmo::BirkhoffBLMO)
+	return collect(1:length(blmo.lower_bounds))
+end
+
+# Get the list of upper bounds.
+function Boscia.get_upper_bound_list(blmo::BirkhoffBLMO)
+	return collect(1:length(blmo.upper_bounds))
+end
+
+# Read bound value for c_idx.
+function Boscia.get_bound(blmo::BirkhoffBLMO, c_idx, sense::Symbol)
+	if sense == :lessthan
+		return blmo.upper_bounds[c_idx]
+	elseif sense == :greaterthan
+		return blmo.lower_bounds[c_idx]
+	else
+		error("Allowed value for sense are :lessthan and :greaterthan!")
+	end
+end
+
+## Changing the bounds constraints.
+
+# Change the value of the bound c_idx.
+function Boscia.set_bound!(blmo::BirkhoffBLMO, c_idx, value, sense::Symbol)
+	if sense == :greaterthan
+		blmo.lower_bounds[c_idx] = value
+	elseif sense == :lessthan
+		blmo.upper_bounds[c_idx] = value
+	else
+		error("Allowed values for sense are :lessthan and :greaterthan.")
+	end
+end
+
+# Delete bounds.
+function Boscia.delete_bounds!(blmo::BirkhoffBLMO, cons_delete)
+	for (d_idx, sense) in cons_delete
+		if sense == :greaterthan
+			blmo.lower_bounds[d_idx] = -Inf
+		else
+			blmo.upper_bounds[d_idx] = Inf
+		end
+	end
+end
+
+# Add bound constraint.
+function Boscia.add_bound_constraint!(blmo::BirkhoffBLMO, key, value, sense::Symbol)
+	idx = findfirst(x -> x == key, blmo.int_vars)
+	if sense == :greaterthan
+		blmo.lower_bounds[idx] = value
+	elseif sense == :lessthan
+		blmo.upper_bounds[idx] = value
+	else
+		error("Allowed value of sense are :lessthan and :greaterthan!")
+	end
+end
+
+## Checks
+
+# Check if the subject of the bound c_idx is an integer variable (recorded in int_vars).
+function Boscia.is_constraint_on_int_var(blmo::BirkhoffBLMO, c_idx, int_vars)
+	return blmo.int_vars[c_idx] in int_vars
+end
+
+# To check if there is bound for the variable in the global or node bounds.
+function Boscia.is_bound_in(blmo::BirkhoffBLMO, c_idx, bounds)
+	return haskey(bounds, blmo.int_vars[c_idx])
+end
+
+# Has variable an integer constraint?
+function Boscia.has_integer_constraint(blmo::BirkhoffBLMO, idx)
+	return idx in blmo.int_vars
+end
+
+## Optional
+
+function Boscia.check_feasibility(blmo::BirkhoffBLMO)
+	for (lb, ub) in zip(blmo.lower_bounds, blmo.upper_bounds)
+		if ub < lb
+			return Boscia.INFEASIBLE
+		end
+	end
+	# For double stochastic matrices, each row and column must sum to 1
+	# We check if the bounds allow for feasible assignments
+	n0 = blmo.dim
+	n = n0^2
+	# Initialize row and column bound tracking
+	row_min_sum = zeros(n)  # minimum possible sum for each row
+	row_max_sum = zeros(n)  # maximum possible sum for each row
+	col_min_sum = zeros(n)  # minimum possible sum for each column
+	col_max_sum = zeros(n)  # maximum possible sum for each column
+
+	# Process each integer variable
+	for idx in eachindex(blmo.int_vars)
+		var_idx = blmo.int_vars[idx]
+
+		# Convert linear index to (row, col) based on storage format
+		if blmo.append_by_column
+			j = ceil(Int, var_idx / n0)  # column index
+			i = Int(var_idx - n0 * (j - 1))  # row index
+		else
+			i = ceil(Int, var_idx / n0)  # row index  
+			j = Int(var_idx - n0 * (i - 1))  # column index
+		end
+
+		# Add bounds to row and column sums
+		row_min_sum[i] += blmo.lower_bounds[idx]
+		row_max_sum[i] += blmo.upper_bounds[idx]
+		col_min_sum[j] += blmo.lower_bounds[idx]
+		col_max_sum[j] += blmo.upper_bounds[idx]
+	end
+
+	# Check feasibility: each row and column must be able to sum to exactly 1
+	for i in 1:n0
+		# Check row sum constraints
+		if row_min_sum[i] > 1 + eps() || row_max_sum[i] < 1 - eps()
+			return Boscia.INFEASIBLE
+		end
+
+		# Check column sum constraints  
+		if col_min_sum[i] > 1 + eps() || col_max_sum[i] < 1 - eps()
+			return Boscia.INFEASIBLE
+		end
+	end
+
+	return Boscia.OPTIMAL
+end
+
+function Boscia.update_model(blmo::BirkhoffBLMO)
+	dim = blmo.dim
+	n = dim^2
+
+	# fixed_to_one_vars = [v for v in blmo.int_vars if blmo.lower_bounds[v] == 1.0]
+
+	fixed_to_one_vars = blmo.int_vars[blmo.lower_bounds .== 1.0]
+	empty!(blmo.fixed_to_one_rows)
+	empty!(blmo.fixed_to_one_cols)
+	for fixed_to_one_var in fixed_to_one_vars
+		q = (fixed_to_one_var - 1) ÷ dim
+		r = (fixed_to_one_var - 1) - q * dim
+		if blmo.append_by_column
+			i = r + 1
+			j = q + 1
+		else
+			i = q + 1
+			j = r + 1
+		end
+		push!(blmo.fixed_to_one_rows, i)
+		push!(blmo.fixed_to_one_cols, j)
+	end
+
+	nfixed = length(blmo.fixed_to_one_rows)
+	nreduced = dim - nfixed
+
+	# stores the indices of the original matrix that are still in the reduced matrix
+	index_map_rows = fill(1, nreduced)
+	index_map_cols = fill(1, nreduced)
+	idx_in_map_row = 1
+	idx_in_map_col = 1
+	for orig_idx in 1:dim
+		if orig_idx ∉ blmo.fixed_to_one_rows
+			index_map_rows[idx_in_map_row] = orig_idx
+			idx_in_map_row += 1
+		end
+		if orig_idx ∉ blmo.fixed_to_one_cols
+			index_map_cols[idx_in_map_col] = orig_idx
+			idx_in_map_col += 1
+		end
+	end
+
+	empty!(blmo.index_map_rows)
+	empty!(blmo.index_map_cols)
+	append!(blmo.index_map_rows, index_map_rows)
+	append!(blmo.index_map_cols, index_map_cols)
+end

From 903d3a4810cd504b2dcce32b66b589744b580a4f Mon Sep 17 00:00:00 2001
From: Wenjie Xiao <116566662+WenjieXiao-2022@users.noreply.github.com>
Date: Mon, 1 Sep 2025 16:04:32 +0200
Subject: [PATCH 02/21] format

---
 src/Birkhoff_polytope.jl | 564 ++++++++++++++++++++-------------------
 1 file changed, 286 insertions(+), 278 deletions(-)

diff --git a/src/Birkhoff_polytope.jl b/src/Birkhoff_polytope.jl
index ebbf32e..a34a324 100644
--- a/src/Birkhoff_polytope.jl
+++ b/src/Birkhoff_polytope.jl
@@ -5,21 +5,33 @@ A bounded LMO for the Birkhoff polytope. This oracle computes an extreme point s
 node-specific bounds on the integer variables.
 """
 struct BirkhoffBLMO <: Boscia.BoundedLinearMinimizationOracle
-	append_by_column::Bool
-	dim::Int
-	lower_bounds::Vector{Float64}
-	upper_bounds::Vector{Float64}
-	int_vars::Vector{Int}
-	fixed_to_one_rows::Vector{Int}
-	fixed_to_one_cols::Vector{Int}
-	index_map_rows::Vector{Int}
-	index_map_cols::Vector{Int}
-	atol::Float64
-	rtol::Float64
+    append_by_column::Bool
+    dim::Int
+    lower_bounds::Vector{Float64}
+    upper_bounds::Vector{Float64}
+    int_vars::Vector{Int}
+    fixed_to_one_rows::Vector{Int}
+    fixed_to_one_cols::Vector{Int}
+    index_map_rows::Vector{Int}
+    index_map_cols::Vector{Int}
+    atol::Float64
+    rtol::Float64
 end
 
 BirkhoffBLMO(dim, lower_bounds, upper_bounds, int_vars; append_by_column = true) =
-	BirkhoffBLMO(append_by_column, dim, lower_bounds, upper_bounds, int_vars, Int[], Int[], collect(1:dim), collect(1:dim), 1e-6, 1e-3)
+    BirkhoffBLMO(
+        append_by_column,
+        dim,
+        lower_bounds,
+        upper_bounds,
+        int_vars,
+        Int[],
+        Int[],
+        collect(1:dim),
+        collect(1:dim),
+        1e-6,
+        1e-3,
+    )
 
 ## Necessary
 
@@ -38,14 +50,14 @@ function Boscia.compute_extreme_point(blmo::BirkhoffBLMO, d; kwargs...)
     fixed_to_one_cols = blmo.fixed_to_one_cols
     index_map_rows = blmo.index_map_rows
     index_map_cols = blmo.index_map_cols
-	ub = blmo.upper_bounds
+    ub = blmo.upper_bounds
 
     nreduced = length(index_map_rows)
     type = typeof(d[1, 1])
     d2 = ones(Union{type,Missing}, nreduced, nreduced)
-    for j in 1:nreduced
-		col_orig = index_map_cols[j]
-        for i in 1:nreduced
+    for j = 1:nreduced
+        col_orig = index_map_cols[j]
+        for i = 1:nreduced
             row_orig = index_map_rows[i]
             if blmo.append_by_column
                 orig_linear_idx = (col_orig - 1) * n + row_orig
@@ -98,110 +110,105 @@ end
 """
 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.compute_inface_extreme_point(
-	blmo::BirkhoffBLMO,
-	direction,
-	x;
-	kwargs...,
-)
-	n = blmo.dim
-
-	if size(direction, 2) == 1
-		direction =
-			blmo.append_by_column ? reshape(direction, (n, n)) :
-			transpose(reshape(direction, (n, n)))
-	end
-
-	if size(x, 2) == 1
-		x = blmo.append_by_column ? reshape(x, (n, n)) : transpose(reshape(x, (n, n)))
-	end
-
-
-	fixed_to_one_rows = copy(blmo.fixed_to_one_rows)
-	fixed_to_one_cols = copy(blmo.fixed_to_one_cols)
-	index_map_rows = copy(blmo.index_map_rows)
-	index_map_cols = copy(blmo.index_map_cols)
-	ub = blmo.upper_bounds
-
-	nreduced = length(blmo.index_map_rows)
-
-	delete_index_map_rows = Int[]
-	delete_index_map_cols = Int[]
-	for j in 1:nreduced
-		for i in 1:nreduced
-			row_orig = index_map_rows[i]
-			col_orig = index_map_cols[j]
-			if x[row_orig, col_orig] >= 1-eps()
-				push!(fixed_to_one_rows, row_orig)
-				push!(fixed_to_one_cols, col_orig)
-
-				push!(delete_index_map_rows, i)
-				push!(delete_index_map_cols, j)
-			end
-		end
-	end
-
-	unique!(delete_index_map_rows)
-	unique!(delete_index_map_cols)
-	sort!(delete_index_map_rows)
-	sort!(delete_index_map_cols)
-	deleteat!(index_map_rows, delete_index_map_rows)
-	deleteat!(index_map_cols, delete_index_map_cols)
-
-	nreduced = length(index_map_rows)
-	type = typeof(direction[1, 1])
-	d2 = ones(Union{type, Missing}, nreduced, nreduced)
-	for j in 1:nreduced
-		for i in 1:nreduced
-			row_orig = index_map_rows[i]
-			col_orig = index_map_cols[j]
-			if blmo.append_by_column
-				orig_linear_idx = (col_orig-1)*n+row_orig
-			else
-				orig_linear_idx = (row_orig-1)*n+col_orig
-			end
-			# interdict arc when fixed to zero
-			if ub[orig_linear_idx] <= eps() || x[row_orig, col_orig] <= eps()
-				if blmo.append_by_column
-					d2[i, j] = missing
-				else
-					d2[j, i] = missing
-				end
-			else
-				if blmo.append_by_column
-					d2[i, j] = direction[row_orig, col_orig]
-				else
-					d2[j, i] = direction[col_orig, row_orig]
-				end
-			end
-		end
-	end
-
-	m = SparseArrays.spzeros(n, n)
-	for (i, j) in zip(fixed_to_one_rows, fixed_to_one_cols)
-		m[i, j] = 1
-	end
-
-	res_mat = Hungarian.munkres(d2)
-	(rows, cols, vals) = SparseArrays.findnz(res_mat)
-	@inbounds for i in eachindex(cols)
-		m[index_map_rows[rows[i]], index_map_cols[cols[i]]] = (vals[i] == 2)
-	end
-
-	m = if blmo.append_by_column
-		# Convert sparse matrix to sparse vector by columns
-		I, J, V = SparseArrays.findnz(m)
-		linear_indices = (J .- 1) .* n .+ I
-		SparseArrays.sparsevec(linear_indices, V, n^2)
-	else
-		# Convert sparse matrix to sparse vector by rows (transpose first)
-		mt = SparseArrays.sparse(LinearAlgebra.transpose(m))
-		I, J, V = SparseArrays.findnz(mt)
-		linear_indices = (J .- 1) .* n .+ I
-		SparseArrays.sparsevec(linear_indices, V, n^2)
-	end
-
-	return m
+function Boscia.compute_inface_extreme_point(blmo::BirkhoffBLMO, direction, x; kwargs...)
+    n = blmo.dim
+
+    if size(direction, 2) == 1
+        direction =
+            blmo.append_by_column ? reshape(direction, (n, n)) :
+            transpose(reshape(direction, (n, n)))
+    end
+
+    if size(x, 2) == 1
+        x = blmo.append_by_column ? reshape(x, (n, n)) : transpose(reshape(x, (n, n)))
+    end
+
+
+    fixed_to_one_rows = copy(blmo.fixed_to_one_rows)
+    fixed_to_one_cols = copy(blmo.fixed_to_one_cols)
+    index_map_rows = copy(blmo.index_map_rows)
+    index_map_cols = copy(blmo.index_map_cols)
+    ub = blmo.upper_bounds
+
+    nreduced = length(blmo.index_map_rows)
+
+    delete_index_map_rows = Int[]
+    delete_index_map_cols = Int[]
+    for j = 1:nreduced
+        for i = 1:nreduced
+            row_orig = index_map_rows[i]
+            col_orig = index_map_cols[j]
+            if x[row_orig, col_orig] >= 1 - eps()
+                push!(fixed_to_one_rows, row_orig)
+                push!(fixed_to_one_cols, col_orig)
+
+                push!(delete_index_map_rows, i)
+                push!(delete_index_map_cols, j)
+            end
+        end
+    end
+
+    unique!(delete_index_map_rows)
+    unique!(delete_index_map_cols)
+    sort!(delete_index_map_rows)
+    sort!(delete_index_map_cols)
+    deleteat!(index_map_rows, delete_index_map_rows)
+    deleteat!(index_map_cols, delete_index_map_cols)
+
+    nreduced = length(index_map_rows)
+    type = typeof(direction[1, 1])
+    d2 = ones(Union{type,Missing}, nreduced, nreduced)
+    for j = 1:nreduced
+        for i = 1:nreduced
+            row_orig = index_map_rows[i]
+            col_orig = index_map_cols[j]
+            if blmo.append_by_column
+                orig_linear_idx = (col_orig - 1) * n + row_orig
+            else
+                orig_linear_idx = (row_orig - 1) * n + col_orig
+            end
+            # interdict arc when fixed to zero
+            if ub[orig_linear_idx] <= eps() || x[row_orig, col_orig] <= eps()
+                if blmo.append_by_column
+                    d2[i, j] = missing
+                else
+                    d2[j, i] = missing
+                end
+            else
+                if blmo.append_by_column
+                    d2[i, j] = direction[row_orig, col_orig]
+                else
+                    d2[j, i] = direction[col_orig, row_orig]
+                end
+            end
+        end
+    end
+
+    m = SparseArrays.spzeros(n, n)
+    for (i, j) in zip(fixed_to_one_rows, fixed_to_one_cols)
+        m[i, j] = 1
+    end
+
+    res_mat = Hungarian.munkres(d2)
+    (rows, cols, vals) = SparseArrays.findnz(res_mat)
+    @inbounds for i in eachindex(cols)
+        m[index_map_rows[rows[i]], index_map_cols[cols[i]]] = (vals[i] == 2)
+    end
+
+    m = if blmo.append_by_column
+        # Convert sparse matrix to sparse vector by columns
+        I, J, V = SparseArrays.findnz(m)
+        linear_indices = (J .- 1) .* n .+ I
+        SparseArrays.sparsevec(linear_indices, V, n^2)
+    else
+        # Convert sparse matrix to sparse vector by rows (transpose first)
+        mt = SparseArrays.sparse(LinearAlgebra.transpose(m))
+        I, J, V = SparseArrays.findnz(mt)
+        linear_indices = (J .- 1) .* n .+ I
+        SparseArrays.sparsevec(linear_indices, V, n^2)
+    end
+
+    return m
 end
 
 """
@@ -209,237 +216,238 @@ LMO-like operation which computes a vertex minimizing in `direction` on the face
 Fixings are maintained by the oracle (or deduced from `x` itself).
 """
 function Boscia.dicg_maximum_step(blmo::BirkhoffBLMO, direction, x; kwargs...)
-	n = blmo.dim
+    n = blmo.dim
 
-	direction =
-		blmo.append_by_column ? reshape(direction, (n, n)) : transpose(reshape(direction, (n, n)))
-	x = blmo.append_by_column ? reshape(x, (n, n)) : transpose(reshape(x, (n, n)))
-	return FrankWolfe.dicg_maximum_step(FrankWolfe.BirkhoffPolytopeLMO(), direction, x)
+    direction =
+        blmo.append_by_column ? reshape(direction, (n, n)) :
+        transpose(reshape(direction, (n, n)))
+    x = blmo.append_by_column ? reshape(x, (n, n)) : transpose(reshape(x, (n, n)))
+    return FrankWolfe.dicg_maximum_step(FrankWolfe.BirkhoffPolytopeLMO(), direction, x)
 end
 
 function Boscia.is_decomposition_invariant_oracle(blmo::BirkhoffBLMO)
-	return true
+    return true
 end
 
 """
 The sum of each row and column has to be equal to 1.
 """
 function Boscia.is_linear_feasible(blmo::BirkhoffBLMO, v::AbstractVector)
-	n = blmo.dim
-	for i in 1:n
-		# append by column ? column sum : row sum 
-		if !isapprox(sum(v[((i-1)*n+1):(i*n)]), 1.0, atol = 1e-6, rtol = 1e-3)
-			@debug "Column sum not 1: $(sum(v[((i-1)*n+1):(i*n)]))"
-			return false
-		end
-		# append by column ? row sum : column sum
-		if !isapprox(sum(v[i: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
+    n = blmo.dim
+    for i = 1:n
+        # append by column ? column sum : row sum 
+        if !isapprox(sum(v[((i-1)*n+1):(i*n)]), 1.0, atol = 1e-6, rtol = 1e-3)
+            @debug "Column sum not 1: $(sum(v[((i-1)*n+1):(i*n)]))"
+            return false
+        end
+        # append by column ? row sum : column sum
+        if !isapprox(sum(v[i: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
 
 # Read global bounds from the problem.
 function Boscia.build_global_bounds(blmo::BirkhoffBLMO, integer_variables)
-	global_bounds = Boscia.IntegerBounds()
-	for (idx, int_var) in enumerate(blmo.int_vars)
-		push!(global_bounds, (int_var, blmo.lower_bounds[idx]), :greaterthan)
-		push!(global_bounds, (int_var, blmo.upper_bounds[idx]), :lessthan)
-	end
-	return global_bounds
+    global_bounds = Boscia.IntegerBounds()
+    for (idx, int_var) in enumerate(blmo.int_vars)
+        push!(global_bounds, (int_var, blmo.lower_bounds[idx]), :greaterthan)
+        push!(global_bounds, (int_var, blmo.upper_bounds[idx]), :lessthan)
+    end
+    return global_bounds
 end
 
 # Get list of variables indices. 
 # If the problem has n variables, they are expected to contiguous and ordered from 1 to n.
 function Boscia.get_list_of_variables(blmo::BirkhoffBLMO)
-	n = blmo.dim^2
-	return n, collect(1:n)
+    n = blmo.dim^2
+    return n, collect(1:n)
 end
 
 # Get list of integer variables
 function Boscia.get_integer_variables(blmo::BirkhoffBLMO)
-	return blmo.int_vars
+    return blmo.int_vars
 end
 
 # Get the index of the integer variable the bound is working on.
 function Boscia.get_int_var(blmo::BirkhoffBLMO, cidx)
-	return blmo.int_vars[cidx]
+    return blmo.int_vars[cidx]
 end
 
 # Get the list of lower bounds.
 function Boscia.get_lower_bound_list(blmo::BirkhoffBLMO)
-	return collect(1:length(blmo.lower_bounds))
+    return collect(1:length(blmo.lower_bounds))
 end
 
 # Get the list of upper bounds.
 function Boscia.get_upper_bound_list(blmo::BirkhoffBLMO)
-	return collect(1:length(blmo.upper_bounds))
+    return collect(1:length(blmo.upper_bounds))
 end
 
 # Read bound value for c_idx.
 function Boscia.get_bound(blmo::BirkhoffBLMO, c_idx, sense::Symbol)
-	if sense == :lessthan
-		return blmo.upper_bounds[c_idx]
-	elseif sense == :greaterthan
-		return blmo.lower_bounds[c_idx]
-	else
-		error("Allowed value for sense are :lessthan and :greaterthan!")
-	end
+    if sense == :lessthan
+        return blmo.upper_bounds[c_idx]
+    elseif sense == :greaterthan
+        return blmo.lower_bounds[c_idx]
+    else
+        error("Allowed value for sense are :lessthan and :greaterthan!")
+    end
 end
 
 ## Changing the bounds constraints.
 
 # Change the value of the bound c_idx.
 function Boscia.set_bound!(blmo::BirkhoffBLMO, c_idx, value, sense::Symbol)
-	if sense == :greaterthan
-		blmo.lower_bounds[c_idx] = value
-	elseif sense == :lessthan
-		blmo.upper_bounds[c_idx] = value
-	else
-		error("Allowed values for sense are :lessthan and :greaterthan.")
-	end
+    if sense == :greaterthan
+        blmo.lower_bounds[c_idx] = value
+    elseif sense == :lessthan
+        blmo.upper_bounds[c_idx] = value
+    else
+        error("Allowed values for sense are :lessthan and :greaterthan.")
+    end
 end
 
 # Delete bounds.
 function Boscia.delete_bounds!(blmo::BirkhoffBLMO, cons_delete)
-	for (d_idx, sense) in cons_delete
-		if sense == :greaterthan
-			blmo.lower_bounds[d_idx] = -Inf
-		else
-			blmo.upper_bounds[d_idx] = Inf
-		end
-	end
+    for (d_idx, sense) in cons_delete
+        if sense == :greaterthan
+            blmo.lower_bounds[d_idx] = -Inf
+        else
+            blmo.upper_bounds[d_idx] = Inf
+        end
+    end
 end
 
 # Add bound constraint.
 function Boscia.add_bound_constraint!(blmo::BirkhoffBLMO, key, value, sense::Symbol)
-	idx = findfirst(x -> x == key, blmo.int_vars)
-	if sense == :greaterthan
-		blmo.lower_bounds[idx] = value
-	elseif sense == :lessthan
-		blmo.upper_bounds[idx] = value
-	else
-		error("Allowed value of sense are :lessthan and :greaterthan!")
-	end
+    idx = findfirst(x -> x == key, blmo.int_vars)
+    if sense == :greaterthan
+        blmo.lower_bounds[idx] = value
+    elseif sense == :lessthan
+        blmo.upper_bounds[idx] = value
+    else
+        error("Allowed value of sense are :lessthan and :greaterthan!")
+    end
 end
 
 ## Checks
 
 # Check if the subject of the bound c_idx is an integer variable (recorded in int_vars).
 function Boscia.is_constraint_on_int_var(blmo::BirkhoffBLMO, c_idx, int_vars)
-	return blmo.int_vars[c_idx] in int_vars
+    return blmo.int_vars[c_idx] in int_vars
 end
 
 # To check if there is bound for the variable in the global or node bounds.
 function Boscia.is_bound_in(blmo::BirkhoffBLMO, c_idx, bounds)
-	return haskey(bounds, blmo.int_vars[c_idx])
+    return haskey(bounds, blmo.int_vars[c_idx])
 end
 
 # Has variable an integer constraint?
 function Boscia.has_integer_constraint(blmo::BirkhoffBLMO, idx)
-	return idx in blmo.int_vars
+    return idx in blmo.int_vars
 end
 
 ## Optional
 
 function Boscia.check_feasibility(blmo::BirkhoffBLMO)
-	for (lb, ub) in zip(blmo.lower_bounds, blmo.upper_bounds)
-		if ub < lb
-			return Boscia.INFEASIBLE
-		end
-	end
-	# For double stochastic matrices, each row and column must sum to 1
-	# We check if the bounds allow for feasible assignments
-	n0 = blmo.dim
-	n = n0^2
-	# Initialize row and column bound tracking
-	row_min_sum = zeros(n)  # minimum possible sum for each row
-	row_max_sum = zeros(n)  # maximum possible sum for each row
-	col_min_sum = zeros(n)  # minimum possible sum for each column
-	col_max_sum = zeros(n)  # maximum possible sum for each column
-
-	# Process each integer variable
-	for idx in eachindex(blmo.int_vars)
-		var_idx = blmo.int_vars[idx]
-
-		# Convert linear index to (row, col) based on storage format
-		if blmo.append_by_column
-			j = ceil(Int, var_idx / n0)  # column index
-			i = Int(var_idx - n0 * (j - 1))  # row index
-		else
-			i = ceil(Int, var_idx / n0)  # row index  
-			j = Int(var_idx - n0 * (i - 1))  # column index
-		end
-
-		# Add bounds to row and column sums
-		row_min_sum[i] += blmo.lower_bounds[idx]
-		row_max_sum[i] += blmo.upper_bounds[idx]
-		col_min_sum[j] += blmo.lower_bounds[idx]
-		col_max_sum[j] += blmo.upper_bounds[idx]
-	end
-
-	# Check feasibility: each row and column must be able to sum to exactly 1
-	for i in 1:n0
-		# Check row sum constraints
-		if row_min_sum[i] > 1 + eps() || row_max_sum[i] < 1 - eps()
-			return Boscia.INFEASIBLE
-		end
-
-		# Check column sum constraints  
-		if col_min_sum[i] > 1 + eps() || col_max_sum[i] < 1 - eps()
-			return Boscia.INFEASIBLE
-		end
-	end
-
-	return Boscia.OPTIMAL
+    for (lb, ub) in zip(blmo.lower_bounds, blmo.upper_bounds)
+        if ub < lb
+            return Boscia.INFEASIBLE
+        end
+    end
+    # For double stochastic matrices, each row and column must sum to 1
+    # We check if the bounds allow for feasible assignments
+    n0 = blmo.dim
+    n = n0^2
+    # Initialize row and column bound tracking
+    row_min_sum = zeros(n)  # minimum possible sum for each row
+    row_max_sum = zeros(n)  # maximum possible sum for each row
+    col_min_sum = zeros(n)  # minimum possible sum for each column
+    col_max_sum = zeros(n)  # maximum possible sum for each column
+
+    # Process each integer variable
+    for idx in eachindex(blmo.int_vars)
+        var_idx = blmo.int_vars[idx]
+
+        # Convert linear index to (row, col) based on storage format
+        if blmo.append_by_column
+            j = ceil(Int, var_idx / n0)  # column index
+            i = Int(var_idx - n0 * (j - 1))  # row index
+        else
+            i = ceil(Int, var_idx / n0)  # row index  
+            j = Int(var_idx - n0 * (i - 1))  # column index
+        end
+
+        # Add bounds to row and column sums
+        row_min_sum[i] += blmo.lower_bounds[idx]
+        row_max_sum[i] += blmo.upper_bounds[idx]
+        col_min_sum[j] += blmo.lower_bounds[idx]
+        col_max_sum[j] += blmo.upper_bounds[idx]
+    end
+
+    # Check feasibility: each row and column must be able to sum to exactly 1
+    for i = 1:n0
+        # Check row sum constraints
+        if row_min_sum[i] > 1 + eps() || row_max_sum[i] < 1 - eps()
+            return Boscia.INFEASIBLE
+        end
+
+        # Check column sum constraints  
+        if col_min_sum[i] > 1 + eps() || col_max_sum[i] < 1 - eps()
+            return Boscia.INFEASIBLE
+        end
+    end
+
+    return Boscia.OPTIMAL
 end
 
 function Boscia.update_model(blmo::BirkhoffBLMO)
-	dim = blmo.dim
-	n = dim^2
-
-	# fixed_to_one_vars = [v for v in blmo.int_vars if blmo.lower_bounds[v] == 1.0]
-
-	fixed_to_one_vars = blmo.int_vars[blmo.lower_bounds .== 1.0]
-	empty!(blmo.fixed_to_one_rows)
-	empty!(blmo.fixed_to_one_cols)
-	for fixed_to_one_var in fixed_to_one_vars
-		q = (fixed_to_one_var - 1) ÷ dim
-		r = (fixed_to_one_var - 1) - q * dim
-		if blmo.append_by_column
-			i = r + 1
-			j = q + 1
-		else
-			i = q + 1
-			j = r + 1
-		end
-		push!(blmo.fixed_to_one_rows, i)
-		push!(blmo.fixed_to_one_cols, j)
-	end
-
-	nfixed = length(blmo.fixed_to_one_rows)
-	nreduced = dim - nfixed
-
-	# stores the indices of the original matrix that are still in the reduced matrix
-	index_map_rows = fill(1, nreduced)
-	index_map_cols = fill(1, nreduced)
-	idx_in_map_row = 1
-	idx_in_map_col = 1
-	for orig_idx in 1:dim
-		if orig_idx ∉ blmo.fixed_to_one_rows
-			index_map_rows[idx_in_map_row] = orig_idx
-			idx_in_map_row += 1
-		end
-		if orig_idx ∉ blmo.fixed_to_one_cols
-			index_map_cols[idx_in_map_col] = orig_idx
-			idx_in_map_col += 1
-		end
-	end
-
-	empty!(blmo.index_map_rows)
-	empty!(blmo.index_map_cols)
-	append!(blmo.index_map_rows, index_map_rows)
-	append!(blmo.index_map_cols, index_map_cols)
+    dim = blmo.dim
+    n = dim^2
+
+    # fixed_to_one_vars = [v for v in blmo.int_vars if blmo.lower_bounds[v] == 1.0]
+
+    fixed_to_one_vars = blmo.int_vars[blmo.lower_bounds.==1.0]
+    empty!(blmo.fixed_to_one_rows)
+    empty!(blmo.fixed_to_one_cols)
+    for fixed_to_one_var in fixed_to_one_vars
+        q = (fixed_to_one_var - 1) ÷ dim
+        r = (fixed_to_one_var - 1) - q * dim
+        if blmo.append_by_column
+            i = r + 1
+            j = q + 1
+        else
+            i = q + 1
+            j = r + 1
+        end
+        push!(blmo.fixed_to_one_rows, i)
+        push!(blmo.fixed_to_one_cols, j)
+    end
+
+    nfixed = length(blmo.fixed_to_one_rows)
+    nreduced = dim - nfixed
+
+    # stores the indices of the original matrix that are still in the reduced matrix
+    index_map_rows = fill(1, nreduced)
+    index_map_cols = fill(1, nreduced)
+    idx_in_map_row = 1
+    idx_in_map_col = 1
+    for orig_idx = 1:dim
+        if orig_idx ∉ blmo.fixed_to_one_rows
+            index_map_rows[idx_in_map_row] = orig_idx
+            idx_in_map_row += 1
+        end
+        if orig_idx ∉ blmo.fixed_to_one_cols
+            index_map_cols[idx_in_map_col] = orig_idx
+            idx_in_map_col += 1
+        end
+    end
+
+    empty!(blmo.index_map_rows)
+    empty!(blmo.index_map_cols)
+    append!(blmo.index_map_rows, index_map_rows)
+    append!(blmo.index_map_cols, index_map_cols)
 end

From 372f0aa3c1ead6e31d0d1ae9d2ba1e7280b306e4 Mon Sep 17 00:00:00 2001
From: Wenjie Xiao <116566662+WenjieXiao-2022@users.noreply.github.com>
Date: Mon, 1 Sep 2025 16:09:28 +0200
Subject: [PATCH 03/21] update dependency

---
 src/CombinatorialLinearOracles.jl | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/src/CombinatorialLinearOracles.jl b/src/CombinatorialLinearOracles.jl
index 343cf4c..9994617 100644
--- a/src/CombinatorialLinearOracles.jl
+++ b/src/CombinatorialLinearOracles.jl
@@ -4,9 +4,12 @@ import FrankWolfe
 using Graphs
 using SparseArrays
 using GraphsMatching
+using Hungarian
+using Boscia
 
 include("matchings.jl")
 include("spanning_tree.jl")
 include("shortest_path.jl")
+include("Birkhoff_polytope.jl")
 
 end

From 595ed15878191c6f047b1b20524d2ddaec9115b2 Mon Sep 17 00:00:00 2001
From: Wenjie Xiao <116566662+WenjieXiao-2022@users.noreply.github.com>
Date: Mon, 8 Sep 2025 12:05:57 +0200
Subject: [PATCH 04/21] make `bounds`  internal arguments

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

diff --git a/src/Birkhoff_polytope.jl b/src/Birkhoff_polytope.jl
index a34a324..89c5d2b 100644
--- a/src/Birkhoff_polytope.jl
+++ b/src/Birkhoff_polytope.jl
@@ -7,8 +7,6 @@ node-specific bounds on the integer variables.
 struct BirkhoffBLMO <: Boscia.BoundedLinearMinimizationOracle
     append_by_column::Bool
     dim::Int
-    lower_bounds::Vector{Float64}
-    upper_bounds::Vector{Float64}
     int_vars::Vector{Int}
     fixed_to_one_rows::Vector{Int}
     fixed_to_one_cols::Vector{Int}
@@ -22,8 +20,8 @@ BirkhoffBLMO(dim, lower_bounds, upper_bounds, int_vars; append_by_column = true)
     BirkhoffBLMO(
         append_by_column,
         dim,
-        lower_bounds,
-        upper_bounds,
+        fill(0.0, length(int_vars)),
+        fill(1.0, length(int_vars)),
         int_vars,
         Int[],
         Int[],

From 85f827eec8ccf21e91b7bbedf7b1a979c9bc5a59 Mon Sep 17 00:00:00 2001
From: Wenjie Xiao <116566662+WenjieXiao-2022@users.noreply.github.com>
Date: Mon, 8 Sep 2025 12:07:55 +0200
Subject: [PATCH 05/21] make `atol` and `rtol` user-defined

---
 src/Birkhoff_polytope.jl | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/Birkhoff_polytope.jl b/src/Birkhoff_polytope.jl
index 89c5d2b..c7a323e 100644
--- a/src/Birkhoff_polytope.jl
+++ b/src/Birkhoff_polytope.jl
@@ -16,7 +16,7 @@ struct BirkhoffBLMO <: Boscia.BoundedLinearMinimizationOracle
     rtol::Float64
 end
 
-BirkhoffBLMO(dim, lower_bounds, upper_bounds, int_vars; append_by_column = true) =
+BirkhoffBLMO(dim, lower_bounds, upper_bounds, int_vars; append_by_column = true, atol=1e-6, rtol=1e-3) =
     BirkhoffBLMO(
         append_by_column,
         dim,
@@ -27,8 +27,8 @@ BirkhoffBLMO(dim, lower_bounds, upper_bounds, int_vars; append_by_column = true)
         Int[],
         collect(1:dim),
         collect(1:dim),
-        1e-6,
-        1e-3,
+        atol,
+        rtol,
     )
 
 ## Necessary

From 9b12464ea651ba1848325e4b7e88d586706f83b3 Mon Sep 17 00:00:00 2001
From: Wenjie Xiao <116566662+WenjieXiao-2022@users.noreply.github.com>
Date: Mon, 8 Sep 2025 12:11:00 +0200
Subject: [PATCH 06/21] simplify `dicg_maximum_step` function

---
 src/Birkhoff_polytope.jl | 29 ++++++++++++++++++++++++-----
 1 file changed, 24 insertions(+), 5 deletions(-)

diff --git a/src/Birkhoff_polytope.jl b/src/Birkhoff_polytope.jl
index c7a323e..20b755d 100644
--- a/src/Birkhoff_polytope.jl
+++ b/src/Birkhoff_polytope.jl
@@ -215,12 +215,31 @@ Fixings are maintained by the oracle (or deduced from `x` itself).
 """
 function Boscia.dicg_maximum_step(blmo::BirkhoffBLMO, direction, x; kwargs...)
     n = blmo.dim
+	T = promote_type(eltype(x), eltype(direction))
+    if size(direction, 2) == 1
+        direction =
+            blmo.append_by_column ? reshape(direction, (n, n)) :
+            transpose(reshape(direction, (n, n)))
+        x = blmo.append_by_column ? reshape(x, (n, n)) : transpose(reshape(x, (n, n)))
+    end
 
-    direction =
-        blmo.append_by_column ? reshape(direction, (n, n)) :
-        transpose(reshape(direction, (n, n)))
-    x = blmo.append_by_column ? reshape(x, (n, n)) : transpose(reshape(x, (n, n)))
-    return FrankWolfe.dicg_maximum_step(FrankWolfe.BirkhoffPolytopeLMO(), direction, x)
+    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(blmo::BirkhoffBLMO)

From eeb68d147f3e2f93a33fcf82bc5b1f4f7799f03c Mon Sep 17 00:00:00 2001
From: Wenjie Xiao <116566662+WenjieXiao-2022@users.noreply.github.com>
Date: Mon, 8 Sep 2025 13:48:04 +0200
Subject: [PATCH 07/21] refactor add/delete functions

---
 src/Birkhoff_polytope.jl | 59 ++++++++++++++++++++++++++++++++++++++++
 1 file changed, 59 insertions(+)

diff --git a/src/Birkhoff_polytope.jl b/src/Birkhoff_polytope.jl
index 20b755d..518368b 100644
--- a/src/Birkhoff_polytope.jl
+++ b/src/Birkhoff_polytope.jl
@@ -320,6 +320,20 @@ end
 function Boscia.set_bound!(blmo::BirkhoffBLMO, c_idx, value, sense::Symbol)
     if sense == :greaterthan
         blmo.lower_bounds[c_idx] = value
+		if value == 1.0
+            n0 = blmo.dim
+            fixed_int_var = blmo.int_vars[c_idx]
+            # Convert linear index to (row, col) based on storage format
+            if blmo.append_by_column
+                j = ceil(Int, fixed_int_var / n0)  # column index
+                i = Int(fixed_int_var - n0 * (j - 1))  # row index
+            else
+                i = ceil(Int, fixed_int_var / n0)  # row index  
+                j = Int(fixed_int_var - n0 * (i - 1))  # column index
+            end
+            push!(blmo.fixed_to_one_rows, i)
+            push!(blmo.fixed_to_one_cols, j)
+        end
     elseif sense == :lessthan
         blmo.upper_bounds[c_idx] = value
     else
@@ -336,6 +350,51 @@ function Boscia.delete_bounds!(blmo::BirkhoffBLMO, cons_delete)
             blmo.upper_bounds[d_idx] = Inf
         end
     end
+
+	# sanity check
+    check_feasibility(blmo, i::Int, j::Int) =
+        blmo.lower_bounds[blmo.append_by_column ? (j-1)*blmo.dim + i : (i-1)*blmo.dim + j] != 0.0
+
+    fixed_to_one_rows = blmo.fixed_to_one_rows
+    fixed_to_one_cols = blmo.fixed_to_one_cols
+
+    feasible_flags  = check_feasibility.(Ref(blmo), fixed_to_one_rows, fixed_to_one_cols)
+    invalid_indices = findall(.!feasible_flags)
+
+    deleteat!(fixed_to_one_rows, invalid_indices)
+    deleteat!(fixed_to_one_cols, invalid_indices)
+
+    # remove the duplicate indices
+    pairs = collect(zip(blmo.fixed_to_one_rows, blmo.fixed_to_one_cols))
+    unique!(pairs)
+    resize!(blmo.fixed_to_one_rows, length(pairs))
+    resize!(blmo.fixed_to_one_cols, length(pairs))
+
+    blmo.fixed_to_one_rows .= first.(pairs)
+    blmo.fixed_to_one_cols .= last.(pairs)
+
+    nfixed = length(blmo.fixed_to_one_rows)
+    nreduced = blmo.dim - nfixed
+
+    # stores the indices of the original matrix that are still in the reduced matrix
+    index_map_rows = fill(1, nreduced)
+    index_map_cols = fill(1, nreduced)
+    idx_in_map_row = 1
+    idx_in_map_col = 1
+    for orig_idx in 1:blmo.dim
+        if orig_idx ∉ blmo.fixed_to_one_rows
+            index_map_rows[idx_in_map_row] = orig_idx
+            idx_in_map_row += 1
+        end
+        if orig_idx ∉ blmo.fixed_to_one_cols
+            index_map_cols[idx_in_map_col] = orig_idx
+            idx_in_map_col += 1
+        end
+    end
+
+    empty!(blmo.index_map_rows)
+    empty!(blmo.index_map_cols)
+    append!(blmo.index_map_rows, index_map_rows)
 end
 
 # Add bound constraint.

From e5067917a0e3809255aa59d47a011513521b3c4a Mon Sep 17 00:00:00 2001
From: Wenjie Xiao <116566662+WenjieXiao-2022@users.noreply.github.com>
Date: Mon, 8 Sep 2025 13:48:29 +0200
Subject: [PATCH 08/21] remove `update_model` function

---
 src/Birkhoff_polytope.jl | 48 ----------------------------------------
 1 file changed, 48 deletions(-)

diff --git a/src/Birkhoff_polytope.jl b/src/Birkhoff_polytope.jl
index 518368b..69a7302 100644
--- a/src/Birkhoff_polytope.jl
+++ b/src/Birkhoff_polytope.jl
@@ -479,51 +479,3 @@ function Boscia.check_feasibility(blmo::BirkhoffBLMO)
 
     return Boscia.OPTIMAL
 end
-
-function Boscia.update_model(blmo::BirkhoffBLMO)
-    dim = blmo.dim
-    n = dim^2
-
-    # fixed_to_one_vars = [v for v in blmo.int_vars if blmo.lower_bounds[v] == 1.0]
-
-    fixed_to_one_vars = blmo.int_vars[blmo.lower_bounds.==1.0]
-    empty!(blmo.fixed_to_one_rows)
-    empty!(blmo.fixed_to_one_cols)
-    for fixed_to_one_var in fixed_to_one_vars
-        q = (fixed_to_one_var - 1) ÷ dim
-        r = (fixed_to_one_var - 1) - q * dim
-        if blmo.append_by_column
-            i = r + 1
-            j = q + 1
-        else
-            i = q + 1
-            j = r + 1
-        end
-        push!(blmo.fixed_to_one_rows, i)
-        push!(blmo.fixed_to_one_cols, j)
-    end
-
-    nfixed = length(blmo.fixed_to_one_rows)
-    nreduced = dim - nfixed
-
-    # stores the indices of the original matrix that are still in the reduced matrix
-    index_map_rows = fill(1, nreduced)
-    index_map_cols = fill(1, nreduced)
-    idx_in_map_row = 1
-    idx_in_map_col = 1
-    for orig_idx = 1:dim
-        if orig_idx ∉ blmo.fixed_to_one_rows
-            index_map_rows[idx_in_map_row] = orig_idx
-            idx_in_map_row += 1
-        end
-        if orig_idx ∉ blmo.fixed_to_one_cols
-            index_map_cols[idx_in_map_col] = orig_idx
-            idx_in_map_col += 1
-        end
-    end
-
-    empty!(blmo.index_map_rows)
-    empty!(blmo.index_map_cols)
-    append!(blmo.index_map_rows, index_map_rows)
-    append!(blmo.index_map_cols, index_map_cols)
-end

From 44d0e54dae1dac3b72075d2949b69fd045d0e02a Mon Sep 17 00:00:00 2001
From: Wenjie Xiao <116566662+WenjieXiao-2022@users.noreply.github.com>
Date: Mon, 8 Sep 2025 13:48:50 +0200
Subject: [PATCH 09/21] Rename Birkhoff_polytope.jl to birkhoff_polytope.jl

---
 src/{Birkhoff_polytope.jl => birkhoff_polytope.jl} | 0
 1 file changed, 0 insertions(+), 0 deletions(-)
 rename src/{Birkhoff_polytope.jl => birkhoff_polytope.jl} (100%)

diff --git a/src/Birkhoff_polytope.jl b/src/birkhoff_polytope.jl
similarity index 100%
rename from src/Birkhoff_polytope.jl
rename to src/birkhoff_polytope.jl

From 9e5fe4444dfec2b7a1bd8c09f8eaf20425f30c9f Mon Sep 17 00:00:00 2001
From: Wenjie Xiao <116566662+WenjieXiao-2022@users.noreply.github.com>
Date: Mon, 8 Sep 2025 13:49:17 +0200
Subject: [PATCH 10/21] rename file

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

diff --git a/src/CombinatorialLinearOracles.jl b/src/CombinatorialLinearOracles.jl
index 9994617..2b8ab0c 100644
--- a/src/CombinatorialLinearOracles.jl
+++ b/src/CombinatorialLinearOracles.jl
@@ -10,6 +10,6 @@ using Boscia
 include("matchings.jl")
 include("spanning_tree.jl")
 include("shortest_path.jl")
-include("Birkhoff_polytope.jl")
+include("birkhoff_polytope.jl")
 
 end

From 0a9acde4797f60478e4902094d80a204f3e0337c Mon Sep 17 00:00:00 2001
From: Wenjie Xiao <116566662+WenjieXiao-2022@users.noreply.github.com>
Date: Tue, 9 Sep 2025 10:22:12 +0200
Subject: [PATCH 11/21] format

---
 src/birkhoff_polytope.jl | 1 -
 1 file changed, 1 deletion(-)

diff --git a/src/birkhoff_polytope.jl b/src/birkhoff_polytope.jl
index 69a7302..d3df39a 100644
--- a/src/birkhoff_polytope.jl
+++ b/src/birkhoff_polytope.jl
@@ -36,7 +36,6 @@ BirkhoffBLMO(dim, lower_bounds, upper_bounds, int_vars; append_by_column = true,
 """
 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.compute_extreme_point(blmo::BirkhoffBLMO, d; kwargs...)
     n = blmo.dim
 

From b89eaa03b703dab0720371ff197ecb8f7b1435f9 Mon Sep 17 00:00:00 2001
From: Wenjie Xiao <116566662+WenjieXiao-2022@users.noreply.github.com>
Date: Tue, 9 Sep 2025 12:27:35 +0200
Subject: [PATCH 12/21] interfaces for both FW and Boscia

---
 src/birkhoff_polytope_interfaces.jl | 639 ++++++++++++++++++++++++++++
 1 file changed, 639 insertions(+)
 create mode 100644 src/birkhoff_polytope_interfaces.jl

diff --git a/src/birkhoff_polytope_interfaces.jl b/src/birkhoff_polytope_interfaces.jl
new file mode 100644
index 0000000..14c01d2
--- /dev/null
+++ b/src/birkhoff_polytope_interfaces.jl
@@ -0,0 +1,639 @@
+"""
+	BirkhoffBLMO
+
+A bounded LMO for the Birkhoff polytope. This oracle computes an extreme point subject to  
+node-specific bounds on the integer variables.
+"""
+struct BirkhoffLMO
+    append_by_column::Bool
+    dim::Int
+    lower_bounds::Vector{Float64}
+    upper_bounds::Vector{Float64}
+    int_vars::Vector{Int}
+    fixed_to_one_rows::Vector{Int}
+    fixed_to_one_cols::Vector{Int}
+    index_map_rows::Vector{Int}
+    index_map_cols::Vector{Int}
+    atol::Float64
+    rtol::Float64
+end
+
+struct BosciaBirkhoffBLMO <: Boscia.BoundedLinearMinimizationOracle
+    birkhoffBLMO::BirkhoffLMO
+end
+
+struct FWBirkhoffLMO <: FrankWolfe.LinearMinimizationOracle
+    birkhoffLMO::BirkhoffLMO
+end
+
+# define the mixed-integer Birkhoff polytope
+function BirkhoffBLMO(dim, int_vars; append_by_column=true, atol=1e-6, rtol=1e-3)
+    birkhoffLMO = BirkhoffLMO(
+        append_by_column,
+        dim,
+        fill(0.0, length(int_vars)),
+        fill(1.0, length(int_vars)),
+        int_vars,
+        Int[],
+        Int[],
+        collect(1:dim),
+        collect(1:dim),
+        atol,
+        rtol,
+    )
+    return BosciaBirkhoffBLMO(birkhoffLMO)
+end
+
+function BirkhoffBLMO(dim; append_by_column=true, use_Boscia=true, atol=1e-6, rtol=1e-3)
+    # define the all-integer Birkhoff polytope
+    if use_Boscia
+        birkhoffLMO = BirkhoffLMO(
+            append_by_column,
+            dim,
+            fill(0.0, length(int_vars)),
+            fill(1.0, length(int_vars)),
+            collect(1:(dim^2)),
+            Int[],
+            Int[],
+            collect(1:dim),
+            collect(1:dim),
+            atol,
+            rtol,
+        )
+        return BosciaBirkhoffBLMO(birkhoffLMO)
+    else
+        # define the continous Birkhoff polytope
+        birkhoffLMO = BirkhoffLMO(
+            append_by_column,
+            dim,
+            [],
+            [],
+            [],
+            Int[],
+            Int[],
+            collect(1:dim),
+            collect(1:dim),
+            atol,
+            rtol,
+        )
+        return FWBirkhoffLMO(birkhoffLMO)
+    end
+
+end
+
+"""
+Computes the extreme point given an direction d, the current lower and upper bounds on the integer variables, and the set of integer variables.
+"""
+_compute_extreme_point_impl(blmo::BirkhoffLMO, d; kwargs...) = begin
+    n = blmo.dim
+
+    # Precompute index mapping to avoid repeated `findfirst` calls,
+    # which would be very costly inside the loop.
+    if length(blmo.int_vars) !== n^2
+        idx_map_ub = zeros(Int, n^2)
+        @inbounds for (c_idx, var) in enumerate(blmo.int_vars)
+            idx_map_ub[var] = c_idx
+        end
+    end
+
+    if size(d, 2) == 1
+        d = blmo.append_by_column ? reshape(d, (n, n)) : transpose(reshape(d, (n, n)))
+    end
+
+    fixed_to_one_rows = blmo.fixed_to_one_rows
+    fixed_to_one_cols = blmo.fixed_to_one_cols
+    index_map_rows = blmo.index_map_rows
+    index_map_cols = blmo.index_map_cols
+    int_vars = blmo.int_vars
+    ub = blmo.upper_bounds
+
+    nreduced = length(index_map_rows)
+    type = typeof(d[1, 1])
+    d2 = ones(Union{type,Missing}, nreduced, nreduced)
+    for j in 1:nreduced
+        col_orig = index_map_cols[j]
+        for i in 1:nreduced
+            row_orig = index_map_rows[i]
+            if blmo.append_by_column
+                orig_linear_idx = (col_orig - 1) * n + row_orig
+            else
+                orig_linear_idx = (row_orig - 1) * n + col_orig
+            end
+            # the problem can only be integer types,
+            # either all-integer or mixed-integer.
+            if orig_linear_idx in int_vars
+                idx = length(int_vars) < n^2 ? idx_map_ub[orig_linear_idx] : orig_linear_idx
+                # interdict arc when fixed to zero
+                if ub[idx] <= eps()
+                    if blmo.append_by_column
+                        d2[i, j] = missing
+                    else
+                        d2[j, i] = missing
+                    end
+                else
+                    if blmo.append_by_column
+                        d2[i, j] = d[row_orig, col_orig]
+                    else
+                        d2[j, i] = d[col_orig, row_orig]
+                    end
+                end
+            else
+                if blmo.append_by_column
+                    d2[i, j] = d[row_orig, col_orig]
+                else
+                    d2[j, i] = d[col_orig, row_orig]
+                end
+            end
+        end
+    end
+
+    m = SparseArrays.spzeros(n, n)
+    for (i, j) in zip(fixed_to_one_rows, fixed_to_one_cols)
+        m[i, j] = 1
+    end
+
+    res_mat = Hungarian.munkres(d2)
+    (rows, cols, vals) = SparseArrays.findnz(res_mat)
+    @inbounds for i in eachindex(cols)
+        m[index_map_rows[rows[i]], index_map_cols[cols[i]]] = (vals[i] == 2)
+    end
+
+    # if size(d, 2) == 1
+    m = if blmo.append_by_column
+        # Convert sparse matrix to sparse vector by columns
+        I, J, V = SparseArrays.findnz(m)
+        linear_indices = (J .- 1) .* n .+ I
+        SparseArrays.sparsevec(linear_indices, V, n^2)
+    else
+        # Convert sparse matrix to sparse vector by rows (transpose first)
+        mt = SparseArrays.sparse(LinearAlgebra.transpose(m))
+        I, J, V = SparseArrays.findnz(mt)
+        linear_indices = (J .- 1) .* n .+ I
+        SparseArrays.sparsevec(linear_indices, V, n^2)
+    end
+    # end
+    return m
+end
+
+"""
+Computes the extreme point given an direction d, the current lower and upper bounds on the integer variables, and the set of integer variables.
+"""
+_compute_inface_extreme_point_impl(blmo::BirkhoffLMO, direction, x; kwargs...) = begin
+    n = blmo.dim
+
+    if size(direction, 2) == 1
+        direction =
+            blmo.append_by_column ? reshape(direction, (n, n)) :
+            transpose(reshape(direction, (n, n)))
+
+        x = blmo.append_by_column ? reshape(x, (n, n)) : transpose(reshape(x, (n, n)))
+    end
+
+    # Precompute index mapping to avoid repeated `findfirst` calls,
+    # which would be very costly inside the loop.
+    if length(blmo.int_vars) !== n^2
+        idx_map_ub = zeros(Int, n^2)
+        @inbounds for (c_idx, var) in enumerate(blmo.int_vars)
+            idx_map_ub[var] = c_idx
+        end
+    end
+
+    fixed_to_one_rows = copy(blmo.fixed_to_one_rows)
+    fixed_to_one_cols = copy(blmo.fixed_to_one_cols)
+    index_map_rows = copy(blmo.index_map_rows)
+    index_map_cols = copy(blmo.index_map_cols)
+    int_vars = blmo.int_vars
+    ub = blmo.upper_bounds
+
+    nreduced = length(blmo.index_map_rows)
+
+    delete_index_map_rows = Int[]
+    delete_index_map_cols = Int[]
+    for j in 1:nreduced
+        for i in 1:nreduced
+            row_orig = index_map_rows[i]
+            col_orig = index_map_cols[j]
+            if x[row_orig, col_orig] >= 1-eps()
+                push!(fixed_to_one_rows, row_orig)
+                push!(fixed_to_one_cols, col_orig)
+
+                push!(delete_index_map_rows, i)
+                push!(delete_index_map_cols, j)
+            end
+        end
+    end
+
+    unique!(delete_index_map_rows)
+    unique!(delete_index_map_cols)
+    sort!(delete_index_map_rows)
+    sort!(delete_index_map_cols)
+    deleteat!(index_map_rows, delete_index_map_rows)
+    deleteat!(index_map_cols, delete_index_map_cols)
+
+    nreduced = length(index_map_rows)
+    type = typeof(direction[1, 1])
+    d2 = ones(Union{type,Missing}, nreduced, nreduced)
+    for j in 1:nreduced
+        for i in 1:nreduced
+            row_orig = index_map_rows[i]
+            col_orig = index_map_cols[j]
+            if blmo.append_by_column
+                orig_linear_idx = (col_orig-1)*n+row_orig
+            else
+                orig_linear_idx = (row_orig-1)*n+col_orig
+            end
+            # idx = findfirst(x -> x == orig_linear_idx, blmo.int_vars)
+
+            if x[row_orig, col_orig] <= eps()
+                if blmo.append_by_column
+                    d2[i, j] = missing
+                else
+                    d2[j, i] = missing
+                end
+            elseif orig_linear_idx in int_vars
+                idx = length(int_vars) < n^2 ? idx_map_ub[orig_linear_idx] : orig_linear_idx
+                # interdict arc when fixed to zero
+                if ub[idx] <= eps()
+                    if blmo.append_by_column
+                        d2[i, j] = missing
+                    else
+                        d2[j, i] = missing
+                    end
+                else
+                    if blmo.append_by_column
+                        d2[i, j] = direction[row_orig, col_orig]
+                    else
+                        d2[j, i] = direction[col_orig, row_orig]
+                    end
+                end
+            else
+                if blmo.append_by_column
+                    d2[i, j] = direction[row_orig, col_orig]
+                else
+                    d2[j, i] = direction[col_orig, row_orig]
+                end
+            end
+        end
+    end
+
+    m = SparseArrays.spzeros(n, n)
+    for (i, j) in zip(fixed_to_one_rows, fixed_to_one_cols)
+        m[i, j] = 1
+    end
+
+    res_mat = Hungarian.munkres(d2)
+    (rows, cols, vals) = SparseArrays.findnz(res_mat)
+    @inbounds for i in eachindex(cols)
+        m[index_map_rows[rows[i]], index_map_cols[cols[i]]] = (vals[i] == 2)
+    end
+
+    # if size(direction, 2) == 1
+    m = if blmo.append_by_column
+        # Convert sparse matrix to sparse vector by columns
+        I, J, V = SparseArrays.findnz(m)
+        linear_indices = (J .- 1) .* n .+ I
+        SparseArrays.sparsevec(linear_indices, V, n^2)
+    else
+        # Convert sparse matrix to sparse vector by rows (transpose first)
+        mt = SparseArrays.sparse(LinearAlgebra.transpose(m))
+        I, J, V = SparseArrays.findnz(mt)
+        linear_indices = (J .- 1) .* n .+ I
+        SparseArrays.sparsevec(linear_indices, V, n^2)
+    end
+    # end
+
+    return m
+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).
+"""
+_dicg_maximum_step_impl(blmo::BirkhoffLMO, direction, x; kwargs...) = begin
+    n = blmo.dim
+
+    T = promote_type(eltype(x), eltype(direction))
+    if size(direction, 2) == 1
+        direction =
+            blmo.append_by_column ? reshape(direction, (n, n)) :
+            transpose(reshape(direction, (n, n)))
+
+        x = blmo.append_by_column ? reshape(x, (n, n)) : transpose(reshape(x, (n, n)))
+    end
+
+    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
+
+_is_decomposition_invariant_oracle_impl(lmo::BirkhoffLMO) = true
+
+
+# Boscia API: method extensions for BirkhoffBLMO
+Boscia.compute_extreme_point(blmo::BosciaBirkhoffBLMO, d; kwargs...) =
+    _compute_extreme_point_impl(blmo.birkhoffBLMO, d; kwargs...)
+
+Boscia.compute_inface_extreme_point(blmo::BosciaBirkhoffBLMO, direction, x; kwargs...) =
+    _compute_inface_extreme_point_impl(blmo.birkhoffBLMO, direction, x; kwargs...)
+
+
+Boscia.dicg_maximum_step(blmo::BosciaBirkhoffBLMO, direction, x; kwargs...) =
+    _dicg_maximum_step_impl(blmo.birkhoffBLMO, direction, x; kwargs...)
+
+Boscia.is_decomposition_invariant_oracle(blmo::BosciaBirkhoffBLMO) =
+    _is_decomposition_invariant_oracle_impl(blmo.birkhoffBLMO)
+
+
+# FrankWolfe API: method extensions for BirkhoffBLMO
+FrankWolfe.compute_extreme_point(lmo::FWBirkhoffLMO, d; kwargs...) =
+    _compute_extreme_point_impl(lmo.birkhoffBLMO, d; kwargs...)
+
+FrankWolfe.compute_inface_extreme_point(lmo::FWBirkhoffLMO, direction, x; kwargs...) =
+    _compute_inface_extreme_point_impl(lmo.birkhoffBLMO, direction, x; kwargs...)
+
+FrankWolfe.dicg_maximum_step(lmo::FWBirkhoffLMO, direction, x; kwargs...) =
+    _dicg_maximum_step_impl(lmo.birkhoffBLMO, direction, x; kwargs...)
+
+FrankWolfe.is_decomposition_invariant_oracle(lmo::FWBirkhoffLMO) =
+    _is_decomposition_invariant_oracle_impl(lmo.birkhoffLMO)
+
+
+## Necessary for Boscia package
+
+"""
+The sum of each row and column has to be equal to 1.
+"""
+function Boscia.is_linear_feasible(lmo::BosciaBirkhoffBLMO, v::AbstractVector)
+    blmo = lmo.birkhoffBLMO
+    n = blmo.dim
+    for i in 1:n
+        # append by column ? column sum : row sum 
+        if !isapprox(sum(v[((i-1)*n+1):(i*n)]), 1.0, atol=1e-6, rtol=1e-3)
+            @debug "Column sum not 1: $(sum(v[((i-1)*n+1):(i*n)]))"
+            return false
+        end
+        # append by column ? row sum : column sum
+        if !isapprox(sum(v[i: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
+
+# Read global bounds from the problem.
+function Boscia.build_global_bounds(lmo::BosciaBirkhoffBLMO, integer_variables)
+    blmo = lmo.birkhoffBLMO
+    global_bounds = Boscia.IntegerBounds()
+    for (idx, int_var) in enumerate(blmo.int_vars)
+        push!(global_bounds, (int_var, blmo.lower_bounds[idx]), :greaterthan)
+        push!(global_bounds, (int_var, blmo.upper_bounds[idx]), :lessthan)
+    end
+    return global_bounds
+end
+
+# Get list of variables indices. 
+# If the problem has n variables, they are expected to contiguous and ordered from 1 to n.
+function Boscia.get_list_of_variables(lmo::BosciaBirkhoffBLMO)
+    blmo = lmo.birkhoffBLMO
+    n = blmo.dim^2
+    return n, collect(1:n)
+end
+
+# Get list of integer variables
+function Boscia.get_integer_variables(lmo::BosciaBirkhoffBLMO)
+    blmo = lmo.birkhoffBLMO
+    return blmo.int_vars
+end
+
+# Get the index of the integer variable the bound is working on.
+function Boscia.get_int_var(lmo::BosciaBirkhoffBLMO, cidx)
+    blmo = lmo.birkhoffBLMO
+    return blmo.int_vars[cidx]
+end
+
+# Get the list of lower bounds.
+function Boscia.get_lower_bound_list(lmo::BosciaBirkhoffBLMO)
+    blmo = lmo.birkhoffBLMO
+    return collect(1:length(blmo.lower_bounds))
+end
+
+# Get the list of upper bounds.
+function Boscia.get_upper_bound_list(lmo::BosciaBirkhoffBLMO)
+    blmo = lmo.birkhoffBLMO
+    return collect(1:length(blmo.upper_bounds))
+end
+
+# Read bound value for c_idx.
+function Boscia.get_bound(lmo::BosciaBirkhoffBLMO, c_idx, sense::Symbol)
+    blmo = lmo.birkhoffBLMO
+    if sense == :lessthan
+        return blmo.upper_bounds[c_idx]
+    elseif sense == :greaterthan
+        return blmo.lower_bounds[c_idx]
+    else
+        error("Allowed value for sense are :lessthan and :greaterthan!")
+    end
+end
+
+## Changing the bounds constraints.
+
+# Change the value of the bound c_idx.
+function Boscia.set_bound!(lmo::BosciaBirkhoffBLMO, c_idx, value, sense::Symbol)
+    blmo = lmo.birkhoffBLMO
+    if sense == :greaterthan
+        blmo.lower_bounds[c_idx] = value
+        if value == 1.0
+            n0 = blmo.dim
+            fixed_int_var = blmo.int_vars[c_idx]
+            # Convert linear index to (row, col) based on storage format
+            if blmo.append_by_column
+                j = ceil(Int, fixed_int_var / n0)  # column index
+                i = Int(fixed_int_var - n0 * (j - 1))  # row index
+            else
+                i = ceil(Int, fixed_int_var / n0)  # row index  
+                j = Int(fixed_int_var - n0 * (i - 1))  # column index
+            end
+            push!(blmo.fixed_to_one_rows, i)
+            push!(blmo.fixed_to_one_cols, j)
+        end
+    elseif sense == :lessthan
+        blmo.upper_bounds[c_idx] = value
+    else
+        error("Allowed values for sense are :lessthan and :greaterthan.")
+    end
+end
+
+# Delete bounds.
+function Boscia.delete_bounds!(lmo::BosciaBirkhoffBLMO, cons_delete)
+    blmo = lmo.birkhoffBLMO
+    for (d_idx, sense) in cons_delete
+        if sense == :greaterthan
+            blmo.lower_bounds[d_idx] = -Inf
+        else
+            blmo.upper_bounds[d_idx] = Inf
+        end
+    end
+
+    # sanity check 
+    check_feasibility(blmo, i::Int, j::Int) =
+        blmo.lower_bounds[blmo.append_by_column ? (j-1)*blmo.dim + i : (i-1)*blmo.dim + j] != 0.0
+
+    fixed_to_one_rows = blmo.fixed_to_one_rows
+    fixed_to_one_cols = blmo.fixed_to_one_cols
+
+    feasible_flags = check_feasibility.(Ref(blmo), fixed_to_one_rows, fixed_to_one_cols)
+    invalid_indices = findall(.!feasible_flags)
+
+    deleteat!(fixed_to_one_rows, invalid_indices)
+    deleteat!(fixed_to_one_cols, invalid_indices)
+
+    # remove the duplicate indices
+    pairs = collect(zip(blmo.fixed_to_one_rows, blmo.fixed_to_one_cols))
+    unique!(pairs)
+    resize!(blmo.fixed_to_one_rows, length(pairs))
+    resize!(blmo.fixed_to_one_cols, length(pairs))
+
+    blmo.fixed_to_one_rows .= first.(pairs)
+    blmo.fixed_to_one_cols .= last.(pairs)
+
+    nfixed = length(blmo.fixed_to_one_rows)
+    nreduced = blmo.dim - nfixed
+
+    # stores the indices of the original matrix that are still in the reduced matrix
+    index_map_rows = fill(1, nreduced)
+    index_map_cols = fill(1, nreduced)
+    idx_in_map_row = 1
+    idx_in_map_col = 1
+    for orig_idx in 1:blmo.dim
+        if orig_idx ∉ blmo.fixed_to_one_rows
+            index_map_rows[idx_in_map_row] = orig_idx
+            idx_in_map_row += 1
+        end
+        if orig_idx ∉ blmo.fixed_to_one_cols
+            index_map_cols[idx_in_map_col] = orig_idx
+            idx_in_map_col += 1
+        end
+    end
+
+    empty!(blmo.index_map_rows)
+    empty!(blmo.index_map_cols)
+    append!(blmo.index_map_rows, index_map_rows)
+    append!(blmo.index_map_cols, index_map_cols)
+
+    return blmo.modify_LMO = false
+end
+
+# Add bound constraint.
+function Boscia.add_bound_constraint!(lmo::BosciaBirkhoffBLMO, key, value, sense::Symbol)
+    blmo = lmo.birkhoffBLMO
+    idx = findfirst(x -> x == key, blmo.int_vars)
+    if sense == :greaterthan
+        blmo.lower_bounds[idx] = value
+
+    elseif sense == :lessthan
+        blmo.upper_bounds[idx] = value
+    else
+        error("Allowed value of sense are :lessthan and :greaterthan!")
+    end
+end
+
+## Checks
+
+# Check if the subject of the bound c_idx is an integer variable (recorded in int_vars).
+function Boscia.is_constraint_on_int_var(lmo::BosciaBirkhoffBLMO, c_idx, int_vars)
+    blmo = lmo.birkhoffBLMO
+    return blmo.int_vars[c_idx] in int_vars
+end
+
+# To check if there is bound for the variable in the global or node bounds.
+function Boscia.is_bound_in(lmo::BosciaBirkhoffBLMO, c_idx, bounds)
+    blmo = lmo.birkhoffBLMO
+    return haskey(bounds, blmo.int_vars[c_idx])
+end
+
+# Has variable an integer constraint?
+function Boscia.has_integer_constraint(lmo::BosciaBirkhoffBLMO, idx)
+    blmo = lmo.birkhoffBLMO
+    return idx in blmo.int_vars
+end
+
+## Optional
+
+function Boscia.check_feasibility(lmo::BosciaBirkhoffBLMO)
+    blmo = lmo.birkhoffBLMO
+    for (lb, ub) in zip(blmo.lower_bounds, blmo.upper_bounds)
+        if ub < lb
+            return Boscia.INFEASIBLE
+        end
+    end
+    # For double stochastic matrices, each row and column must sum to 1
+    # We check if the bounds allow for feasible assignments
+    n0 = blmo.dim
+    n = n0^2
+    int_vars = blmo.int_vars
+    # Initialize row and column bound tracking
+    row_min_sum = zeros(n0)  # minimum possible sum for each row
+    row_max_sum = zeros(n0)  # maximum possible sum for each row
+    col_min_sum = zeros(n0)  # minimum possible sum for each column
+    col_max_sum = zeros(n0)  # maximum possible sum for each column
+
+    rows_with_integer_variables = Int[]
+    cols_with_integer_variables = Int[]
+
+    # Process each integer variable
+    for idx in eachindex(int_vars)
+        var_idx = int_vars[idx]
+
+        # Convert linear index to (row, col) based on storage format
+        if blmo.append_by_column
+            j = ceil(Int, var_idx / n0)  # column index
+            i = Int(var_idx - n0 * (j - 1))  # row index
+        else
+            i = ceil(Int, var_idx / n0)  # row index  
+            j = Int(var_idx - n0 * (i - 1))  # column index
+        end
+
+        # Add bounds to row and column sums
+        row_min_sum[i] += blmo.lower_bounds[idx]
+        row_max_sum[i] += blmo.upper_bounds[idx]
+        col_min_sum[j] += blmo.lower_bounds[idx]
+        col_max_sum[j] += blmo.upper_bounds[idx]
+
+        push!(rows_with_integer_variables, i)
+        push!(cols_with_integer_variables, j)
+    end
+
+    # Check feasibility: each row and column must be able to sum to exactly 1
+    for i in 1:n0
+        if i in rows_with_integer_variables
+            # Check row sum constraints
+            if row_min_sum[i] > 1 + eps() || row_max_sum[i] < 1 - eps()
+                return Boscia.INFEASIBLE
+            end
+        end
+
+        if i in cols_with_integer_variables
+            # Check column sum constraints  
+            if col_min_sum[i] > 1 + eps() || col_max_sum[i] < 1 - eps()
+                return Boscia.INFEASIBLE
+            end
+        end
+    end
+
+    return Boscia.OPTIMAL
+end

From c3cd313b34ea00622c9167a26ac9417629c6e843 Mon Sep 17 00:00:00 2001
From: Wenjie Xiao <116566662+WenjieXiao-2022@users.noreply.github.com>
Date: Tue, 9 Sep 2025 12:40:23 +0200
Subject: [PATCH 13/21] minor fix

---
 src/birkhoff_polytope_interfaces.jl | 9 +++++----
 1 file changed, 5 insertions(+), 4 deletions(-)

diff --git a/src/birkhoff_polytope_interfaces.jl b/src/birkhoff_polytope_interfaces.jl
index 14c01d2..f76e5cf 100644
--- a/src/birkhoff_polytope_interfaces.jl
+++ b/src/birkhoff_polytope_interfaces.jl
@@ -358,15 +358,15 @@ Boscia.is_decomposition_invariant_oracle(blmo::BosciaBirkhoffBLMO) =
     _is_decomposition_invariant_oracle_impl(blmo.birkhoffBLMO)
 
 
-# FrankWolfe API: method extensions for BirkhoffBLMO
+# FrankWolfe API: method extensions for BirkhoffLMO
 FrankWolfe.compute_extreme_point(lmo::FWBirkhoffLMO, d; kwargs...) =
-    _compute_extreme_point_impl(lmo.birkhoffBLMO, d; kwargs...)
+    _compute_extreme_point_impl(lmo.birkhoffLMO, d; kwargs...)
 
 FrankWolfe.compute_inface_extreme_point(lmo::FWBirkhoffLMO, direction, x; kwargs...) =
-    _compute_inface_extreme_point_impl(lmo.birkhoffBLMO, direction, x; kwargs...)
+    _compute_inface_extreme_point_impl(lmo.birkhoffLMO, direction, x; kwargs...)
 
 FrankWolfe.dicg_maximum_step(lmo::FWBirkhoffLMO, direction, x; kwargs...) =
-    _dicg_maximum_step_impl(lmo.birkhoffBLMO, direction, x; kwargs...)
+    _dicg_maximum_step_impl(lmo.birkhoffLMO, direction, x; kwargs...)
 
 FrankWolfe.is_decomposition_invariant_oracle(lmo::FWBirkhoffLMO) =
     _is_decomposition_invariant_oracle_impl(lmo.birkhoffLMO)
@@ -637,3 +637,4 @@ function Boscia.check_feasibility(lmo::BosciaBirkhoffBLMO)
 
     return Boscia.OPTIMAL
 end
+

From f640bd6295b279d34645c2d6f096476c3c262262 Mon Sep 17 00:00:00 2001
From: Wenjie Xiao <116566662+WenjieXiao-2022@users.noreply.github.com>
Date: Wed, 10 Sep 2025 21:18:34 +0200
Subject: [PATCH 14/21] Delete src/birkhoff_polytope_interfaces.jl

---
 src/birkhoff_polytope_interfaces.jl | 640 ----------------------------
 1 file changed, 640 deletions(-)
 delete mode 100644 src/birkhoff_polytope_interfaces.jl

diff --git a/src/birkhoff_polytope_interfaces.jl b/src/birkhoff_polytope_interfaces.jl
deleted file mode 100644
index f76e5cf..0000000
--- a/src/birkhoff_polytope_interfaces.jl
+++ /dev/null
@@ -1,640 +0,0 @@
-"""
-	BirkhoffBLMO
-
-A bounded LMO for the Birkhoff polytope. This oracle computes an extreme point subject to  
-node-specific bounds on the integer variables.
-"""
-struct BirkhoffLMO
-    append_by_column::Bool
-    dim::Int
-    lower_bounds::Vector{Float64}
-    upper_bounds::Vector{Float64}
-    int_vars::Vector{Int}
-    fixed_to_one_rows::Vector{Int}
-    fixed_to_one_cols::Vector{Int}
-    index_map_rows::Vector{Int}
-    index_map_cols::Vector{Int}
-    atol::Float64
-    rtol::Float64
-end
-
-struct BosciaBirkhoffBLMO <: Boscia.BoundedLinearMinimizationOracle
-    birkhoffBLMO::BirkhoffLMO
-end
-
-struct FWBirkhoffLMO <: FrankWolfe.LinearMinimizationOracle
-    birkhoffLMO::BirkhoffLMO
-end
-
-# define the mixed-integer Birkhoff polytope
-function BirkhoffBLMO(dim, int_vars; append_by_column=true, atol=1e-6, rtol=1e-3)
-    birkhoffLMO = BirkhoffLMO(
-        append_by_column,
-        dim,
-        fill(0.0, length(int_vars)),
-        fill(1.0, length(int_vars)),
-        int_vars,
-        Int[],
-        Int[],
-        collect(1:dim),
-        collect(1:dim),
-        atol,
-        rtol,
-    )
-    return BosciaBirkhoffBLMO(birkhoffLMO)
-end
-
-function BirkhoffBLMO(dim; append_by_column=true, use_Boscia=true, atol=1e-6, rtol=1e-3)
-    # define the all-integer Birkhoff polytope
-    if use_Boscia
-        birkhoffLMO = BirkhoffLMO(
-            append_by_column,
-            dim,
-            fill(0.0, length(int_vars)),
-            fill(1.0, length(int_vars)),
-            collect(1:(dim^2)),
-            Int[],
-            Int[],
-            collect(1:dim),
-            collect(1:dim),
-            atol,
-            rtol,
-        )
-        return BosciaBirkhoffBLMO(birkhoffLMO)
-    else
-        # define the continous Birkhoff polytope
-        birkhoffLMO = BirkhoffLMO(
-            append_by_column,
-            dim,
-            [],
-            [],
-            [],
-            Int[],
-            Int[],
-            collect(1:dim),
-            collect(1:dim),
-            atol,
-            rtol,
-        )
-        return FWBirkhoffLMO(birkhoffLMO)
-    end
-
-end
-
-"""
-Computes the extreme point given an direction d, the current lower and upper bounds on the integer variables, and the set of integer variables.
-"""
-_compute_extreme_point_impl(blmo::BirkhoffLMO, d; kwargs...) = begin
-    n = blmo.dim
-
-    # Precompute index mapping to avoid repeated `findfirst` calls,
-    # which would be very costly inside the loop.
-    if length(blmo.int_vars) !== n^2
-        idx_map_ub = zeros(Int, n^2)
-        @inbounds for (c_idx, var) in enumerate(blmo.int_vars)
-            idx_map_ub[var] = c_idx
-        end
-    end
-
-    if size(d, 2) == 1
-        d = blmo.append_by_column ? reshape(d, (n, n)) : transpose(reshape(d, (n, n)))
-    end
-
-    fixed_to_one_rows = blmo.fixed_to_one_rows
-    fixed_to_one_cols = blmo.fixed_to_one_cols
-    index_map_rows = blmo.index_map_rows
-    index_map_cols = blmo.index_map_cols
-    int_vars = blmo.int_vars
-    ub = blmo.upper_bounds
-
-    nreduced = length(index_map_rows)
-    type = typeof(d[1, 1])
-    d2 = ones(Union{type,Missing}, nreduced, nreduced)
-    for j in 1:nreduced
-        col_orig = index_map_cols[j]
-        for i in 1:nreduced
-            row_orig = index_map_rows[i]
-            if blmo.append_by_column
-                orig_linear_idx = (col_orig - 1) * n + row_orig
-            else
-                orig_linear_idx = (row_orig - 1) * n + col_orig
-            end
-            # the problem can only be integer types,
-            # either all-integer or mixed-integer.
-            if orig_linear_idx in int_vars
-                idx = length(int_vars) < n^2 ? idx_map_ub[orig_linear_idx] : orig_linear_idx
-                # interdict arc when fixed to zero
-                if ub[idx] <= eps()
-                    if blmo.append_by_column
-                        d2[i, j] = missing
-                    else
-                        d2[j, i] = missing
-                    end
-                else
-                    if blmo.append_by_column
-                        d2[i, j] = d[row_orig, col_orig]
-                    else
-                        d2[j, i] = d[col_orig, row_orig]
-                    end
-                end
-            else
-                if blmo.append_by_column
-                    d2[i, j] = d[row_orig, col_orig]
-                else
-                    d2[j, i] = d[col_orig, row_orig]
-                end
-            end
-        end
-    end
-
-    m = SparseArrays.spzeros(n, n)
-    for (i, j) in zip(fixed_to_one_rows, fixed_to_one_cols)
-        m[i, j] = 1
-    end
-
-    res_mat = Hungarian.munkres(d2)
-    (rows, cols, vals) = SparseArrays.findnz(res_mat)
-    @inbounds for i in eachindex(cols)
-        m[index_map_rows[rows[i]], index_map_cols[cols[i]]] = (vals[i] == 2)
-    end
-
-    # if size(d, 2) == 1
-    m = if blmo.append_by_column
-        # Convert sparse matrix to sparse vector by columns
-        I, J, V = SparseArrays.findnz(m)
-        linear_indices = (J .- 1) .* n .+ I
-        SparseArrays.sparsevec(linear_indices, V, n^2)
-    else
-        # Convert sparse matrix to sparse vector by rows (transpose first)
-        mt = SparseArrays.sparse(LinearAlgebra.transpose(m))
-        I, J, V = SparseArrays.findnz(mt)
-        linear_indices = (J .- 1) .* n .+ I
-        SparseArrays.sparsevec(linear_indices, V, n^2)
-    end
-    # end
-    return m
-end
-
-"""
-Computes the extreme point given an direction d, the current lower and upper bounds on the integer variables, and the set of integer variables.
-"""
-_compute_inface_extreme_point_impl(blmo::BirkhoffLMO, direction, x; kwargs...) = begin
-    n = blmo.dim
-
-    if size(direction, 2) == 1
-        direction =
-            blmo.append_by_column ? reshape(direction, (n, n)) :
-            transpose(reshape(direction, (n, n)))
-
-        x = blmo.append_by_column ? reshape(x, (n, n)) : transpose(reshape(x, (n, n)))
-    end
-
-    # Precompute index mapping to avoid repeated `findfirst` calls,
-    # which would be very costly inside the loop.
-    if length(blmo.int_vars) !== n^2
-        idx_map_ub = zeros(Int, n^2)
-        @inbounds for (c_idx, var) in enumerate(blmo.int_vars)
-            idx_map_ub[var] = c_idx
-        end
-    end
-
-    fixed_to_one_rows = copy(blmo.fixed_to_one_rows)
-    fixed_to_one_cols = copy(blmo.fixed_to_one_cols)
-    index_map_rows = copy(blmo.index_map_rows)
-    index_map_cols = copy(blmo.index_map_cols)
-    int_vars = blmo.int_vars
-    ub = blmo.upper_bounds
-
-    nreduced = length(blmo.index_map_rows)
-
-    delete_index_map_rows = Int[]
-    delete_index_map_cols = Int[]
-    for j in 1:nreduced
-        for i in 1:nreduced
-            row_orig = index_map_rows[i]
-            col_orig = index_map_cols[j]
-            if x[row_orig, col_orig] >= 1-eps()
-                push!(fixed_to_one_rows, row_orig)
-                push!(fixed_to_one_cols, col_orig)
-
-                push!(delete_index_map_rows, i)
-                push!(delete_index_map_cols, j)
-            end
-        end
-    end
-
-    unique!(delete_index_map_rows)
-    unique!(delete_index_map_cols)
-    sort!(delete_index_map_rows)
-    sort!(delete_index_map_cols)
-    deleteat!(index_map_rows, delete_index_map_rows)
-    deleteat!(index_map_cols, delete_index_map_cols)
-
-    nreduced = length(index_map_rows)
-    type = typeof(direction[1, 1])
-    d2 = ones(Union{type,Missing}, nreduced, nreduced)
-    for j in 1:nreduced
-        for i in 1:nreduced
-            row_orig = index_map_rows[i]
-            col_orig = index_map_cols[j]
-            if blmo.append_by_column
-                orig_linear_idx = (col_orig-1)*n+row_orig
-            else
-                orig_linear_idx = (row_orig-1)*n+col_orig
-            end
-            # idx = findfirst(x -> x == orig_linear_idx, blmo.int_vars)
-
-            if x[row_orig, col_orig] <= eps()
-                if blmo.append_by_column
-                    d2[i, j] = missing
-                else
-                    d2[j, i] = missing
-                end
-            elseif orig_linear_idx in int_vars
-                idx = length(int_vars) < n^2 ? idx_map_ub[orig_linear_idx] : orig_linear_idx
-                # interdict arc when fixed to zero
-                if ub[idx] <= eps()
-                    if blmo.append_by_column
-                        d2[i, j] = missing
-                    else
-                        d2[j, i] = missing
-                    end
-                else
-                    if blmo.append_by_column
-                        d2[i, j] = direction[row_orig, col_orig]
-                    else
-                        d2[j, i] = direction[col_orig, row_orig]
-                    end
-                end
-            else
-                if blmo.append_by_column
-                    d2[i, j] = direction[row_orig, col_orig]
-                else
-                    d2[j, i] = direction[col_orig, row_orig]
-                end
-            end
-        end
-    end
-
-    m = SparseArrays.spzeros(n, n)
-    for (i, j) in zip(fixed_to_one_rows, fixed_to_one_cols)
-        m[i, j] = 1
-    end
-
-    res_mat = Hungarian.munkres(d2)
-    (rows, cols, vals) = SparseArrays.findnz(res_mat)
-    @inbounds for i in eachindex(cols)
-        m[index_map_rows[rows[i]], index_map_cols[cols[i]]] = (vals[i] == 2)
-    end
-
-    # if size(direction, 2) == 1
-    m = if blmo.append_by_column
-        # Convert sparse matrix to sparse vector by columns
-        I, J, V = SparseArrays.findnz(m)
-        linear_indices = (J .- 1) .* n .+ I
-        SparseArrays.sparsevec(linear_indices, V, n^2)
-    else
-        # Convert sparse matrix to sparse vector by rows (transpose first)
-        mt = SparseArrays.sparse(LinearAlgebra.transpose(m))
-        I, J, V = SparseArrays.findnz(mt)
-        linear_indices = (J .- 1) .* n .+ I
-        SparseArrays.sparsevec(linear_indices, V, n^2)
-    end
-    # end
-
-    return m
-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).
-"""
-_dicg_maximum_step_impl(blmo::BirkhoffLMO, direction, x; kwargs...) = begin
-    n = blmo.dim
-
-    T = promote_type(eltype(x), eltype(direction))
-    if size(direction, 2) == 1
-        direction =
-            blmo.append_by_column ? reshape(direction, (n, n)) :
-            transpose(reshape(direction, (n, n)))
-
-        x = blmo.append_by_column ? reshape(x, (n, n)) : transpose(reshape(x, (n, n)))
-    end
-
-    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
-
-_is_decomposition_invariant_oracle_impl(lmo::BirkhoffLMO) = true
-
-
-# Boscia API: method extensions for BirkhoffBLMO
-Boscia.compute_extreme_point(blmo::BosciaBirkhoffBLMO, d; kwargs...) =
-    _compute_extreme_point_impl(blmo.birkhoffBLMO, d; kwargs...)
-
-Boscia.compute_inface_extreme_point(blmo::BosciaBirkhoffBLMO, direction, x; kwargs...) =
-    _compute_inface_extreme_point_impl(blmo.birkhoffBLMO, direction, x; kwargs...)
-
-
-Boscia.dicg_maximum_step(blmo::BosciaBirkhoffBLMO, direction, x; kwargs...) =
-    _dicg_maximum_step_impl(blmo.birkhoffBLMO, direction, x; kwargs...)
-
-Boscia.is_decomposition_invariant_oracle(blmo::BosciaBirkhoffBLMO) =
-    _is_decomposition_invariant_oracle_impl(blmo.birkhoffBLMO)
-
-
-# FrankWolfe API: method extensions for BirkhoffLMO
-FrankWolfe.compute_extreme_point(lmo::FWBirkhoffLMO, d; kwargs...) =
-    _compute_extreme_point_impl(lmo.birkhoffLMO, d; kwargs...)
-
-FrankWolfe.compute_inface_extreme_point(lmo::FWBirkhoffLMO, direction, x; kwargs...) =
-    _compute_inface_extreme_point_impl(lmo.birkhoffLMO, direction, x; kwargs...)
-
-FrankWolfe.dicg_maximum_step(lmo::FWBirkhoffLMO, direction, x; kwargs...) =
-    _dicg_maximum_step_impl(lmo.birkhoffLMO, direction, x; kwargs...)
-
-FrankWolfe.is_decomposition_invariant_oracle(lmo::FWBirkhoffLMO) =
-    _is_decomposition_invariant_oracle_impl(lmo.birkhoffLMO)
-
-
-## Necessary for Boscia package
-
-"""
-The sum of each row and column has to be equal to 1.
-"""
-function Boscia.is_linear_feasible(lmo::BosciaBirkhoffBLMO, v::AbstractVector)
-    blmo = lmo.birkhoffBLMO
-    n = blmo.dim
-    for i in 1:n
-        # append by column ? column sum : row sum 
-        if !isapprox(sum(v[((i-1)*n+1):(i*n)]), 1.0, atol=1e-6, rtol=1e-3)
-            @debug "Column sum not 1: $(sum(v[((i-1)*n+1):(i*n)]))"
-            return false
-        end
-        # append by column ? row sum : column sum
-        if !isapprox(sum(v[i: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
-
-# Read global bounds from the problem.
-function Boscia.build_global_bounds(lmo::BosciaBirkhoffBLMO, integer_variables)
-    blmo = lmo.birkhoffBLMO
-    global_bounds = Boscia.IntegerBounds()
-    for (idx, int_var) in enumerate(blmo.int_vars)
-        push!(global_bounds, (int_var, blmo.lower_bounds[idx]), :greaterthan)
-        push!(global_bounds, (int_var, blmo.upper_bounds[idx]), :lessthan)
-    end
-    return global_bounds
-end
-
-# Get list of variables indices. 
-# If the problem has n variables, they are expected to contiguous and ordered from 1 to n.
-function Boscia.get_list_of_variables(lmo::BosciaBirkhoffBLMO)
-    blmo = lmo.birkhoffBLMO
-    n = blmo.dim^2
-    return n, collect(1:n)
-end
-
-# Get list of integer variables
-function Boscia.get_integer_variables(lmo::BosciaBirkhoffBLMO)
-    blmo = lmo.birkhoffBLMO
-    return blmo.int_vars
-end
-
-# Get the index of the integer variable the bound is working on.
-function Boscia.get_int_var(lmo::BosciaBirkhoffBLMO, cidx)
-    blmo = lmo.birkhoffBLMO
-    return blmo.int_vars[cidx]
-end
-
-# Get the list of lower bounds.
-function Boscia.get_lower_bound_list(lmo::BosciaBirkhoffBLMO)
-    blmo = lmo.birkhoffBLMO
-    return collect(1:length(blmo.lower_bounds))
-end
-
-# Get the list of upper bounds.
-function Boscia.get_upper_bound_list(lmo::BosciaBirkhoffBLMO)
-    blmo = lmo.birkhoffBLMO
-    return collect(1:length(blmo.upper_bounds))
-end
-
-# Read bound value for c_idx.
-function Boscia.get_bound(lmo::BosciaBirkhoffBLMO, c_idx, sense::Symbol)
-    blmo = lmo.birkhoffBLMO
-    if sense == :lessthan
-        return blmo.upper_bounds[c_idx]
-    elseif sense == :greaterthan
-        return blmo.lower_bounds[c_idx]
-    else
-        error("Allowed value for sense are :lessthan and :greaterthan!")
-    end
-end
-
-## Changing the bounds constraints.
-
-# Change the value of the bound c_idx.
-function Boscia.set_bound!(lmo::BosciaBirkhoffBLMO, c_idx, value, sense::Symbol)
-    blmo = lmo.birkhoffBLMO
-    if sense == :greaterthan
-        blmo.lower_bounds[c_idx] = value
-        if value == 1.0
-            n0 = blmo.dim
-            fixed_int_var = blmo.int_vars[c_idx]
-            # Convert linear index to (row, col) based on storage format
-            if blmo.append_by_column
-                j = ceil(Int, fixed_int_var / n0)  # column index
-                i = Int(fixed_int_var - n0 * (j - 1))  # row index
-            else
-                i = ceil(Int, fixed_int_var / n0)  # row index  
-                j = Int(fixed_int_var - n0 * (i - 1))  # column index
-            end
-            push!(blmo.fixed_to_one_rows, i)
-            push!(blmo.fixed_to_one_cols, j)
-        end
-    elseif sense == :lessthan
-        blmo.upper_bounds[c_idx] = value
-    else
-        error("Allowed values for sense are :lessthan and :greaterthan.")
-    end
-end
-
-# Delete bounds.
-function Boscia.delete_bounds!(lmo::BosciaBirkhoffBLMO, cons_delete)
-    blmo = lmo.birkhoffBLMO
-    for (d_idx, sense) in cons_delete
-        if sense == :greaterthan
-            blmo.lower_bounds[d_idx] = -Inf
-        else
-            blmo.upper_bounds[d_idx] = Inf
-        end
-    end
-
-    # sanity check 
-    check_feasibility(blmo, i::Int, j::Int) =
-        blmo.lower_bounds[blmo.append_by_column ? (j-1)*blmo.dim + i : (i-1)*blmo.dim + j] != 0.0
-
-    fixed_to_one_rows = blmo.fixed_to_one_rows
-    fixed_to_one_cols = blmo.fixed_to_one_cols
-
-    feasible_flags = check_feasibility.(Ref(blmo), fixed_to_one_rows, fixed_to_one_cols)
-    invalid_indices = findall(.!feasible_flags)
-
-    deleteat!(fixed_to_one_rows, invalid_indices)
-    deleteat!(fixed_to_one_cols, invalid_indices)
-
-    # remove the duplicate indices
-    pairs = collect(zip(blmo.fixed_to_one_rows, blmo.fixed_to_one_cols))
-    unique!(pairs)
-    resize!(blmo.fixed_to_one_rows, length(pairs))
-    resize!(blmo.fixed_to_one_cols, length(pairs))
-
-    blmo.fixed_to_one_rows .= first.(pairs)
-    blmo.fixed_to_one_cols .= last.(pairs)
-
-    nfixed = length(blmo.fixed_to_one_rows)
-    nreduced = blmo.dim - nfixed
-
-    # stores the indices of the original matrix that are still in the reduced matrix
-    index_map_rows = fill(1, nreduced)
-    index_map_cols = fill(1, nreduced)
-    idx_in_map_row = 1
-    idx_in_map_col = 1
-    for orig_idx in 1:blmo.dim
-        if orig_idx ∉ blmo.fixed_to_one_rows
-            index_map_rows[idx_in_map_row] = orig_idx
-            idx_in_map_row += 1
-        end
-        if orig_idx ∉ blmo.fixed_to_one_cols
-            index_map_cols[idx_in_map_col] = orig_idx
-            idx_in_map_col += 1
-        end
-    end
-
-    empty!(blmo.index_map_rows)
-    empty!(blmo.index_map_cols)
-    append!(blmo.index_map_rows, index_map_rows)
-    append!(blmo.index_map_cols, index_map_cols)
-
-    return blmo.modify_LMO = false
-end
-
-# Add bound constraint.
-function Boscia.add_bound_constraint!(lmo::BosciaBirkhoffBLMO, key, value, sense::Symbol)
-    blmo = lmo.birkhoffBLMO
-    idx = findfirst(x -> x == key, blmo.int_vars)
-    if sense == :greaterthan
-        blmo.lower_bounds[idx] = value
-
-    elseif sense == :lessthan
-        blmo.upper_bounds[idx] = value
-    else
-        error("Allowed value of sense are :lessthan and :greaterthan!")
-    end
-end
-
-## Checks
-
-# Check if the subject of the bound c_idx is an integer variable (recorded in int_vars).
-function Boscia.is_constraint_on_int_var(lmo::BosciaBirkhoffBLMO, c_idx, int_vars)
-    blmo = lmo.birkhoffBLMO
-    return blmo.int_vars[c_idx] in int_vars
-end
-
-# To check if there is bound for the variable in the global or node bounds.
-function Boscia.is_bound_in(lmo::BosciaBirkhoffBLMO, c_idx, bounds)
-    blmo = lmo.birkhoffBLMO
-    return haskey(bounds, blmo.int_vars[c_idx])
-end
-
-# Has variable an integer constraint?
-function Boscia.has_integer_constraint(lmo::BosciaBirkhoffBLMO, idx)
-    blmo = lmo.birkhoffBLMO
-    return idx in blmo.int_vars
-end
-
-## Optional
-
-function Boscia.check_feasibility(lmo::BosciaBirkhoffBLMO)
-    blmo = lmo.birkhoffBLMO
-    for (lb, ub) in zip(blmo.lower_bounds, blmo.upper_bounds)
-        if ub < lb
-            return Boscia.INFEASIBLE
-        end
-    end
-    # For double stochastic matrices, each row and column must sum to 1
-    # We check if the bounds allow for feasible assignments
-    n0 = blmo.dim
-    n = n0^2
-    int_vars = blmo.int_vars
-    # Initialize row and column bound tracking
-    row_min_sum = zeros(n0)  # minimum possible sum for each row
-    row_max_sum = zeros(n0)  # maximum possible sum for each row
-    col_min_sum = zeros(n0)  # minimum possible sum for each column
-    col_max_sum = zeros(n0)  # maximum possible sum for each column
-
-    rows_with_integer_variables = Int[]
-    cols_with_integer_variables = Int[]
-
-    # Process each integer variable
-    for idx in eachindex(int_vars)
-        var_idx = int_vars[idx]
-
-        # Convert linear index to (row, col) based on storage format
-        if blmo.append_by_column
-            j = ceil(Int, var_idx / n0)  # column index
-            i = Int(var_idx - n0 * (j - 1))  # row index
-        else
-            i = ceil(Int, var_idx / n0)  # row index  
-            j = Int(var_idx - n0 * (i - 1))  # column index
-        end
-
-        # Add bounds to row and column sums
-        row_min_sum[i] += blmo.lower_bounds[idx]
-        row_max_sum[i] += blmo.upper_bounds[idx]
-        col_min_sum[j] += blmo.lower_bounds[idx]
-        col_max_sum[j] += blmo.upper_bounds[idx]
-
-        push!(rows_with_integer_variables, i)
-        push!(cols_with_integer_variables, j)
-    end
-
-    # Check feasibility: each row and column must be able to sum to exactly 1
-    for i in 1:n0
-        if i in rows_with_integer_variables
-            # Check row sum constraints
-            if row_min_sum[i] > 1 + eps() || row_max_sum[i] < 1 - eps()
-                return Boscia.INFEASIBLE
-            end
-        end
-
-        if i in cols_with_integer_variables
-            # Check column sum constraints  
-            if col_min_sum[i] > 1 + eps() || col_max_sum[i] < 1 - eps()
-                return Boscia.INFEASIBLE
-            end
-        end
-    end
-
-    return Boscia.OPTIMAL
-end
-

From f1dc3ee586d549cd166e1b88191509a89cde63dc Mon Sep 17 00:00:00 2001
From: Wenjie Xiao <116566662+WenjieXiao-2022@users.noreply.github.com>
Date: Wed, 10 Sep 2025 21:25:34 +0200
Subject: [PATCH 15/21] Support mixed-integer and continous problem

---
 src/birkhoff_polytope.jl | 364 ++++++++++++++++++++++++---------------
 1 file changed, 225 insertions(+), 139 deletions(-)

diff --git a/src/birkhoff_polytope.jl b/src/birkhoff_polytope.jl
index d3df39a..4cf30bc 100644
--- a/src/birkhoff_polytope.jl
+++ b/src/birkhoff_polytope.jl
@@ -1,75 +1,121 @@
 """
-	BirkhoffBLMO
+	BirkhoffLMO
 
 A bounded LMO for the Birkhoff polytope. This oracle computes an extreme point subject to  
 node-specific bounds on the integer variables.
 """
-struct BirkhoffBLMO <: Boscia.BoundedLinearMinimizationOracle
+mutable struct BirkhoffLMO <: Boscia.BoundedLinearMinimizationOracle
     append_by_column::Bool
     dim::Int
+    lower_bounds::Vector{Float64}
+    upper_bounds::Vector{Float64}
     int_vars::Vector{Int}
     fixed_to_one_rows::Vector{Int}
     fixed_to_one_cols::Vector{Int}
     index_map_rows::Vector{Int}
     index_map_cols::Vector{Int}
+    updated_lmo::Bool
     atol::Float64
     rtol::Float64
+    boscia_use::Bool
 end
 
-BirkhoffBLMO(dim, lower_bounds, upper_bounds, int_vars; append_by_column = true, atol=1e-6, rtol=1e-3) =
-    BirkhoffBLMO(
-        append_by_column,
-        dim,
-        fill(0.0, length(int_vars)),
-        fill(1.0, length(int_vars)),
-        int_vars,
-        Int[],
-        Int[],
-        collect(1:dim),
-        collect(1:dim),
-        atol,
-        rtol,
-    )
+# return an integer-type BirkhoffLMO for Boscia
+BirkhoffLMO(dim, int_vars; append_by_column=true, atol=1e-6, rtol=1e-3) = BirkhoffLMO(
+    append_by_column,
+    dim,
+    fill(0.0, length(int_vars)),
+    fill(1.0, length(int_vars)),
+    int_vars,
+    Int[],
+    Int[],
+    collect(1:dim),
+    collect(1:dim),
+    false,
+    atol,
+    rtol,
+    true,
+)
+
+# return a continuous BirkhoffLMO for FrankWolfe
+BirkhoffLMO(dim; append_by_column=true, atol=1e-6, rtol=1e-3) = BirkhoffLMO(
+    append_by_column,
+    dim,
+    Float64[],
+    Float64[],
+    Int[],
+    Int[],
+    Int[],
+    collect(1:dim),
+    collect(1:dim),
+    false,
+    atol,
+    rtol,
+    false,
+)
 
 ## Necessary
 
 """
 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.compute_extreme_point(blmo::BirkhoffBLMO, d; kwargs...)
-    n = blmo.dim
+
+function Boscia.compute_extreme_point(lmo::BirkhoffLMO, d; kwargs...)
+    n = lmo.dim
+
+    # Precompute index mapping to avoid repeated `findfirst` calls,
+    # which would be very costly inside the loop.
+    if length(lmo.int_vars) != n^2
+        idx_map_ub = zeros(Int, n^2)
+        @inbounds for (c_idx, var) in enumerate(lmo.int_vars)
+            idx_map_ub[var] = c_idx
+        end
+    end
 
     if size(d, 2) == 1
-        d = blmo.append_by_column ? reshape(d, (n, n)) : transpose(reshape(d, (n, n)))
+        d = lmo.append_by_column ? reshape(d, (n, n)) : transpose(reshape(d, (n, n)))
     end
 
-    fixed_to_one_rows = blmo.fixed_to_one_rows
-    fixed_to_one_cols = blmo.fixed_to_one_cols
-    index_map_rows = blmo.index_map_rows
-    index_map_cols = blmo.index_map_cols
-    ub = blmo.upper_bounds
+    fixed_to_one_rows = lmo.fixed_to_one_rows
+    fixed_to_one_cols = lmo.fixed_to_one_cols
+    index_map_rows = lmo.index_map_rows
+    index_map_cols = lmo.index_map_cols
+    int_vars = lmo.int_vars
+    ub = lmo.upper_bounds
 
     nreduced = length(index_map_rows)
     type = typeof(d[1, 1])
     d2 = ones(Union{type,Missing}, nreduced, nreduced)
-    for j = 1:nreduced
+
+    for j in 1:nreduced
         col_orig = index_map_cols[j]
-        for i = 1:nreduced
+        for i in 1:nreduced
             row_orig = index_map_rows[i]
-            if blmo.append_by_column
+            if lmo.append_by_column
                 orig_linear_idx = (col_orig - 1) * n + row_orig
             else
                 orig_linear_idx = (row_orig - 1) * n + col_orig
             end
-            # interdict arc when fixed to zero
-            if ub[orig_linear_idx] <= eps()
-                if blmo.append_by_column
-                    d2[i, j] = missing
+            # the problem can only be integer types,
+            # either all-integer or mixed-integer.
+            if length(int_vars) == n^2 || idx_map_ub[orig_linear_idx] != 0
+                idx = length(int_vars) < n^2 ? idx_map_ub[orig_linear_idx] : orig_linear_idx
+                # interdict arc when fixed to zero
+                if ub[idx] <= eps()
+                    if lmo.append_by_column
+                        d2[i, j] = missing
+                    else
+                        d2[j, i] = missing
+                    end
                 else
-                    d2[j, i] = missing
+                    if lmo.append_by_column
+                        d2[i, j] = d[row_orig, col_orig]
+                    else
+                        d2[j, i] = d[col_orig, row_orig]
+                    end
                 end
             else
-                if blmo.append_by_column
+                if lmo.append_by_column
                     d2[i, j] = d[row_orig, col_orig]
                 else
                     d2[j, i] = d[col_orig, row_orig]
@@ -89,53 +135,63 @@ function Boscia.compute_extreme_point(blmo::BirkhoffBLMO, d; kwargs...)
         m[index_map_rows[rows[i]], index_map_cols[cols[i]]] = (vals[i] == 2)
     end
 
-    m = if blmo.append_by_column
-        # Convert sparse matrix to sparse vector by columns
-        I, J, V = SparseArrays.findnz(m)
-        linear_indices = (J .- 1) .* n .+ I
-        SparseArrays.sparsevec(linear_indices, V, n^2)
-    else
-        # Convert sparse matrix to sparse vector by rows (transpose first)
-        mt = SparseArrays.sparse(LinearAlgebra.transpose(m))
-        I, J, V = SparseArrays.findnz(mt)
-        linear_indices = (J .- 1) .* n .+ I
-        SparseArrays.sparsevec(linear_indices, V, n^2)
+    if size(d, 2) == 1 || lmo.boscia_use
+        m = if lmo.append_by_column
+            # Convert sparse matrix to sparse vector by columns
+            I, J, V = SparseArrays.findnz(m)
+            linear_indices = (J .- 1) .* n .+ I
+            SparseArrays.sparsevec(linear_indices, V, n^2)
+        else
+            # Convert sparse matrix to sparse vector by rows (transpose first)
+            mt = SparseArrays.sparse(LinearAlgebra.transpose(m))
+            I, J, V = SparseArrays.findnz(mt)
+            linear_indices = (J .- 1) .* n .+ I
+            SparseArrays.sparsevec(linear_indices, V, n^2)
+        end
     end
     return m
 end
 
+
 """
 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.compute_inface_extreme_point(blmo::BirkhoffBLMO, direction, x; kwargs...)
-    n = blmo.dim
+function Boscia.compute_inface_extreme_point(lmo::BirkhoffLMO, direction, x; kwargs...)
+    n = lmo.dim
 
     if size(direction, 2) == 1
         direction =
-            blmo.append_by_column ? reshape(direction, (n, n)) :
+            lmo.append_by_column ? reshape(direction, (n, n)) :
             transpose(reshape(direction, (n, n)))
-    end
 
-    if size(x, 2) == 1
-        x = blmo.append_by_column ? reshape(x, (n, n)) : transpose(reshape(x, (n, n)))
+        x = lmo.append_by_column ? reshape(x, (n, n)) : transpose(reshape(x, (n, n)))
     end
 
+    # Precompute index mapping to avoid repeated `findfirst` calls,
+    # which would be very costly inside the loop.
+    if length(lmo.int_vars) !== n^2
+        idx_map_ub = zeros(Int, n^2)
+        @inbounds for (c_idx, var) in enumerate(lmo.int_vars)
+            idx_map_ub[var] = c_idx
+        end
+    end
 
-    fixed_to_one_rows = copy(blmo.fixed_to_one_rows)
-    fixed_to_one_cols = copy(blmo.fixed_to_one_cols)
-    index_map_rows = copy(blmo.index_map_rows)
-    index_map_cols = copy(blmo.index_map_cols)
-    ub = blmo.upper_bounds
+    fixed_to_one_rows = copy(lmo.fixed_to_one_rows)
+    fixed_to_one_cols = copy(lmo.fixed_to_one_cols)
+    index_map_rows = copy(lmo.index_map_rows)
+    index_map_cols = copy(lmo.index_map_cols)
+    int_vars = lmo.int_vars
+    ub = lmo.upper_bounds
 
-    nreduced = length(blmo.index_map_rows)
+    nreduced = length(lmo.index_map_rows)
 
     delete_index_map_rows = Int[]
     delete_index_map_cols = Int[]
-    for j = 1:nreduced
-        for i = 1:nreduced
+    for j in 1:nreduced
+        for i in 1:nreduced
             row_orig = index_map_rows[i]
             col_orig = index_map_cols[j]
-            if x[row_orig, col_orig] >= 1 - eps()
+            if x[row_orig, col_orig] >= 1-eps()
                 push!(fixed_to_one_rows, row_orig)
                 push!(fixed_to_one_cols, col_orig)
 
@@ -155,24 +211,39 @@ function Boscia.compute_inface_extreme_point(blmo::BirkhoffBLMO, direction, x; k
     nreduced = length(index_map_rows)
     type = typeof(direction[1, 1])
     d2 = ones(Union{type,Missing}, nreduced, nreduced)
-    for j = 1:nreduced
-        for i = 1:nreduced
+    for j in 1:nreduced
+        for i in 1:nreduced
             row_orig = index_map_rows[i]
             col_orig = index_map_cols[j]
-            if blmo.append_by_column
-                orig_linear_idx = (col_orig - 1) * n + row_orig
+            if lmo.append_by_column
+                orig_linear_idx = (col_orig-1)*n+row_orig
             else
-                orig_linear_idx = (row_orig - 1) * n + col_orig
+                orig_linear_idx = (row_orig-1)*n+col_orig
             end
-            # interdict arc when fixed to zero
-            if ub[orig_linear_idx] <= eps() || x[row_orig, col_orig] <= eps()
-                if blmo.append_by_column
+            if x[row_orig, col_orig] <= eps()
+                if lmo.append_by_column
                     d2[i, j] = missing
                 else
                     d2[j, i] = missing
                 end
+            elseif length(int_vars) == n^2 || idx_map_ub[orig_linear_idx] != 0
+                idx = length(int_vars) < n^2 ? idx_map_ub[orig_linear_idx] : orig_linear_idx
+                # interdict arc when fixed to zero
+                if ub[idx] <= eps()
+                    if lmo.append_by_column
+                        d2[i, j] = missing
+                    else
+                        d2[j, i] = missing
+                    end
+                else
+                    if lmo.append_by_column
+                        d2[i, j] = direction[row_orig, col_orig]
+                    else
+                        d2[j, i] = direction[col_orig, row_orig]
+                    end
+                end
             else
-                if blmo.append_by_column
+                if lmo.append_by_column
                     d2[i, j] = direction[row_orig, col_orig]
                 else
                     d2[j, i] = direction[col_orig, row_orig]
@@ -192,17 +263,19 @@ function Boscia.compute_inface_extreme_point(blmo::BirkhoffBLMO, direction, x; k
         m[index_map_rows[rows[i]], index_map_cols[cols[i]]] = (vals[i] == 2)
     end
 
-    m = if blmo.append_by_column
-        # Convert sparse matrix to sparse vector by columns
-        I, J, V = SparseArrays.findnz(m)
-        linear_indices = (J .- 1) .* n .+ I
-        SparseArrays.sparsevec(linear_indices, V, n^2)
-    else
-        # Convert sparse matrix to sparse vector by rows (transpose first)
-        mt = SparseArrays.sparse(LinearAlgebra.transpose(m))
-        I, J, V = SparseArrays.findnz(mt)
-        linear_indices = (J .- 1) .* n .+ I
-        SparseArrays.sparsevec(linear_indices, V, n^2)
+    if size(direction, 2) == 1 || lmo.boscia_use
+        m = if lmo.append_by_column
+            # Convert sparse matrix to sparse vector by columns
+            I, J, V = SparseArrays.findnz(m)
+            linear_indices = (J .- 1) .* n .+ I
+            SparseArrays.sparsevec(linear_indices, V, n^2)
+        else
+            # Convert sparse matrix to sparse vector by rows (transpose first)
+            mt = SparseArrays.sparse(LinearAlgebra.transpose(m))
+            I, J, V = SparseArrays.findnz(mt)
+            linear_indices = (J .- 1) .* n .+ I
+            SparseArrays.sparsevec(linear_indices, V, n^2)
+        end
     end
 
     return m
@@ -212,14 +285,15 @@ 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.dicg_maximum_step(blmo::BirkhoffBLMO, direction, x; kwargs...)
-    n = blmo.dim
-	T = promote_type(eltype(x), eltype(direction))
+function Boscia.dicg_maximum_step(lmo::BirkhoffLMO, direction, x; kwargs...)
+    n = lmo.dim
+    T = promote_type(eltype(x), eltype(direction))
     if size(direction, 2) == 1
         direction =
-            blmo.append_by_column ? reshape(direction, (n, n)) :
+            lmo.append_by_column ? reshape(direction, (n, n)) :
             transpose(reshape(direction, (n, n)))
-        x = blmo.append_by_column ? reshape(x, (n, n)) : transpose(reshape(x, (n, n)))
+
+        x = lmo.append_by_column ? reshape(x, (n, n)) : transpose(reshape(x, (n, n)))
     end
 
     gamma_max = one(T)
@@ -239,25 +313,26 @@ function Boscia.dicg_maximum_step(blmo::BirkhoffBLMO, direction, x; kwargs...)
         end
     end
     return gamma_max
+
 end
 
-function Boscia.is_decomposition_invariant_oracle(blmo::BirkhoffBLMO)
+function Boscia.is_decomposition_invariant_oracle(lmo::BirkhoffLMO)
     return true
 end
 
 """
 The sum of each row and column has to be equal to 1.
 """
-function Boscia.is_linear_feasible(blmo::BirkhoffBLMO, v::AbstractVector)
+function Boscia.is_linear_feasible(blmo::BirkhoffLMO, v::AbstractVector)
     n = blmo.dim
-    for i = 1:n
+    for i in 1:n
         # append by column ? column sum : row sum 
-        if !isapprox(sum(v[((i-1)*n+1):(i*n)]), 1.0, atol = 1e-6, rtol = 1e-3)
+        if !isapprox(sum(v[((i-1)*n+1):(i*n)]), 1.0, atol=1e-6, rtol=1e-3)
             @debug "Column sum not 1: $(sum(v[((i-1)*n+1):(i*n)]))"
             return false
         end
         # append by column ? row sum : column sum
-        if !isapprox(sum(v[i:n:(n^2)]), 1.0, atol = 1e-6, rtol = 1e-3)
+        if !isapprox(sum(v[i: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
@@ -266,7 +341,7 @@ function Boscia.is_linear_feasible(blmo::BirkhoffBLMO, v::AbstractVector)
 end
 
 # Read global bounds from the problem.
-function Boscia.build_global_bounds(blmo::BirkhoffBLMO, integer_variables)
+function Boscia.build_global_bounds(blmo::BirkhoffLMO, integer_variables)
     global_bounds = Boscia.IntegerBounds()
     for (idx, int_var) in enumerate(blmo.int_vars)
         push!(global_bounds, (int_var, blmo.lower_bounds[idx]), :greaterthan)
@@ -277,33 +352,33 @@ end
 
 # Get list of variables indices. 
 # If the problem has n variables, they are expected to contiguous and ordered from 1 to n.
-function Boscia.get_list_of_variables(blmo::BirkhoffBLMO)
+function Boscia.get_list_of_variables(blmo::BirkhoffLMO)
     n = blmo.dim^2
     return n, collect(1:n)
 end
 
 # Get list of integer variables
-function Boscia.get_integer_variables(blmo::BirkhoffBLMO)
+function Boscia.get_integer_variables(blmo::BirkhoffLMO)
     return blmo.int_vars
 end
 
 # Get the index of the integer variable the bound is working on.
-function Boscia.get_int_var(blmo::BirkhoffBLMO, cidx)
+function Boscia.get_int_var(blmo::BirkhoffLMO, cidx)
     return blmo.int_vars[cidx]
 end
 
 # Get the list of lower bounds.
-function Boscia.get_lower_bound_list(blmo::BirkhoffBLMO)
+function Boscia.get_lower_bound_list(blmo::BirkhoffLMO)
     return collect(1:length(blmo.lower_bounds))
 end
 
 # Get the list of upper bounds.
-function Boscia.get_upper_bound_list(blmo::BirkhoffBLMO)
+function Boscia.get_upper_bound_list(blmo::BirkhoffLMO)
     return collect(1:length(blmo.upper_bounds))
 end
 
 # Read bound value for c_idx.
-function Boscia.get_bound(blmo::BirkhoffBLMO, c_idx, sense::Symbol)
+function Boscia.get_bound(blmo::BirkhoffLMO, c_idx, sense::Symbol)
     if sense == :lessthan
         return blmo.upper_bounds[c_idx]
     elseif sense == :greaterthan
@@ -316,10 +391,16 @@ end
 ## Changing the bounds constraints.
 
 # Change the value of the bound c_idx.
-function Boscia.set_bound!(blmo::BirkhoffBLMO, c_idx, value, sense::Symbol)
+function Boscia.set_bound!(blmo::BirkhoffLMO, c_idx, value, sense::Symbol)
+    # Reset the lmo if necessary
+    if !blmo.updated_lmo
+        empty!(blmo.fixed_to_one_rows)
+        empty!(blmo.fixed_to_one_cols)
+        blmo.updated_lmo = true
+    end
     if sense == :greaterthan
         blmo.lower_bounds[c_idx] = value
-		if value == 1.0
+        if value == 1.0
             n0 = blmo.dim
             fixed_int_var = blmo.int_vars[c_idx]
             # Convert linear index to (row, col) based on storage format
@@ -341,41 +422,19 @@ function Boscia.set_bound!(blmo::BirkhoffBLMO, c_idx, value, sense::Symbol)
 end
 
 # Delete bounds.
-function Boscia.delete_bounds!(blmo::BirkhoffBLMO, cons_delete)
+function Boscia.delete_bounds!(blmo::BirkhoffLMO, cons_delete)
     for (d_idx, sense) in cons_delete
         if sense == :greaterthan
-            blmo.lower_bounds[d_idx] = -Inf
+            blmo.lower_bounds[d_idx] = 0.0
         else
-            blmo.upper_bounds[d_idx] = Inf
+            blmo.upper_bounds[d_idx] = 1.0
         end
     end
 
-	# sanity check
-    check_feasibility(blmo, i::Int, j::Int) =
-        blmo.lower_bounds[blmo.append_by_column ? (j-1)*blmo.dim + i : (i-1)*blmo.dim + j] != 0.0
-
-    fixed_to_one_rows = blmo.fixed_to_one_rows
-    fixed_to_one_cols = blmo.fixed_to_one_cols
-
-    feasible_flags  = check_feasibility.(Ref(blmo), fixed_to_one_rows, fixed_to_one_cols)
-    invalid_indices = findall(.!feasible_flags)
-
-    deleteat!(fixed_to_one_rows, invalid_indices)
-    deleteat!(fixed_to_one_cols, invalid_indices)
-
-    # remove the duplicate indices
-    pairs = collect(zip(blmo.fixed_to_one_rows, blmo.fixed_to_one_cols))
-    unique!(pairs)
-    resize!(blmo.fixed_to_one_rows, length(pairs))
-    resize!(blmo.fixed_to_one_cols, length(pairs))
-
-    blmo.fixed_to_one_rows .= first.(pairs)
-    blmo.fixed_to_one_cols .= last.(pairs)
-
     nfixed = length(blmo.fixed_to_one_rows)
     nreduced = blmo.dim - nfixed
 
-    # stores the indices of the original matrix that are still in the reduced matrix
+    # Store the indices of the original matrix that are still in the reduced matrix
     index_map_rows = fill(1, nreduced)
     index_map_cols = fill(1, nreduced)
     idx_in_map_row = 1
@@ -394,13 +453,30 @@ function Boscia.delete_bounds!(blmo::BirkhoffBLMO, cons_delete)
     empty!(blmo.index_map_rows)
     empty!(blmo.index_map_cols)
     append!(blmo.index_map_rows, index_map_rows)
+    append!(blmo.index_map_cols, index_map_cols)
+    blmo.updated_lmo = false
+    return true
 end
 
 # Add bound constraint.
-function Boscia.add_bound_constraint!(blmo::BirkhoffBLMO, key, value, sense::Symbol)
+function Boscia.add_bound_constraint!(blmo::BirkhoffLMO, key, value, sense::Symbol)
     idx = findfirst(x -> x == key, blmo.int_vars)
     if sense == :greaterthan
         blmo.lower_bounds[idx] = value
+        if value == 1.0
+            n0 = blmo.dim
+            fixed_int_var = blmo.int_vars[c_idx]
+            # Convert linear index to (row, col) based on storage format
+            if blmo.append_by_column
+                j = ceil(Int, fixed_int_var / n0)  # column index
+                i = Int(fixed_int_var - n0 * (j - 1))  # row index
+            else
+                i = ceil(Int, fixed_int_var / n0)  # row index  
+                j = Int(fixed_int_var - n0 * (i - 1))  # column index
+            end
+            push!(blmo.fixed_to_one_rows, i)
+            push!(blmo.fixed_to_one_cols, j)
+        end
     elseif sense == :lessthan
         blmo.upper_bounds[idx] = value
     else
@@ -411,23 +487,23 @@ end
 ## Checks
 
 # Check if the subject of the bound c_idx is an integer variable (recorded in int_vars).
-function Boscia.is_constraint_on_int_var(blmo::BirkhoffBLMO, c_idx, int_vars)
+function Boscia.is_constraint_on_int_var(blmo::BirkhoffLMO, c_idx, int_vars)
     return blmo.int_vars[c_idx] in int_vars
 end
 
 # To check if there is bound for the variable in the global or node bounds.
-function Boscia.is_bound_in(blmo::BirkhoffBLMO, c_idx, bounds)
+function Boscia.is_bound_in(blmo::BirkhoffLMO, c_idx, bounds)
     return haskey(bounds, blmo.int_vars[c_idx])
 end
 
 # Has variable an integer constraint?
-function Boscia.has_integer_constraint(blmo::BirkhoffBLMO, idx)
+function Boscia.has_integer_constraint(blmo::BirkhoffLMO, idx)
     return idx in blmo.int_vars
 end
 
 ## Optional
 
-function Boscia.check_feasibility(blmo::BirkhoffBLMO)
+function Boscia.check_feasibility(blmo::BirkhoffLMO)
     for (lb, ub) in zip(blmo.lower_bounds, blmo.upper_bounds)
         if ub < lb
             return Boscia.INFEASIBLE
@@ -437,15 +513,19 @@ function Boscia.check_feasibility(blmo::BirkhoffBLMO)
     # We check if the bounds allow for feasible assignments
     n0 = blmo.dim
     n = n0^2
+    int_vars = blmo.int_vars
     # Initialize row and column bound tracking
-    row_min_sum = zeros(n)  # minimum possible sum for each row
-    row_max_sum = zeros(n)  # maximum possible sum for each row
-    col_min_sum = zeros(n)  # minimum possible sum for each column
-    col_max_sum = zeros(n)  # maximum possible sum for each column
+    row_min_sum = zeros(n0)  # minimum possible sum for each row
+    row_max_sum = zeros(n0)  # maximum possible sum for each row
+    col_min_sum = zeros(n0)  # minimum possible sum for each column
+    col_max_sum = zeros(n0)  # maximum possible sum for each column
+
+    rows_with_integer_variables = Int[]
+    cols_with_integer_variables = Int[]
 
     # Process each integer variable
-    for idx in eachindex(blmo.int_vars)
-        var_idx = blmo.int_vars[idx]
+    for idx in eachindex(int_vars)
+        var_idx = int_vars[idx]
 
         # Convert linear index to (row, col) based on storage format
         if blmo.append_by_column
@@ -461,17 +541,23 @@ function Boscia.check_feasibility(blmo::BirkhoffBLMO)
         row_max_sum[i] += blmo.upper_bounds[idx]
         col_min_sum[j] += blmo.lower_bounds[idx]
         col_max_sum[j] += blmo.upper_bounds[idx]
+
+        push!(rows_with_integer_variables, i)
+        push!(cols_with_integer_variables, j)
     end
 
+    rows_with_integer_variables = unique(rows_with_integer_variables)
+    cols_with_integer_variables = unique(cols_with_integer_variables)
+
     # Check feasibility: each row and column must be able to sum to exactly 1
-    for i = 1:n0
-        # Check row sum constraints
+    for i in rows_with_integer_variables
         if row_min_sum[i] > 1 + eps() || row_max_sum[i] < 1 - eps()
             return Boscia.INFEASIBLE
         end
+    end
 
-        # Check column sum constraints  
-        if col_min_sum[i] > 1 + eps() || col_max_sum[i] < 1 - eps()
+    for j in cols_with_integer_variables
+        if col_min_sum[j] > 1 + eps() || col_max_sum[j] < 1 - eps()
             return Boscia.INFEASIBLE
         end
     end

From 84f767c4aaf85438bec0bbccad5edf5da2d202e3 Mon Sep 17 00:00:00 2001
From: Wenjie Xiao <116566662+WenjieXiao-2022@users.noreply.github.com>
Date: Wed, 10 Sep 2025 21:29:51 +0200
Subject: [PATCH 16/21] format

---
 src/birkhoff_polytope.jl | 1 -
 1 file changed, 1 deletion(-)

diff --git a/src/birkhoff_polytope.jl b/src/birkhoff_polytope.jl
index 4cf30bc..da9b9d4 100644
--- a/src/birkhoff_polytope.jl
+++ b/src/birkhoff_polytope.jl
@@ -59,7 +59,6 @@ BirkhoffLMO(dim; append_by_column=true, atol=1e-6, rtol=1e-3) = BirkhoffLMO(
 """
 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.compute_extreme_point(lmo::BirkhoffLMO, d; kwargs...)
     n = lmo.dim
 

From 62bf94a1b561e31ca2728ec2ae28241c4d049358 Mon Sep 17 00:00:00 2001
From: Wenjie Xiao <116566662+WenjieXiao-2022@users.noreply.github.com>
Date: Mon, 15 Sep 2025 10:27:49 +0200
Subject: [PATCH 17/21] support mixed-integer and continous problem

---
 src/birkhoff_polytope.jl | 78 ++++++++++++++++++++++++++++------------
 1 file changed, 56 insertions(+), 22 deletions(-)

diff --git a/src/birkhoff_polytope.jl b/src/birkhoff_polytope.jl
index da9b9d4..c696db8 100644
--- a/src/birkhoff_polytope.jl
+++ b/src/birkhoff_polytope.jl
@@ -20,7 +20,7 @@ mutable struct BirkhoffLMO <: Boscia.BoundedLinearMinimizationOracle
     boscia_use::Bool
 end
 
-# return an integer-type BirkhoffLMO for Boscia
+# return an integer-type BirkhoffLMO
 BirkhoffLMO(dim, int_vars; append_by_column=true, atol=1e-6, rtol=1e-3) = BirkhoffLMO(
     append_by_column,
     dim,
@@ -31,13 +31,13 @@ BirkhoffLMO(dim, int_vars; append_by_column=true, atol=1e-6, rtol=1e-3) = Birkho
     Int[],
     collect(1:dim),
     collect(1:dim),
-    false,
+    true,
     atol,
     rtol,
     true,
 )
 
-# return a continuous BirkhoffLMO for FrankWolfe
+# return a continuous BirkhoffLMO
 BirkhoffLMO(dim; append_by_column=true, atol=1e-6, rtol=1e-3) = BirkhoffLMO(
     append_by_column,
     dim,
@@ -48,7 +48,7 @@ BirkhoffLMO(dim; append_by_column=true, atol=1e-6, rtol=1e-3) = BirkhoffLMO(
     Int[],
     collect(1:dim),
     collect(1:dim),
-    false,
+    true,
     atol,
     rtol,
     false,
@@ -59,18 +59,10 @@ BirkhoffLMO(dim; append_by_column=true, atol=1e-6, rtol=1e-3) = BirkhoffLMO(
 """
 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.compute_extreme_point(lmo::BirkhoffLMO, d; kwargs...)
     n = lmo.dim
 
-    # Precompute index mapping to avoid repeated `findfirst` calls,
-    # which would be very costly inside the loop.
-    if length(lmo.int_vars) != n^2
-        idx_map_ub = zeros(Int, n^2)
-        @inbounds for (c_idx, var) in enumerate(lmo.int_vars)
-            idx_map_ub[var] = c_idx
-        end
-    end
-
     if size(d, 2) == 1
         d = lmo.append_by_column ? reshape(d, (n, n)) : transpose(reshape(d, (n, n)))
     end
@@ -82,6 +74,16 @@ function Boscia.compute_extreme_point(lmo::BirkhoffLMO, d; kwargs...)
     int_vars = lmo.int_vars
     ub = lmo.upper_bounds
 
+    is_full_integer = length(int_vars) == n^2 ? true : false
+    # Precompute index mapping to avoid repeated `findfirst` calls,
+    # which would be very costly inside the loop.
+    if !is_full_integer
+        idx_map_ub = zeros(Int, n^2)
+        @inbounds for (c_idx, var) in enumerate(lmo.int_vars)
+            idx_map_ub[var] = c_idx
+        end
+    end
+
     nreduced = length(index_map_rows)
     type = typeof(d[1, 1])
     d2 = ones(Union{type,Missing}, nreduced, nreduced)
@@ -96,9 +98,9 @@ function Boscia.compute_extreme_point(lmo::BirkhoffLMO, d; kwargs...)
                 orig_linear_idx = (row_orig - 1) * n + col_orig
             end
             # the problem can only be integer types,
-            # either all-integer or mixed-integer.
-            if length(int_vars) == n^2 || idx_map_ub[orig_linear_idx] != 0
-                idx = length(int_vars) < n^2 ? idx_map_ub[orig_linear_idx] : orig_linear_idx
+            # either full-integer or mixed-integer.
+            if is_full_integer || idx_map_ub[orig_linear_idx] != 0
+                idx = is_full_integer ? orig_linear_idx : idx_map_ub[orig_linear_idx]
                 # interdict arc when fixed to zero
                 if ub[idx] <= eps()
                     if lmo.append_by_column
@@ -186,7 +188,7 @@ function Boscia.compute_inface_extreme_point(lmo::BirkhoffLMO, direction, x; kwa
 
     delete_index_map_rows = Int[]
     delete_index_map_cols = Int[]
-    for j in 1:nreduced
+    delete_reducedUB = for j in 1:nreduced
         for i in 1:nreduced
             row_orig = index_map_rows[i]
             col_orig = index_map_cols[j]
@@ -211,9 +213,9 @@ function Boscia.compute_inface_extreme_point(lmo::BirkhoffLMO, direction, x; kwa
     type = typeof(direction[1, 1])
     d2 = ones(Union{type,Missing}, nreduced, nreduced)
     for j in 1:nreduced
+        col_orig = index_map_cols[j]
         for i in 1:nreduced
             row_orig = index_map_rows[i]
-            col_orig = index_map_cols[j]
             if lmo.append_by_column
                 orig_linear_idx = (col_orig-1)*n+row_orig
             else
@@ -225,6 +227,8 @@ function Boscia.compute_inface_extreme_point(lmo::BirkhoffLMO, direction, x; kwa
                 else
                     d2[j, i] = missing
                 end
+                # the problem can only be integer types,
+                # either full-integer or mixed-integer.
             elseif length(int_vars) == n^2 || idx_map_ub[orig_linear_idx] != 0
                 idx = length(int_vars) < n^2 ? idx_map_ub[orig_linear_idx] : orig_linear_idx
                 # interdict arc when fixed to zero
@@ -276,7 +280,6 @@ function Boscia.compute_inface_extreme_point(lmo::BirkhoffLMO, direction, x; kwa
             SparseArrays.sparsevec(linear_indices, V, n^2)
         end
     end
-
     return m
 end
 
@@ -323,6 +326,16 @@ end
 The sum of each row and column has to be equal to 1.
 """
 function Boscia.is_linear_feasible(blmo::BirkhoffLMO, v::AbstractVector)
+    for (i, int_var) in enumerate(blmo.int_vars)
+        if !(
+            blmo.lower_bounds[i] ≤ v[int_var] + 1e-6 || !(v[int_var] - 1e-6 ≤ blmo.upper_bounds[i])
+        )
+            @debug(
+                "Variable: $(int_var) Vertex entry: $(v[int_var]) Lower bound: $(blmo.lower_bounds[i]) Upper bound: $(blmo.upper_bounds[i]))"
+            )
+            return false
+        end
+    end
     n = blmo.dim
     for i in 1:n
         # append by column ? column sum : row sum 
@@ -392,10 +405,10 @@ end
 # Change the value of the bound c_idx.
 function Boscia.set_bound!(blmo::BirkhoffLMO, c_idx, value, sense::Symbol)
     # Reset the lmo if necessary
-    if !blmo.updated_lmo
+    if blmo.updated_lmo
         empty!(blmo.fixed_to_one_rows)
         empty!(blmo.fixed_to_one_cols)
-        blmo.updated_lmo = true
+        blmo.updated_lmo = false
     end
     if sense == :greaterthan
         blmo.lower_bounds[c_idx] = value
@@ -453,7 +466,7 @@ function Boscia.delete_bounds!(blmo::BirkhoffLMO, cons_delete)
     empty!(blmo.index_map_cols)
     append!(blmo.index_map_rows, index_map_rows)
     append!(blmo.index_map_cols, index_map_cols)
-    blmo.updated_lmo = false
+    blmo.updated_lmo = true
     return true
 end
 
@@ -500,6 +513,26 @@ function Boscia.has_integer_constraint(blmo::BirkhoffLMO, idx)
     return idx in blmo.int_vars
 end
 
+## Safety Functions
+
+# Check if the bounds were set correctly in build_LMO.
+# Safety check only.
+function Boscia.build_LMO_correct(blmo::BirkhoffLMO, node_bounds)
+    for key in keys(node_bounds.lower_bounds)
+        idx = findfirst(x -> x == key, blmo.int_vars)
+        if idx === nothing || blmo.lower_bounds[idx] != node_bounds[key, :greaterthan]
+            return false
+        end
+    end
+    for key in keys(node_bounds.upper_bounds)
+        idx = findfirst(x -> x == key, blmo.int_vars)
+        if idx === nothing || blmo.upper_bounds[idx] != node_bounds[key, :lessthan]
+            return false
+        end
+    end
+    return true
+end
+
 ## Optional
 
 function Boscia.check_feasibility(blmo::BirkhoffLMO)
@@ -563,3 +596,4 @@ function Boscia.check_feasibility(blmo::BirkhoffLMO)
 
     return Boscia.OPTIMAL
 end
+

From 591643fb2011b3b75241be6252b62d68e1de6b67 Mon Sep 17 00:00:00 2001
From: Wenjie Xiao <116566662+WenjieXiao-2022@users.noreply.github.com>
Date: Mon, 15 Sep 2025 10:28:09 +0200
Subject: [PATCH 18/21] format

---
 src/birkhoff_polytope.jl | 1 -
 1 file changed, 1 deletion(-)

diff --git a/src/birkhoff_polytope.jl b/src/birkhoff_polytope.jl
index c696db8..7bbe284 100644
--- a/src/birkhoff_polytope.jl
+++ b/src/birkhoff_polytope.jl
@@ -59,7 +59,6 @@ BirkhoffLMO(dim; append_by_column=true, atol=1e-6, rtol=1e-3) = BirkhoffLMO(
 """
 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.compute_extreme_point(lmo::BirkhoffLMO, d; kwargs...)
     n = lmo.dim
 

From b5a8fd3417d417ec4a439198b23081dbd7980ef7 Mon Sep 17 00:00:00 2001
From: Wenjie Xiao <116566662+WenjieXiao-2022@users.noreply.github.com>
Date: Mon, 15 Sep 2025 10:28:45 +0200
Subject: [PATCH 19/21] form

---
 src/birkhoff_polytope.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/birkhoff_polytope.jl b/src/birkhoff_polytope.jl
index 7bbe284..943fd58 100644
--- a/src/birkhoff_polytope.jl
+++ b/src/birkhoff_polytope.jl
@@ -226,8 +226,8 @@ function Boscia.compute_inface_extreme_point(lmo::BirkhoffLMO, direction, x; kwa
                 else
                     d2[j, i] = missing
                 end
-                # the problem can only be integer types,
-                # either full-integer or mixed-integer.
+            # the problem can only be integer types,
+            # either full-integer or mixed-integer.
             elseif length(int_vars) == n^2 || idx_map_ub[orig_linear_idx] != 0
                 idx = length(int_vars) < n^2 ? idx_map_ub[orig_linear_idx] : orig_linear_idx
                 # interdict arc when fixed to zero

From 9c2686023a662873d39d1ddc35c63332c3a573b4 Mon Sep 17 00:00:00 2001
From: WenjieXiao-2022 <shawwj2020@126.com>
Date: Tue, 7 Oct 2025 16:09:04 +0200
Subject: [PATCH 20/21] implemented vector-to-matrix conversion

---
 src/birkhoff_polytope.jl | 51 +++++++++++++++++-----------------------
 1 file changed, 21 insertions(+), 30 deletions(-)

diff --git a/src/birkhoff_polytope.jl b/src/birkhoff_polytope.jl
index 943fd58..a997bf9 100644
--- a/src/birkhoff_polytope.jl
+++ b/src/birkhoff_polytope.jl
@@ -59,13 +59,9 @@ BirkhoffLMO(dim; append_by_column=true, atol=1e-6, rtol=1e-3) = BirkhoffLMO(
 """
 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.compute_extreme_point(lmo::BirkhoffLMO, d; kwargs...)
+function Boscia.compute_extreme_point(lmo::BirkhoffLMO, d::AbstractMatrix{T}; kwargs...) where {T}
     n = lmo.dim
 
-    if size(d, 2) == 1
-        d = lmo.append_by_column ? reshape(d, (n, n)) : transpose(reshape(d, (n, n)))
-    end
-
     fixed_to_one_rows = lmo.fixed_to_one_rows
     fixed_to_one_cols = lmo.fixed_to_one_cols
     index_map_rows = lmo.index_map_rows
@@ -135,8 +131,14 @@ function Boscia.compute_extreme_point(lmo::BirkhoffLMO, d; kwargs...)
         m[index_map_rows[rows[i]], index_map_cols[cols[i]]] = (vals[i] == 2)
     end
 
-    if size(d, 2) == 1 || lmo.boscia_use
-        m = if lmo.append_by_column
+    return m
+end
+
+function Boscia.compute_extreme_point(lmo::BirkhoffLMO, d::AbstractVector{T}; kwargs...) where {T}
+    n = lmo.dim
+    d = lmo.append_by_column ? reshape(d, (n, n)) : transpose(reshape(d, (n, n)))
+    m = Boscia.compute_extreme_point(lmo, d; kwargs...)
+    m = if lmo.append_by_column
             # Convert sparse matrix to sparse vector by columns
             I, J, V = SparseArrays.findnz(m)
             linear_indices = (J .- 1) .* n .+ I
@@ -148,25 +150,14 @@ function Boscia.compute_extreme_point(lmo::BirkhoffLMO, d; kwargs...)
             linear_indices = (J .- 1) .* n .+ I
             SparseArrays.sparsevec(linear_indices, V, n^2)
         end
-    end
     return m
 end
 
-
 """
 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.compute_inface_extreme_point(lmo::BirkhoffLMO, direction, x; kwargs...)
+function Boscia.compute_inface_extreme_point(llmo::BirkhoffLMO, direction::AbstractMatrix{T}, x::AbstractMatrix{T}; kwargs...) where {T}
     n = lmo.dim
-
-    if size(direction, 2) == 1
-        direction =
-            lmo.append_by_column ? reshape(direction, (n, n)) :
-            transpose(reshape(direction, (n, n)))
-
-        x = lmo.append_by_column ? reshape(x, (n, n)) : transpose(reshape(x, (n, n)))
-    end
-
     # Precompute index mapping to avoid repeated `findfirst` calls,
     # which would be very costly inside the loop.
     if length(lmo.int_vars) !== n^2
@@ -265,8 +256,15 @@ function Boscia.compute_inface_extreme_point(lmo::BirkhoffLMO, direction, x; kwa
         m[index_map_rows[rows[i]], index_map_cols[cols[i]]] = (vals[i] == 2)
     end
 
-    if size(direction, 2) == 1 || lmo.boscia_use
-        m = if lmo.append_by_column
+    return m
+end
+
+function Boscia.compute_inface_extreme_point(lmo::BirkhoffLMO, direction::AbstractVector{T}, x::AbstractVector{T}; kwargs...) where {T}
+    n = lmo.dim
+    direction = lmo.append_by_column ? reshape(direction, (n, n)) : transpose(reshape(direction, (n, n)))
+    x = lmo.append_by_column ? reshape(x, (n, n)) : transpose(reshape(x, (n, n)))
+    m = Boscia.compute_inface_extreme_point(lmo, direction, x; kwargs...)
+    m = if lmo.append_by_column
             # Convert sparse matrix to sparse vector by columns
             I, J, V = SparseArrays.findnz(m)
             linear_indices = (J .- 1) .* n .+ I
@@ -278,7 +276,7 @@ function Boscia.compute_inface_extreme_point(lmo::BirkhoffLMO, direction, x; kwa
             linear_indices = (J .- 1) .* n .+ I
             SparseArrays.sparsevec(linear_indices, V, n^2)
         end
-    end
+    
     return m
 end
 
@@ -289,15 +287,8 @@ Fixings are maintained by the oracle (or deduced from `x` itself).
 function Boscia.dicg_maximum_step(lmo::BirkhoffLMO, direction, x; kwargs...)
     n = lmo.dim
     T = promote_type(eltype(x), eltype(direction))
-    if size(direction, 2) == 1
-        direction =
-            lmo.append_by_column ? reshape(direction, (n, n)) :
-            transpose(reshape(direction, (n, n)))
-
-        x = lmo.append_by_column ? reshape(x, (n, n)) : transpose(reshape(x, (n, n)))
-    end
-
     gamma_max = one(T)
+    
     for idx in eachindex(x)
         if direction[idx] != 0.0
             # iterate already on the boundary

From 21f65d77495cc18449b56c048990d356145c629a Mon Sep 17 00:00:00 2001
From: Wenjie Xiao <116566662+WenjieXiao-2022@users.noreply.github.com>
Date: Mon, 20 Oct 2025 19:03:13 +0200
Subject: [PATCH 21/21] removed `type` and renamed DICG functions

---
 src/birkhoff_polytope.jl | 80 ++++++++++++++++++++++------------------
 1 file changed, 44 insertions(+), 36 deletions(-)

diff --git a/src/birkhoff_polytope.jl b/src/birkhoff_polytope.jl
index a997bf9..d92b82b 100644
--- a/src/birkhoff_polytope.jl
+++ b/src/birkhoff_polytope.jl
@@ -80,8 +80,7 @@ function Boscia.compute_extreme_point(lmo::BirkhoffLMO, d::AbstractMatrix{T}; kw
     end
 
     nreduced = length(index_map_rows)
-    type = typeof(d[1, 1])
-    d2 = ones(Union{type,Missing}, nreduced, nreduced)
+    d2 = ones(Union{T,Missing}, nreduced, nreduced)
 
     for j in 1:nreduced
         col_orig = index_map_cols[j]
@@ -139,24 +138,29 @@ function Boscia.compute_extreme_point(lmo::BirkhoffLMO, d::AbstractVector{T}; kw
     d = lmo.append_by_column ? reshape(d, (n, n)) : transpose(reshape(d, (n, n)))
     m = Boscia.compute_extreme_point(lmo, d; kwargs...)
     m = if lmo.append_by_column
-            # Convert sparse matrix to sparse vector by columns
-            I, J, V = SparseArrays.findnz(m)
-            linear_indices = (J .- 1) .* n .+ I
-            SparseArrays.sparsevec(linear_indices, V, n^2)
-        else
-            # Convert sparse matrix to sparse vector by rows (transpose first)
-            mt = SparseArrays.sparse(LinearAlgebra.transpose(m))
-            I, J, V = SparseArrays.findnz(mt)
-            linear_indices = (J .- 1) .* n .+ I
-            SparseArrays.sparsevec(linear_indices, V, n^2)
-        end
+        # Convert sparse matrix to sparse vector by columns
+        I, J, V = SparseArrays.findnz(m)
+        linear_indices = (J .- 1) .* n .+ I
+        SparseArrays.sparsevec(linear_indices, V, n^2)
+    else
+        # Convert sparse matrix to sparse vector by rows (transpose first)
+        mt = SparseArrays.sparse(LinearAlgebra.transpose(m))
+        I, J, V = SparseArrays.findnz(mt)
+        linear_indices = (J .- 1) .* n .+ I
+        SparseArrays.sparsevec(linear_indices, V, n^2)
+    end
     return m
 end
 
 """
 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.compute_inface_extreme_point(llmo::BirkhoffLMO, direction::AbstractMatrix{T}, x::AbstractMatrix{T}; kwargs...) where {T}
+function FrankWolfe.compute_inface_extreme_point(
+    lmo::BirkhoffLMO,
+    direction::AbstractMatrix{T},
+    x::AbstractMatrix{T};
+    kwargs...,
+) where {T}
     n = lmo.dim
     # Precompute index mapping to avoid repeated `findfirst` calls,
     # which would be very costly inside the loop.
@@ -200,8 +204,7 @@ function Boscia.compute_inface_extreme_point(llmo::BirkhoffLMO, direction::Abstr
     deleteat!(index_map_cols, delete_index_map_cols)
 
     nreduced = length(index_map_rows)
-    type = typeof(direction[1, 1])
-    d2 = ones(Union{type,Missing}, nreduced, nreduced)
+    d2 = ones(Union{T,Missing}, nreduced, nreduced)
     for j in 1:nreduced
         col_orig = index_map_cols[j]
         for i in 1:nreduced
@@ -217,8 +220,8 @@ function Boscia.compute_inface_extreme_point(llmo::BirkhoffLMO, direction::Abstr
                 else
                     d2[j, i] = missing
                 end
-            # the problem can only be integer types,
-            # either full-integer or mixed-integer.
+                # the problem can only be integer types,
+                # either full-integer or mixed-integer.
             elseif length(int_vars) == n^2 || idx_map_ub[orig_linear_idx] != 0
                 idx = length(int_vars) < n^2 ? idx_map_ub[orig_linear_idx] : orig_linear_idx
                 # interdict arc when fixed to zero
@@ -259,24 +262,30 @@ function Boscia.compute_inface_extreme_point(llmo::BirkhoffLMO, direction::Abstr
     return m
 end
 
-function Boscia.compute_inface_extreme_point(lmo::BirkhoffLMO, direction::AbstractVector{T}, x::AbstractVector{T}; kwargs...) where {T}
+function FrankWolfe.compute_inface_extreme_point(
+    lmo::BirkhoffLMO,
+    direction::AbstractVector{T},
+    x::AbstractVector{T};
+    kwargs...,
+) where {T}
     n = lmo.dim
-    direction = lmo.append_by_column ? reshape(direction, (n, n)) : transpose(reshape(direction, (n, n)))
+    direction =
+        lmo.append_by_column ? reshape(direction, (n, n)) : transpose(reshape(direction, (n, n)))
     x = lmo.append_by_column ? reshape(x, (n, n)) : transpose(reshape(x, (n, n)))
     m = Boscia.compute_inface_extreme_point(lmo, direction, x; kwargs...)
     m = if lmo.append_by_column
-            # Convert sparse matrix to sparse vector by columns
-            I, J, V = SparseArrays.findnz(m)
-            linear_indices = (J .- 1) .* n .+ I
-            SparseArrays.sparsevec(linear_indices, V, n^2)
-        else
-            # Convert sparse matrix to sparse vector by rows (transpose first)
-            mt = SparseArrays.sparse(LinearAlgebra.transpose(m))
-            I, J, V = SparseArrays.findnz(mt)
-            linear_indices = (J .- 1) .* n .+ I
-            SparseArrays.sparsevec(linear_indices, V, n^2)
-        end
-    
+        # Convert sparse matrix to sparse vector by columns
+        I, J, V = SparseArrays.findnz(m)
+        linear_indices = (J .- 1) .* n .+ I
+        SparseArrays.sparsevec(linear_indices, V, n^2)
+    else
+        # Convert sparse matrix to sparse vector by rows (transpose first)
+        mt = SparseArrays.sparse(LinearAlgebra.transpose(m))
+        I, J, V = SparseArrays.findnz(mt)
+        linear_indices = (J .- 1) .* n .+ I
+        SparseArrays.sparsevec(linear_indices, V, n^2)
+    end
+
     return m
 end
 
@@ -284,11 +293,11 @@ 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.dicg_maximum_step(lmo::BirkhoffLMO, direction, x; kwargs...)
+function FrankWolfe.dicg_maximum_step(lmo::BirkhoffLMO, direction, x; kwargs...)
     n = lmo.dim
     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
@@ -308,7 +317,7 @@ function Boscia.dicg_maximum_step(lmo::BirkhoffLMO, direction, x; kwargs...)
 
 end
 
-function Boscia.is_decomposition_invariant_oracle(lmo::BirkhoffLMO)
+function FrankWolfe.is_decomposition_invariant_oracle(lmo::BirkhoffLMO)
     return true
 end
 
@@ -586,4 +595,3 @@ function Boscia.check_feasibility(blmo::BirkhoffLMO)
 
     return Boscia.OPTIMAL
 end
-