From 6090be27234ce0baac60f61e7b33938648799412 Mon Sep 17 00:00:00 2001
From: dhendryc <92737336+dhendryc@users.noreply.github.com>
Date: Fri, 20 Sep 2024 13:41:52 +0200
Subject: [PATCH] Lower bound improvement based on sharpness of the objective
 (#195)
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

* Move tightening functions into a separate file.

* Tighten bounds using sharpness.

* Add sharpness constant and exponent to the interface.

* Minor fixes.

* Tests for strong convexity and sharpness on the Normbox example.

* Use the node dual gap instead of root dual gap.

* Make sure that we collect the statistics, even if in case of timeout.

* We only need this if the key is not in the result dictionary.

* Move the BnB Callback into the appriopiate file.

* We don't need the root dual gap.

* Update src/custom_bonobo.jl

Co-authored-by: Mathieu Besançon <mathieu.besancon@gmail.com>

* Add a bit of give to the sharpness constant to avoid cutting of optimal solutions.

* Added small definition of sharpness.

* Update src/callbacks.jl

* Update src/callbacks.jl

* Don't create child node if the domain oracle is infeasible.

* Set version up after bug fix.

* Optimal Experiment Design example.

* Fix syntax issues.

* Apply suggestions from code review

* Add corrected sharpness constants.

* In case the vertex is domain infeasible, use x to get the data type of f(x).

* Syntax issue fix.

* Corrections.

* Print statement in case assert fails.

* Print logs for better bug hunting.

* Another try to reproduce error.

* Test lower bound against tree.incumbent plus current fw_epsilon.

* Have the same seed as in runtests.jl.

* Rename A to Ex_mat to avoid constant redefinition.

* Clean up another constant redefinition.

* Add test wrapper.

* Debug statements.

* More debug statements.

* If heuristic return infeasible solution, ignore it.

* Minor change.

* Clean up example.

* Show solutions.

* Don't use vertex storage for the heuristics.

* Final clean up.

* Run tests with one specific seed.

* Added time limit.

* Corrected sharpness constant.

* Corrected sharpness constant and slight relaxation of test for A-Optimal.

* Update src/tightenings.jl

---------

Co-authored-by: Hendrych <dhendryc@mob-009366.zib.de>
Co-authored-by: Mathieu Besançon <mathieu.besancon@gmail.com>
---
 examples/birkhoff.jl                  |   6 -
 examples/int_sparse_reg.jl            |   6 -
 examples/mps-example.jl               |   4 -
 examples/mps-examples/mip-examples.jl |   4 -
 examples/nonlinear.jl                 |   6 +-
 examples/oed_utils.jl                 | 293 +++++++++++++++++++++++++
 examples/optimal_experiment_design.jl | 129 +++++++++++
 examples/poisson_reg.jl               |   6 -
 examples/portfolio.jl                 |   4 -
 src/Boscia.jl                         |   1 +
 src/MOI_bounded_oracle.jl             |   4 +
 src/callbacks.jl                      | 259 +++++++++++++++++++++-
 src/custom_bonobo.jl                  |  19 ++
 src/heuristics.jl                     |   7 +-
 src/interface.jl                      | 267 ++---------------------
 src/managed_blmo.jl                   |   4 +
 src/node.jl                           | 253 +---------------------
 src/tightenings.jl                    | 300 ++++++++++++++++++++++++++
 test/interface_test.jl                |  57 +++++
 19 files changed, 1085 insertions(+), 544 deletions(-)
 create mode 100644 examples/oed_utils.jl
 create mode 100644 examples/optimal_experiment_design.jl
 create mode 100644 src/tightenings.jl

diff --git a/examples/birkhoff.jl b/examples/birkhoff.jl
index 7c7025b48..655ef7ccf 100644
--- a/examples/birkhoff.jl
+++ b/examples/birkhoff.jl
@@ -12,12 +12,6 @@ import HiGHS
 # https://arxiv.org/pdf/2011.02752.pdf
 # https://www.sciencedirect.com/science/article/pii/S0024379516001257
 
-# For bug hunting:
-seed = rand(UInt64)
-@show seed
-seed = 0x3eb09305cecf69f0
-Random.seed!(seed)
-
 
 # min_{X, θ} 1/2 * || ∑_{i in [k]} θ_i X_i - Xhat ||^2
 # θ ∈ Δ_k (simplex)
diff --git a/examples/int_sparse_reg.jl b/examples/int_sparse_reg.jl
index 1d7a0a3b9..9d6a037d8 100644
--- a/examples/int_sparse_reg.jl
+++ b/examples/int_sparse_reg.jl
@@ -21,12 +21,6 @@ const MOI = MathOptInterface
 # r - controls how often we have to maximal split on a index.
 # k - is the sparsity parameter. We only want a few non zero entries.
 
-# For bug hunting:
-seed = rand(UInt64)
-@show seed
-#seed = 0xeadb922ca734998b  
-Random.seed!(seed)
-
 n = 10
 m = 30
 l = 5
diff --git a/examples/mps-example.jl b/examples/mps-example.jl
index 6ea84f0ed..e60208919 100644
--- a/examples/mps-example.jl
+++ b/examples/mps-example.jl
@@ -16,10 +16,6 @@ const MOI = MathOptInterface
 # Number of binaries    231
 # Number of constraints 198
 
-seed = rand(UInt64)
-@show seed
-Random.seed!(seed)
-
 src = MOI.FileFormats.Model(filename="22433.mps")
 MOI.read_from_file(src, joinpath(@__DIR__, "mps-examples/mps-files/22433.mps"))
 
diff --git a/examples/mps-examples/mip-examples.jl b/examples/mps-examples/mip-examples.jl
index 173dfd089..aba7854be 100644
--- a/examples/mps-examples/mip-examples.jl
+++ b/examples/mps-examples/mip-examples.jl
@@ -21,10 +21,6 @@ import Ipopt
 # ran14x18-disj-8     https://miplib.zib.de/instance_details_ran14x18-disj-8.html
 # timtab1             https://miplib.zib.de/instance_details_timtab1.html   (Takes LONG!)
 
-seed = rand(UInt64)
-@show seed
-Random.seed!(seed)
-
 # To see debug statements
 #ENV["JULIA_DEBUG"] = "Boscia"
 
diff --git a/examples/nonlinear.jl b/examples/nonlinear.jl
index 4fac4e8f9..fdc27539d 100644
--- a/examples/nonlinear.jl
+++ b/examples/nonlinear.jl
@@ -11,10 +11,6 @@ using Dates
 # using SCIP
 # const MOI = MathOptInterface
 
-n = 40
-seed = 10
-
-Random.seed!(seed)
 
 ################################################################
 # alternative implementation of LMO using MOI and SCIP
@@ -113,6 +109,6 @@ heu2 = Boscia.Heuristic(Boscia.rounding_lmo_01_heuristic, 0.8, :lmo_rounding)
 heuristics = [heu, heu2]
 # heuristics = []
 
-x, _, _ = Boscia.solve(f, grad!, lmo, verbose=true, print_iter=500, custom_heuristics=heuristics)
+x, _, _ = Boscia.solve(f, grad!, lmo, verbose=true, print_iter=500, custom_heuristics=heuristics, time_limit=300)
 
 @show x
diff --git a/examples/oed_utils.jl b/examples/oed_utils.jl
new file mode 100644
index 000000000..32064d071
--- /dev/null
+++ b/examples/oed_utils.jl
@@ -0,0 +1,293 @@
+"""
+    build_data
+
+seed - for the Random functions.
+m    - number of experiments.
+fusion - boolean deiciding whether we build the fusion or standard problem.
+corr - boolean deciding whether we build the independent or correlated data.   
+"""
+function build_data(m)
+    # set up
+    n = Int(floor(m/10))
+    
+    B = rand(m,n)
+    B = B'*B
+    @assert isposdef(B)
+    D = MvNormal(randn(n),B)
+    
+    A = rand(D, m)'
+    @assert rank(A) == n 
+
+    N = floor(1.5*n)
+    u = floor(N/3)
+    ub = rand(1.0:u, m)
+        
+    return A, n, N, ub 
+end
+
+"""
+Build Probability Simplex BLMO for Boscia
+"""
+function build_blmo(m, N, ub)
+    simplex_lmo = Boscia.ProbabilitySimplexSimpleBLMO(N)
+    blmo = Boscia.ManagedBoundedLMO(simplex_lmo, fill(0.0, m), ub, collect(1:m), m)
+    return blmo
+end
+
+"""
+Build function and the gradient for the A-criterion. 
+There is the option to build a safe version of the objective and gradient. If x does not
+satisfy the domain oracle, infinity is returned. 
+If one builds the unsafe version, a FrankWolfe line search must be chosen that can take 
+a domain oracle as an input like Secant or MonotonicGenericStepSize. 
+"""
+function build_a_criterion(A; μ=1e-4, build_safe=false)
+    m, n = size(A) 
+    a = m
+    domain_oracle = build_domain_oracle(A, n)
+
+    function f_a(x)
+        X = transpose(A)*diagm(x)*A + Matrix(μ *I, n, n)
+        X = Symmetric(X)
+        U = cholesky(X)
+        X_inv = U \ I
+        return LinearAlgebra.tr(X_inv)/a 
+    end
+
+    function grad_a!(storage, x)
+        X = transpose(A)*diagm(x)*A + Matrix(μ *I, n, n)
+        X = Symmetric(X*X)
+        F = cholesky(X)
+        for i in 1:length(x)
+            storage[i] = LinearAlgebra.tr(- (F \ A[i,:]) * transpose(A[i,:]))/a
+        end
+        return storage 
+    end 
+
+    function f_a_safe(x)
+        if !domain_oracle(x)
+            return Inf
+        end
+        X = transpose(A)*diagm(x)*A + Matrix(μ *I, n, n)
+        X = Symmetric(X)
+        X_inv = LinearAlgebra.inv(X)
+        return LinearAlgebra.tr(X_inv)/a 
+    end
+
+    function grad_a_safe!(storage, x)
+        if !domain_oracle(x)
+            return fill(Inf, length(x))        
+        end
+        X = transpose(A)*diagm(x)*A + Matrix(μ *I, n, n)
+        X = Symmetric(X*X)
+        F = cholesky(X)
+        for i in 1:length(x)
+            storage[i] = LinearAlgebra.tr(- (F \ A[i,:]) * transpose(A[i,:]))/a
+        end
+        return storage
+    end
+
+    if build_safe
+        return f_a_safe, grad_a_safe!
+    end
+
+    return f_a, grad_a!
+end
+
+"""
+Build function and gradient for the D-criterion.
+There is the option to build a safe version of the objective and gradient. If x does not
+satisfy the domain oracle, infinity is returned. 
+If one builds the unsafe version, a FrankWolfe line search must be chosen that can take 
+a domain oracle as an input like Secant or MonotonicGenericStepSize. 
+"""
+function build_d_criterion(A; μ =0.0, build_safe=false)
+    m, n = size(A)
+    a = 1#m
+    domain_oracle = build_domain_oracle(A, n)
+
+    function f_d(x)
+        X = transpose(A)*diagm(x)*A + Matrix(μ *I, n, n)
+        X = Symmetric(X)
+        return float(-log(det(X))/a)
+    end
+
+    function grad_d!(storage, x)
+        X = transpose(A)*diagm(x)*A + Matrix(μ *I, n, n)
+        X = Symmetric(X)
+        F = cholesky(X) 
+        for i in 1:length(x)        
+            storage[i] = 1/a * LinearAlgebra.tr(-(F \ A[i,:] )*transpose(A[i,:])) 
+        end
+        return storage
+    end
+
+    function f_d_safe(x)
+        if !domain_oracle(x)
+            return Inf
+        end
+        X = transpose(A)*diagm(x)*A + Matrix(μ *I, n, n)
+        X = Symmetric(X)
+        return float(-log(det(X))/a)
+    end
+
+    function grad_d_safe!(storage, x)
+        if !domain_oracle(x)
+            return fill(Inf, length(x))
+        end
+        X = transpose(A)*diagm(x)*A + Matrix(μ *I, n, n)
+        X = Symmetric(X)
+        F = cholesky(X) 
+        for i in 1:length(x)        
+            storage[i] = 1/a * LinearAlgebra.tr(-(F \ A[i,:] )*transpose(A[i,:])) 
+        end
+        # https://stackoverflow.com/questions/46417005/exclude-elements-of-array-based-on-index-julia
+        return storage
+    end
+
+    if build_safe
+        return f_d_safe, grad_d_safe!
+    end
+
+    return f_d, grad_d!
+end
+
+"""
+Find n linearly independent rows of A to build the starting point.
+"""
+function linearly_independent_rows(A)
+    S = []
+    m, n = size(A)
+    for i in 1:m
+        S_i= vcat(S, i)
+        if rank(A[S_i,:])==length(S_i)
+            S=S_i
+        end
+        if length(S) == n # we only n linearly independent points
+            return S
+        end
+    end 
+    return S # then x= zeros(m) and x[S] = 1
+end
+
+function add_to_min(x, u)
+    perm = sortperm(x)
+    j = findfirst(x->x != 0, x[perm])
+    
+    for i in j:length(x)
+        if x[perm[i]] < u[perm[i]]
+            x[perm[i]] += 1
+            break
+        else
+            continue
+        end
+    end
+    return x
+end
+
+function remove_from_max(x)
+    perm = sortperm(x, rev = true)
+    j = findlast(x->x != 0, x[perm])
+    
+    for i in 1:j
+        if x[perm[i]] > 1
+            x[perm[i]] -= 1
+            break
+        else
+            continue
+        end
+    end
+    return x
+end
+
+"""
+Build start point used in FrankWolfe and Boscia for the A-Optimal and D-Optimal Design Problem.
+The functions are self concordant and so not every point in the feasible region
+is in the domain of f and grad!.
+""" 
+function build_start_point(A, N, ub)
+    # Get n linearly independent rows of A
+    m, n = size(A)
+    S = linearly_independent_rows(A)
+    @assert length(S) == n
+    
+    x = zeros(m)
+    E = []
+    V = Vector{Float64}[]
+
+    while !isempty(setdiff(S, findall(x-> !(iszero(x)),x)))
+        v = zeros(m)
+        while sum(v) < N
+            idx = isempty(setdiff(S, findall(x-> !(iszero(x)),v))) ? rand(setdiff(collect(1:m), S)) : rand(setdiff(S, findall(x-> !(iszero(x)),v)))
+            if !isapprox(v[idx], 0.0)
+                @debug "Index $(idx) already picked"
+                continue
+            end
+            v[idx] = min(ub[idx], N - sum(v))
+            push!(E, idx)
+        end
+        push!(V,v)
+        x = sum(V .* 1/length(V)) 
+    end
+    unique!(V)
+    a = length(V)
+    x = sum(V .* 1/a)
+    active_set= FrankWolfe.ActiveSet(fill(1/a, a), V, x)
+
+    return x, active_set, S
+end
+
+"""
+Create first incumbent for Boscia and custom BB in a greedy fashion.
+"""
+function greedy_incumbent(A, N, ub)
+    # Get n linearly independent rows of A
+    m, n = size(A)
+    S = linearly_independent_rows(A)
+    @assert length(S) == n
+
+    # set entries to their upper bound
+    x = zeros(m)
+    x[S] .= ub[S]
+
+    if isapprox(sum(x), N; atol=1e-4, rtol=1e-2)
+        return x
+    elseif sum(x) > N
+        while sum(x) > N
+            remove_from_max(x)
+        end
+    elseif sum(x) < N
+        S1 = S
+        while sum(x) < N
+            jdx = rand(setdiff(collect(1:m), S1))
+            x[jdx] = min(N-sum(x), ub[jdx])
+            push!(S1,jdx)
+            sort!(S1)
+        end
+    end
+    @assert isapprox(sum(x), N; atol=1e-4, rtol=1e-2)
+    @assert sum(ub - x .>= 0) == m 
+    return x
+end
+
+"""
+Check if given point is in the domain of f, i.e. X = transpose(A) * diagm(x) * A 
+positive definite.
+
+(a) Check the rank of A restricted to the rows activated by x.
+(b) Check if the resulting information matrix A' * diagm(x) * A is psd.
+
+(b) is a bit faster for smaller dimensions (< 100). For larger (> 200) (a) is faster.
+"""
+function build_domain_oracle(A, n)
+    return function domain_oracle(x)
+        S = findall(x-> !iszero(x),x)
+        return length(S) >= n && rank(A[S,:]) == n 
+    end
+end
+
+function build_domain_oracle2(A)
+    return function domain_oracle2(x)
+        return isposdef(Symmetric(A' * diagm(x) * A))
+    end 
+end
diff --git a/examples/optimal_experiment_design.jl b/examples/optimal_experiment_design.jl
new file mode 100644
index 000000000..236c40939
--- /dev/null
+++ b/examples/optimal_experiment_design.jl
@@ -0,0 +1,129 @@
+using Boscia
+using Random
+using Distributions
+using LinearAlgebra
+using FrankWolfe
+using Statistics
+using Test 
+
+# The function building the problem data and other structures is in a separate file.
+include("oed_utils.jl")
+
+"""
+    The Optimal Experiment Design Problem consists of choosing a subset of experiments
+    maximising the information gain. 
+    We generate the Experiment matrix A ∈ R^{mxn} randomly. The information matrix is a linear
+    map X(x) = A' * diag(x) * A. There exists different information measures Φ. 
+    We concentrate on the A-criterion and D-criterion, i.e. 
+
+        Trace(X^{-1})       (A-Criterion)
+    and
+        -logdet(X)          (D-Criterion).
+
+    Consequently, the optimization problems we want to solve are
+
+    min_x  Trace( (A' * diag(x) * A)^{-1} )
+    s.t.   ∑ x_i = N                            (A-Optimal Design Problem)
+           0 ≤ x ≤ u
+
+           min_x  -logdet(A' * diag(x) * A)
+    s.t.   ∑ x_i = N                            (D-Optimal Design Problem)
+           0 ≤ x ≤ u
+
+    where N is our bugdet for the experiments, i.e. this is the amount of experiments
+    we can perform. We set N = 3/2 * n. The upperbounds u are randomly generated. 
+    
+    Also, check this paper: https://arxiv.org/abs/2312.11200 and the corresponding 
+    respository https://github.com/ZIB-IOL/OptimalDesignWithBoscia.
+
+    A continuous version of the problem can be found in the examples in FrankWolfe.jl:
+    https://github.com/ZIB-IOL/FrankWolfe.jl/blob/master/examples/optimal_experiment_design.jl
+"""
+seed = rand(UInt64)
+
+@show seed
+Random.seed!(seed)
+m = 40
+verbose = true
+
+## A-Optimal Design Problem
+
+@testset "A-Optimal Design" begin
+
+    Ex_mat, n, N, ub = build_data(m)
+
+    @show Ex_mat
+
+    # sharpness constants
+    σ = minimum(Ex_mat' * Ex_mat)
+    λ_max = maximum(ub) * maximum([norm(Ex_mat[i,:])^2 for i=1:size(Ex_mat,1)])
+    θ = 1/2
+    M = sqrt(λ_max^3/  n * σ^4)
+
+
+    g, grad! = build_a_criterion(Ex_mat, build_safe=true)
+    blmo = build_blmo(m, N, ub)
+    x0, active_set = build_start_point(Ex_mat, N, ub)
+    z = greedy_incumbent(Ex_mat, N, ub)
+    domain_oracle = build_domain_oracle(Ex_mat, n)
+    line_search = FrankWolfe.MonotonicGenericStepsize(FrankWolfe.Adaptive(), domain_oracle)
+    heu = Boscia.Heuristic(Boscia.rounding_hyperplane_heuristic, 0.7, :hyperplane_aware_rounding)
+
+    x, _, result = Boscia.solve(g, grad!, blmo, active_set=active_set, start_solution=z, verbose=verbose, domain_oracle=domain_oracle, custom_heuristics=[heu], sharpness_exponent=θ, sharpness_constant=M, line_search=line_search) #sharpness_exponent=θ, sharpness_constant=M,
+
+    gradient = similar(x)
+    grad!(gradient, x)
+    v = Boscia.compute_extreme_point(blmo, gradient)
+    dual_gap = FrankWolfe.dot(gradient, x) - FrankWolfe.dot(gradient, v)
+    @show dual_gap
+
+    blmo = build_blmo(m, N, ub)
+    x0, active_set = build_start_point(Ex_mat, N, ub)
+    z = greedy_incumbent(Ex_mat, N, ub)
+    domain_oracle = build_domain_oracle(Ex_mat, n)
+    heu = Boscia.Heuristic(Boscia.rounding_hyperplane_heuristic, 0.7, :hyperplane_aware_rounding)
+
+    x_s, _, result = Boscia.solve(g, grad!, blmo, active_set=active_set, start_solution=z, verbose=verbose, line_search=FrankWolfe.Secant(40, 1e-8, domain_oracle), domain_oracle=domain_oracle, sharpness_exponent=θ, sharpness_constant=M, custom_heuristics=[heu]) #sharpness_exponent=θ, sharpness_constant=M, 
+
+    @show x
+    @show x_s
+    #@show f(x), f(x_s)
+    @show g(x), g(x_s)
+    @test isapprox(g(x), g(x_s), atol=1e-3, rtol=5e-2)
+end 
+
+## D-Optimal Design Problem
+
+@testset "D-optimal Design" begin
+    Ex_mat, n, N, ub = build_data(m)
+
+    # sharpness constants
+    σ = minimum(Ex_mat' * Ex_mat)
+    λ_max = maximum(ub) * maximum([norm(Ex_mat[i,:])^2 for i=1:size(Ex_mat,1)])
+    θ = 1/2
+    M = sqrt(2 * λ_max^2/ n * σ^4 )
+
+
+    g, grad! = build_d_criterion(Ex_mat, build_safe=true)
+    blmo = build_blmo(m, N, ub)
+    x0, active_set = build_start_point(Ex_mat, N, ub)
+    z = greedy_incumbent(Ex_mat, N, ub)
+    domain_oracle = build_domain_oracle(Ex_mat, n)
+    heu = Boscia.Heuristic(Boscia.rounding_hyperplane_heuristic, 0.7, :hyperplane_aware_rounding)
+
+    x, _, result = Boscia.solve(g, grad!, blmo, active_set=active_set, start_solution=z, verbose=verbose,  domain_oracle=domain_oracle, sharpness_exponent=θ, sharpness_constant=M, custom_heuristics=[heu]) 
+
+    blmo = build_blmo(m, N, ub)
+    x0, active_set = build_start_point(Ex_mat, N, ub)
+    z = greedy_incumbent(Ex_mat, N, ub)
+    domain_oracle = build_domain_oracle(Ex_mat, n)
+    heu = Boscia.Heuristic(Boscia.rounding_hyperplane_heuristic, 0.7, :hyperplane_aware_rounding)
+
+    x_s, _, result = Boscia.solve(g, grad!, blmo, active_set=active_set, start_solution=z, verbose=verbose, line_search=FrankWolfe.Secant(40, 1e-8, domain_oracle), domain_oracle=domain_oracle, sharpness_exponent=θ, sharpness_constant=M, custom_heuristics=[heu]) 
+
+    @show x 
+    @show x_s
+    @show g(x), g(x_s)
+    @test isapprox(g(x), g(x_s), atol=1e-4, rtol=1e-2)
+end
+
diff --git a/examples/poisson_reg.jl b/examples/poisson_reg.jl
index 1f87527f1..8bde3f90e 100644
--- a/examples/poisson_reg.jl
+++ b/examples/poisson_reg.jl
@@ -11,12 +11,6 @@ const MOI = MathOptInterface
 
 # Poisson sparse regression
 
-# For bug hunting:
-seed = rand(UInt64)
-@show seed
-#seed = 0xfe03ee83ca373eab   
-Random.seed!(seed)
-
 # min_{w, b, z} ∑_i exp(w x_i + b) - y_i (w x_i + b) + α norm(w)^2
 # s.t. -N z_i <= w_i <= N z_i
 # b ∈ [-N, N]
diff --git a/examples/portfolio.jl b/examples/portfolio.jl
index 8b0221384..22b4c3c22 100644
--- a/examples/portfolio.jl
+++ b/examples/portfolio.jl
@@ -7,10 +7,6 @@ using LinearAlgebra
 import MathOptInterface
 const MOI = MathOptInterface
 
-# seed = 0x946d4b7835e92ffa takes 90 minutes to solve! -> not anymore
-seed = 0x946d4b7835e92ffa
-Random.seed!(seed)
-
 n = 30
 const ri = rand(n)
 const ai = rand(n)
diff --git a/src/Boscia.jl b/src/Boscia.jl
index 5e92fc6ca..6853fa541 100644
--- a/src/Boscia.jl
+++ b/src/Boscia.jl
@@ -17,6 +17,7 @@ include("blmo_interface.jl")
 include("time_tracking_lmo.jl")
 include("frank_wolfe_variants.jl")
 include("build_lmo.jl")
+include("tightenings.jl")
 include("node.jl")
 include("custom_bonobo.jl")
 include("callbacks.jl")
diff --git a/src/MOI_bounded_oracle.jl b/src/MOI_bounded_oracle.jl
index 6039cee10..39e883401 100644
--- a/src/MOI_bounded_oracle.jl
+++ b/src/MOI_bounded_oracle.jl
@@ -629,6 +629,8 @@ function solve(
     global_dual_tightening=true,
     bnb_callback=nothing,
     strong_convexity=0.0,
+    sharpness_constant = 0.0,
+    sharpness_exponent = Inf,
     domain_oracle=x -> true,
     start_solution=nothing,
     fw_verbose=false,
@@ -666,6 +668,8 @@ function solve(
         global_dual_tightening=global_dual_tightening,
         bnb_callback=bnb_callback,
         strong_convexity=strong_convexity,
+        sharpness_constant=sharpness_constant,
+        sharpness_exponent=sharpness_exponent,
         domain_oracle=domain_oracle,
         start_solution=start_solution,
         fw_verbose=fw_verbose,
diff --git a/src/callbacks.jl b/src/callbacks.jl
index a1968d028..67d9062b8 100644
--- a/src/callbacks.jl
+++ b/src/callbacks.jl
@@ -1,4 +1,11 @@
-# FW callback
+"""
+Frank-Wolfe Callback.
+
+Is called in every Frank-Wolfe iteration.
+Node evaluation can be dynamically stopped here.
+Time limit is checked.
+If the vertex is providing a better incumbent, it is added as solution.
+"""
 function build_FW_callback(
     tree,
     min_number_lower,
@@ -84,3 +91,253 @@ function build_FW_callback(
         return true
     end
 end
+
+"""
+Branch-and-Bound Callback.
+Collects statistics and prints logs if verbose is turned on.
+
+Output of Boscia:
+    iter :          current iteration of Boscia
+    node id :       current node id
+    lower bound :   tree_lb(tree)
+    incumbent :     tree.incumbent
+    gap :           tree.incumbent-tree_lb(tree)
+    rel. gap :      dual_gap/tree.incumbent
+    time :          total time of Boscia
+    time/nodes :    average time per node
+    FW time :       time spent in FW 
+    LMO time :      time used by LMO
+    LMO calls :     number of compute_extreme_point calls in FW
+    FW iterations : number of iterations in FW
+"""
+function build_bnb_callback(
+    tree,
+    time_ref,
+    list_lb_cb,
+    list_ub_cb,
+    list_time_cb,
+    list_num_nodes_cb,
+    list_lmo_calls_cb,
+    verbose,
+    fw_iterations,
+    list_active_set_size_cb,
+    list_discarded_set_size_cb,
+    result,
+    lmo_calls_per_layer,
+    active_set_size_per_layer,
+    discarded_set_size_per_layer,
+    node_level,
+    baseline_callback,
+    local_tightenings,
+    global_tightenings,
+    local_potential_tightenings,
+    num_bin,
+    num_int,
+)
+    iteration = 0
+
+    headers = [
+        " ",
+        "Iter",
+        "Open",
+        "Bound",
+        "Incumbent",
+        "Gap (abs)",
+        "Gap (rel)",
+        "Time (s)",
+        "Nodes/sec",
+        "FW (ms)",
+        "LMO (ms)",
+        "LMO (calls c)",
+        "FW (its)",
+        "#activeset",
+        "#shadow",
+    ]
+    format_string = "%1s %5i %5i %14e %14e %14e %14e %14e %14e %12i %10i %14i %10i %8i %8i\n"
+    print_iter = get(tree.root.options, :print_iter, 100)
+
+    if verbose
+        FrankWolfe.print_callback(headers, format_string, print_header=true)
+    end
+    return function callback(
+        tree,
+        node;
+        worse_than_incumbent=false,
+        node_infeasible=false,
+        lb_update=false,
+    )
+        if baseline_callback !== nothing
+            baseline_callback(
+                tree,
+                node,
+                worse_than_incumbent=worse_than_incumbent,
+                node_infeasible=node_infeasible,
+                lb_update=lb_update,
+            )
+        end
+        if !node_infeasible
+            #update lower bound
+            if lb_update == true
+                tree.node_queue[node.id] = (node.lb, node.id)
+                _, prio = peek(tree.node_queue)
+                @assert tree.lb <= prio[1]
+                tree.lb = min(minimum([prio[2][1] for prio in tree.node_queue]), tree.incumbent)
+            end
+            push!(list_ub_cb, tree.incumbent)
+            push!(list_num_nodes_cb, tree.num_nodes)
+            push!(node_level, node.level)
+            iteration += 1
+            if tree.lb == -Inf && isempty(tree.nodes)
+                tree.lb = node.lb
+            end
+
+            time = float(Dates.value(Dates.now() - time_ref))
+            push!(list_time_cb, time)
+
+            if tree.root.options[:time_limit] < Inf
+                if time / 1000.0 ≥ tree.root.options[:time_limit]
+                    @assert tree.root.problem.solving_stage == SOLVING
+                    tree.root.problem.solving_stage = TIME_LIMIT_REACHED
+                end
+            end
+
+            fw_time = Dates.value(node.fw_time)
+            fw_iter = if !isempty(fw_iterations)
+                fw_iterations[end]
+            else
+                0
+            end
+            if !isempty(tree.root.problem.tlmo.optimizing_times)
+                LMO_time = sum(1000 * tree.root.problem.tlmo.optimizing_times)
+                empty!(tree.root.problem.tlmo.optimizing_times)
+            else
+                LMO_time = 0
+            end
+            LMO_calls_c = tree.root.problem.tlmo.ncalls
+            push!(list_lmo_calls_cb, copy(LMO_calls_c))
+
+            if !isempty(tree.node_queue)
+                p_lb = tree.lb
+                tree.lb = min(minimum([prio[2][1] for prio in tree.node_queue]), tree.incumbent)
+                @assert p_lb <= tree.lb + tree.root.options[:dual_gap] "p_lb <= tree.lb + tree.root.options[:dual_gap] $(p_lb) <= $(tree.lb + tree.root.options[:dual_gap])"
+            end
+            # correct lower bound if necessary
+            tree.lb = tree_lb(tree)
+            dual_gap = tree.incumbent - tree_lb(tree)
+            push!(list_lb_cb, tree_lb(tree))
+            active_set_size = length(node.active_set)
+            discarded_set_size = length(node.discarded_vertices.storage)
+            push!(list_active_set_size_cb, active_set_size)
+            push!(list_discarded_set_size_cb, discarded_set_size)
+            nodes_left = length(tree.nodes)
+            if tree.root.updated_incumbent[]
+                incumbent_updated = "*"
+            else
+                incumbent_updated = " "
+            end
+            if verbose && (
+                mod(iteration, print_iter) == 0 ||
+                iteration == 1 ||
+                Bonobo.terminated(tree) ||
+                tree.root.updated_incumbent[]
+            )
+                if (mod(iteration, print_iter * 40) == 0)
+                    FrankWolfe.print_callback(headers, format_string, print_header=true)
+                end
+                FrankWolfe.print_callback(
+                    (
+                        incumbent_updated,
+                        iteration,
+                        nodes_left,
+                        tree_lb(tree),
+                        tree.incumbent,
+                        dual_gap,
+                        relative_gap(tree.incumbent, tree_lb(tree)),
+                        time / 1000.0,
+                        tree.num_nodes / time * 1000.0,
+                        fw_time,
+                        LMO_time,
+                        tree.root.problem.tlmo.ncalls,
+                        fw_iter,
+                        active_set_size,
+                        discarded_set_size,
+                    ),
+                    format_string,
+                    print_header=false,
+                )
+                tree.root.updated_incumbent[] = false
+            end
+            # lmo calls per layer
+            if length(list_lmo_calls_cb) > 1
+                LMO_calls = list_lmo_calls_cb[end] - list_lmo_calls_cb[end-1]
+            else
+                LMO_calls = list_lmo_calls_cb[end]
+            end
+            if length(lmo_calls_per_layer) < node.level
+                push!(lmo_calls_per_layer, [LMO_calls])
+                push!(active_set_size_per_layer, [active_set_size])
+                push!(discarded_set_size_per_layer, [discarded_set_size])
+            else
+                push!(lmo_calls_per_layer[node.level], LMO_calls)
+                push!(active_set_size_per_layer[node.level], active_set_size)
+                push!(discarded_set_size_per_layer[node.level], discarded_set_size)
+            end
+
+            # add tightenings
+            push!(global_tightenings, node.global_tightenings)
+            push!(local_tightenings, node.local_tightenings)
+            push!(local_potential_tightenings, node.local_potential_tightenings)
+            @assert node.local_potential_tightenings <= num_bin + num_int
+            @assert node.local_tightenings <= num_bin + num_int
+            @assert node.global_tightenings <= num_bin + num_int
+        end
+        # update current_node_id
+        if !Bonobo.terminated(tree)
+            tree.root.current_node_id[] =
+                Bonobo.get_next_node(tree, tree.options.traverse_strategy).id
+        end
+
+        if Bonobo.terminated(tree)
+            Bonobo.sort_solutions!(tree.solutions, tree.sense)
+            x = Bonobo.get_solution(tree)
+            # x can be nothing if the user supplied a custom domain oracle and the time limit is reached
+            if x === nothing
+                @assert tree.root.problem.solving_stage == TIME_LIMIT_REACHED
+            end
+            primal_value = x !== nothing ? tree.root.problem.f(x) : Inf
+            # deactivate postsolve if there is no solution
+            tree.root.options[:usePostsolve] =
+                x === nothing ? false : tree.root.options[:usePostsolve]
+
+            # TODO: here we need to calculate the actual state
+
+            # If the tree is empty, incumbent and solution should be the same!
+            if isempty(tree.nodes)
+                @assert isapprox(tree.incumbent, primal_value)
+            end
+
+            result[:number_nodes] = tree.num_nodes
+            result[:lmo_calls] = tree.root.problem.tlmo.ncalls
+            result[:heu_lmo_calls] = tree.root.options[:heu_ncalls]
+            result[:list_num_nodes] = list_num_nodes_cb
+            result[:list_lmo_calls_acc] = list_lmo_calls_cb
+            result[:list_active_set_size] = list_active_set_size_cb
+            result[:list_discarded_set_size] = list_discarded_set_size_cb
+            result[:list_lb] = list_lb_cb
+            result[:list_ub] = list_ub_cb
+            result[:list_time] = list_time_cb
+            result[:lmo_calls_per_layer] = lmo_calls_per_layer
+            result[:active_set_size_per_layer] = active_set_size_per_layer
+            result[:discarded_set_size_per_layer] = discarded_set_size_per_layer
+            result[:node_level] = node_level
+            result[:global_tightenings] = global_tightenings
+            result[:local_tightenings] = local_tightenings
+            result[:local_potential_tightenings] = local_potential_tightenings
+
+            if verbose
+                FrankWolfe.print_callback(headers, format_string, print_footer=true)
+                println()
+            end
+        end
+    end
+end
diff --git a/src/custom_bonobo.jl b/src/custom_bonobo.jl
index dcce29444..d50b3b24b 100644
--- a/src/custom_bonobo.jl
+++ b/src/custom_bonobo.jl
@@ -74,6 +74,25 @@ function Bonobo.optimize!(
         Bonobo.branch!(tree, node)
         callback(tree, node)
     end
+    # To make sure that we collect the statistics in case the time limit is reached.
+    if !haskey(tree.root.result, :global_tightenings)
+        y = Bonobo.get_solution(tree)
+        vertex_storage = FrankWolfe.DeletedVertexStorage(typeof(y)[], 1)
+        dummy_node = FrankWolfeNode(
+            NodeInfo(-1, Inf, Inf),
+            FrankWolfe.ActiveSet([(1.0, y)]),
+            vertex_storage,
+            IntegerBounds(),
+            1,
+            1e-3,
+            Millisecond(0),
+            0,
+            0,
+            0,
+            0.0,
+        )
+        callback(tree, dummy_node, node_infeasible=true)
+    end
     return Bonobo.sort_solutions!(tree.solutions, tree.sense)
 end
 
diff --git a/src/heuristics.jl b/src/heuristics.jl
index 3a5d7b3a6..c30392a8c 100644
--- a/src/heuristics.jl
+++ b/src/heuristics.jl
@@ -53,7 +53,7 @@ function run_heuristics(tree, x, heuristic_list; rng=Random.GLOBAL_RNG)
                 min_val = Inf
                 min_idx = -1
                 for (i, x_heu) in enumerate(list_x_heu)
-                    feasible = check_feasibility ? is_linear_feasible(tree.root.problem.tlmo, x_heu) && is_integer_feasible(tree, x_heu) : true
+                    feasible = check_feasibility ? is_linear_feasible(tree.root.problem.tlmo, x_heu) && is_integer_feasible(tree, x_heu) : false
                     if feasible
                         val = tree.root.problem.f(x_heu)
                         if val < min_val
@@ -159,13 +159,12 @@ function probability_rounding(tree::Bonobo.BnBTree, tlmo::Boscia.TimeTrackingLMO
         line_search=tree.root.options[:lineSearch],
         lazy=tree.root.options[:lazy],
         lazy_tolerance=tree.root.options[:lazy_tolerance],
-        add_dropped_vertices=tree.root.options[:use_shadow_set],
-        use_extra_vertex_storage=tree.root.options[:use_shadow_set],
-        extra_vertex_storage=node.discarded_vertices,
         callback=tree.root.options[:callback],
         verbose=tree.root.options[:fwVerbose],
     )
 
+    @assert sum(isapprox.(x_rounded[tlmo.blmo.int_vars], round.(x_rounded[tlmo.blmo.int_vars]))) == length(tlmo.blmo.int_vars) "$(sum(isapprox.(x_rounded[tlmo.blmo.int_vars], round.(x_rounded[tlmo.blmo.int_vars])))) == $(length(tlmo.blmo.int_vars)) $(x_rounded[tlmo.blmo.int_vars])"
+
     # reset LMO to node state
     build_LMO(tlmo, tree.root.problem.integer_variable_bounds, original_bounds, tlmo.blmo.int_vars)
     
diff --git a/src/interface.jl b/src/interface.jl
index 74d0f30fe..772d9c9a0 100644
--- a/src/interface.jl
+++ b/src/interface.jl
@@ -43,6 +43,11 @@ dual_tightening        - whether to use dual tightening techniques (make sure yo
 global_dual_tightening - dual tightening maintained globally valid (when new solutions are found)
 bnb_callback           - an optional callback called at every node of the tree, for example for heuristics
 strong_convexity       - strong convexity of the function, used for tighter dual bound at every node
+sharpness_constant     - the constant M > 0 for (θ, M)-sharpness.
+                         f is (θ, M)-sharpness: f satisfies
+                         min_{x* ∈ X*} || x - x* || ≤ M (f(x) - f*)
+                         where X* is the set minimizer of f. 
+sharpness_exponent     - the exponent θ ∈ [0, 1/2] for (θ, M)-sharpness.
 domain_oracle          - For a point x: returns true if x is in the domain of f, else false. Per default is true.
                          In case of the non trivial domain oracle, the starting point has to be feasible for f. Also, depending 
                          on the Line Search method, you might have to provide the domain oracle to it, too.
@@ -91,6 +96,8 @@ function solve(
     global_dual_tightening=true,
     bnb_callback=nothing,
     strong_convexity=0.0,
+    sharpness_constant = 0.0,
+    sharpness_exponent = Inf,
     domain_oracle=x -> true,
     start_solution=nothing,
     fw_verbose=false,
@@ -152,8 +159,8 @@ function solve(
         for a in active_set.atoms
             @assert is_linear_feasible(blmo, a)
         end
-        v = active_set.atoms[1]
         x = FrankWolfe.compute_active_set_iterate!(active_set)
+        v = x
         @assert isfinite(f(x))
     end
     vertex_storage = FrankWolfe.DeletedVertexStorage(typeof(v)[], 1)
@@ -206,6 +213,8 @@ function solve(
                 :max_fw_iter => max_fw_iter,
                 :print_iter => print_iter,
                 :strong_convexity => strong_convexity,
+                :sharpness_constant => sharpness_constant,
+                :sharpness_exponent => sharpness_exponent,
                 :time_limit => time_limit,
                 :usePostsolve => use_postsolve,
                 :variant => variant,
@@ -215,6 +224,7 @@ function solve(
                 :max_clean_iter => max_clean_iter,
                 :clean_solutions => clean_solutions,
             ),
+            result=Dict{Symbol,Any}(),
         ),
         branch_strategy=branching_strategy,
         dual_gap_limit=rel_dual_gap,
@@ -259,7 +269,6 @@ function solve(
     list_discarded_set_size_cb = Int[]
     fw_iterations = Int[]
     node_level = Int[]
-    result = Dict{Symbol,Any}()
     lmo_calls_per_layer = Vector{Vector{Int}}()
     active_set_size_per_layer = Vector{Vector{Int}}()
     discarded_set_size_per_layer = Vector{Vector{Int}}()
@@ -280,7 +289,7 @@ function solve(
         fw_iterations,
         list_active_set_size_cb,
         list_discarded_set_size_cb,
-        result,
+        tree.root.result,
         lmo_calls_per_layer,
         active_set_size_per_layer,
         discarded_set_size_per_layer,
@@ -300,7 +309,7 @@ function solve(
 
     Bonobo.optimize!(tree; callback=bnb_callback)
 
-    x = postsolve(tree, result, time_ref, verbose, max_iteration_post)
+    x = postsolve(tree, tree.root.result, time_ref, verbose, max_iteration_post)
 
     # Check solution and polish
     x_polished = x
@@ -323,255 +332,7 @@ function solve(
     end
     println() # cleaner output
 
-    return x, tree.root.problem.tlmo, result
-end
-
-"""
-Output of Boscia
-
-    iter :          current iteration of Boscia
-    node id :       current node id
-    lower bound :   tree_lb(tree)
-    incumbent :     tree.incumbent
-    gap :           tree.incumbent-tree_lb(tree)
-    rel. gap :      dual_gap/tree.incumbent
-    time :          total time of Boscia
-    time/nodes :    average time per node
-    FW time :       time spent in FW 
-    LMO time :      time used by LMO
-    LMO calls :     number of compute_extreme_point calls in FW
-    FW iterations : number of iterations in FW
-"""
-function build_bnb_callback(
-    tree,
-    time_ref,
-    list_lb_cb,
-    list_ub_cb,
-    list_time_cb,
-    list_num_nodes_cb,
-    list_lmo_calls_cb,
-    verbose,
-    fw_iterations,
-    list_active_set_size_cb,
-    list_discarded_set_size_cb,
-    result,
-    lmo_calls_per_layer,
-    active_set_size_per_layer,
-    discarded_set_size_per_layer,
-    node_level,
-    baseline_callback,
-    local_tightenings,
-    global_tightenings,
-    local_potential_tightenings,
-    num_bin,
-    num_int,
-)
-    iteration = 0
-
-    headers = [
-        " ",
-        "Iter",
-        "Open",
-        "Bound",
-        "Incumbent",
-        "Gap (abs)",
-        "Gap (rel)",
-        "Time (s)",
-        "Nodes/sec",
-        "FW (ms)",
-        "LMO (ms)",
-        "LMO (calls c)",
-        "FW (its)",
-        "#activeset",
-        "#shadow",
-    ]
-    format_string = "%1s %5i %5i %14e %14e %14e %14e %14e %14e %12i %10i %14i %10i %8i %8i\n"
-    print_iter = get(tree.root.options, :print_iter, 100)
-
-    if verbose
-        FrankWolfe.print_callback(headers, format_string, print_header=true)
-    end
-    return function callback(
-        tree,
-        node;
-        worse_than_incumbent=false,
-        node_infeasible=false,
-        lb_update=false,
-    )
-        if baseline_callback !== nothing
-            baseline_callback(
-                tree,
-                node,
-                worse_than_incumbent=worse_than_incumbent,
-                node_infeasible=node_infeasible,
-                lb_update=lb_update,
-            )
-        end
-        if !node_infeasible
-            #update lower bound
-            if lb_update == true
-                tree.node_queue[node.id] = (node.lb, node.id)
-                _, prio = peek(tree.node_queue)
-                @assert tree.lb <= prio[1]
-                tree.lb = min(minimum([prio[2][1] for prio in tree.node_queue]), tree.incumbent)
-            end
-            push!(list_ub_cb, tree.incumbent)
-            push!(list_num_nodes_cb, tree.num_nodes)
-            push!(node_level, node.level)
-            iteration += 1
-            if tree.lb == -Inf && isempty(tree.nodes)
-                tree.lb = node.lb
-            end
-
-            time = float(Dates.value(Dates.now() - time_ref))
-            push!(list_time_cb, time)
-
-            if tree.root.options[:time_limit] < Inf
-                if time / 1000.0 ≥ tree.root.options[:time_limit]
-                    @assert tree.root.problem.solving_stage == SOLVING
-                    tree.root.problem.solving_stage = TIME_LIMIT_REACHED
-                end
-            end
-
-            fw_time = Dates.value(node.fw_time)
-            fw_iter = if !isempty(fw_iterations)
-                fw_iterations[end]
-            else
-                0
-            end
-            if !isempty(tree.root.problem.tlmo.optimizing_times)
-                LMO_time = sum(1000 * tree.root.problem.tlmo.optimizing_times)
-                empty!(tree.root.problem.tlmo.optimizing_times)
-            else
-                LMO_time = 0
-            end
-            LMO_calls_c = tree.root.problem.tlmo.ncalls
-            push!(list_lmo_calls_cb, copy(LMO_calls_c))
-
-            if !isempty(tree.node_queue)
-                p_lb = tree.lb
-                tree.lb = min(minimum([prio[2][1] for prio in tree.node_queue]), tree.incumbent)
-                @assert p_lb <= tree.lb + tree.root.options[:dual_gap] "p_lb <= tree.lb + tree.root.options[:dual_gap] $(p_lb) <= $(tree.lb + tree.root.options[:dual_gap])"
-            end
-            # correct lower bound if necessary
-            tree.lb = tree_lb(tree)
-            dual_gap = tree.incumbent - tree_lb(tree)
-            push!(list_lb_cb, tree_lb(tree))
-            active_set_size = length(node.active_set)
-            discarded_set_size = length(node.discarded_vertices.storage)
-            push!(list_active_set_size_cb, active_set_size)
-            push!(list_discarded_set_size_cb, discarded_set_size)
-            nodes_left = length(tree.nodes)
-            if tree.root.updated_incumbent[]
-                incumbent_updated = "*"
-            else
-                incumbent_updated = " "
-            end
-            if verbose && (
-                mod(iteration, print_iter) == 0 ||
-                iteration == 1 ||
-                Bonobo.terminated(tree) ||
-                tree.root.updated_incumbent[]
-            )
-                if (mod(iteration, print_iter * 40) == 0)
-                    FrankWolfe.print_callback(headers, format_string, print_header=true)
-                end
-                FrankWolfe.print_callback(
-                    (
-                        incumbent_updated,
-                        iteration,
-                        nodes_left,
-                        tree_lb(tree),
-                        tree.incumbent,
-                        dual_gap,
-                        relative_gap(tree.incumbent, tree_lb(tree)),
-                        time / 1000.0,
-                        tree.num_nodes / time * 1000.0,
-                        fw_time,
-                        LMO_time,
-                        tree.root.problem.tlmo.ncalls,
-                        fw_iter,
-                        active_set_size,
-                        discarded_set_size,
-                    ),
-                    format_string,
-                    print_header=false,
-                )
-                tree.root.updated_incumbent[] = false
-            end
-            # lmo calls per layer
-            if length(list_lmo_calls_cb) > 1
-                LMO_calls = list_lmo_calls_cb[end] - list_lmo_calls_cb[end-1]
-            else
-                LMO_calls = list_lmo_calls_cb[end]
-            end
-            if length(lmo_calls_per_layer) < node.level
-                push!(lmo_calls_per_layer, [LMO_calls])
-                push!(active_set_size_per_layer, [active_set_size])
-                push!(discarded_set_size_per_layer, [discarded_set_size])
-            else
-                push!(lmo_calls_per_layer[node.level], LMO_calls)
-                push!(active_set_size_per_layer[node.level], active_set_size)
-                push!(discarded_set_size_per_layer[node.level], discarded_set_size)
-            end
-
-            # add tightenings
-            push!(global_tightenings, node.global_tightenings)
-            push!(local_tightenings, node.local_tightenings)
-            push!(local_potential_tightenings, node.local_potential_tightenings)
-            @assert node.local_potential_tightenings <= num_bin + num_int
-            @assert node.local_tightenings <= num_bin + num_int
-            @assert node.global_tightenings <= num_bin + num_int
-        end
-        # update current_node_id
-        if !Bonobo.terminated(tree)
-            tree.root.current_node_id[] =
-                Bonobo.get_next_node(tree, tree.options.traverse_strategy).id
-        end
-
-        if Bonobo.terminated(tree)
-            Bonobo.sort_solutions!(tree.solutions, tree.sense)
-            x = Bonobo.get_solution(tree)
-            # x can be nothing if the user supplied a custom domain oracle and the time limit is reached
-            if x === nothing
-                @assert tree.root.problem.solving_stage == TIME_LIMIT_REACHED
-            end
-            primal_value = x !== nothing ? tree.root.problem.f(x) : Inf
-            # deactivate postsolve if there is no solution
-            tree.root.options[:usePostsolve] =
-                x === nothing ? false : tree.root.options[:usePostsolve]
-
-            # TODO: here we need to calculate the actual state
-
-            # If the tree is empty, incumbent and solution should be the same!
-            if isempty(tree.nodes)
-                @assert isapprox(tree.incumbent, primal_value)
-            end
-
-            result[:number_nodes] = tree.num_nodes
-            result[:lmo_calls] = tree.root.problem.tlmo.ncalls
-            result[:heu_lmo_calls] = tree.root.options[:heu_ncalls]
-            result[:list_num_nodes] = list_num_nodes_cb
-            result[:list_lmo_calls_acc] = list_lmo_calls_cb
-            result[:list_active_set_size] = list_active_set_size_cb
-            result[:list_discarded_set_size] = list_discarded_set_size_cb
-            result[:list_lb] = list_lb_cb
-            result[:list_ub] = list_ub_cb
-            result[:list_time] = list_time_cb
-            result[:lmo_calls_per_layer] = lmo_calls_per_layer
-            result[:active_set_size_per_layer] = active_set_size_per_layer
-            result[:discarded_set_size_per_layer] = discarded_set_size_per_layer
-            result[:node_level] = node_level
-            result[:global_tightenings] = global_tightenings
-            result[:local_tightenings] = local_tightenings
-            result[:local_potential_tightenings] = local_potential_tightenings
-
-            if verbose
-                FrankWolfe.print_callback(headers, format_string, print_footer=true)
-                println()
-            end
-        end
-    end
+    return x, tree.root.problem.tlmo, tree.root.result
 end
 
 """
diff --git a/src/managed_blmo.jl b/src/managed_blmo.jl
index e13267355..c64967e28 100644
--- a/src/managed_blmo.jl
+++ b/src/managed_blmo.jl
@@ -269,6 +269,8 @@ function solve(
     global_dual_tightening=true,
     bnb_callback=nothing,
     strong_convexity=0.0,
+    sharpness_constant = 0.0,
+    sharpness_exponent = Inf,
     domain_oracle=x -> true,
     start_solution=nothing,
     fw_verbose=false,
@@ -308,6 +310,8 @@ function solve(
         global_dual_tightening=global_dual_tightening,
         bnb_callback=bnb_callback,
         strong_convexity=strong_convexity,
+        sharpness_constant=sharpness_constant,
+        sharpness_exponent=sharpness_exponent,
         domain_oracle=domain_oracle,
         start_solution=start_solution,
         fw_verbose=fw_verbose,
diff --git a/src/node.jl b/src/node.jl
index 170960ebd..e7eccc8e0 100644
--- a/src/node.jl
+++ b/src/node.jl
@@ -175,52 +175,11 @@ function Bonobo.get_branching_nodes_info(tree::Bonobo.BnBTree, node::FrankWolfeN
     elseif domain_left # x_left in domain
         [node_info_left]
     else
-        warn("No childern nodes can be created.")
+        @warn "No childern nodes can be created."
     end
     return nodes
 end
 
-"""
-Use strong convexity to potentially remove one of the children nodes 
-"""
-function prune_children(tree, node, lower_bound_base, x, vidx)
-    prune_left = false
-    prune_right = false
-
-    μ = tree.root.options[:strong_convexity]
-    if μ > 0
-        @debug "Using strong convexity $μ"
-        for j in tree.root.problem.integer_variables
-            if vidx == j
-                continue
-            end
-            lower_bound_base += μ / 2 * min((x[j] - floor(x[j]))^2, (ceil(x[j]) - x[j])^2)
-        end
-        new_bound_left = lower_bound_base + μ / 2 * (x[vidx] - floor(x[vidx]))^2
-        new_bound_right = lower_bound_base + μ / 2 * (ceil(x[vidx]) - x[vidx])^2
-        if new_bound_left > tree.incumbent
-            @debug "prune left, from $(node.lb) -> $new_bound_left, ub $(tree.incumbent), lb $(node.lb)"
-            prune_left = true
-        end
-        if new_bound_right > tree.incumbent
-            @debug "prune right, from $(node.lb) -> $new_bound_right, ub $(tree.incumbent), lb $(node.lb)"
-            prune_right = true
-        end
-        @assert !(
-            (new_bound_left > tree.incumbent + tree.root.options[:dual_gap]) &&
-            (new_bound_right > tree.incumbent + tree.root.options[:dual_gap])
-        ) "both sides should not be pruned"
-    end
-
-    # If both nodes are pruned, when one of them has to be equal to the incumbent.
-    # Thus, we have proof of optimality by strong convexity.
-    if prune_left && prune_right
-        tree.lb = min(new_bound_left, new_bound_right)
-    end
-
-    return prune_left, prune_right
-end
-
 """
 Computes the relaxation at that node
 """
@@ -311,7 +270,7 @@ function Bonobo.evaluate_node!(tree::Bonobo.BnBTree, node::FrankWolfeNode)
 
     lower_bound = primal - dual_gap
     # improvement of the lower bound using strong convexity
-    lower_bound = tightening_strong_convexity(tree, x, lower_bound)
+    lower_bound = tightening_lowerbound(tree, node, x, lower_bound)
 
     # Found an upper bound
     if is_integer_feasible(tree, x)
@@ -332,214 +291,6 @@ function Bonobo.evaluate_node!(tree::Bonobo.BnBTree, node::FrankWolfeNode)
     return lower_bound, NaN
 end
 
-"""
-Tightening of the bounds at node level. Children node inherit the updated bounds.
-"""
-function dual_tightening(tree, node, x, dual_gap)
-    if tree.root.options[:dual_tightening] && isfinite(tree.incumbent)
-        grad = similar(x)
-        tree.root.problem.g(grad, x)
-        num_tightenings = 0
-        num_potential_tightenings = 0
-        μ = tree.root.options[:strong_convexity]
-        for j in tree.root.problem.integer_variables
-            lb_global = get(tree.root.problem.integer_variable_bounds, (j, :greaterthan), -Inf)
-            ub_global = get(tree.root.problem.integer_variable_bounds, (j, :lessthan), Inf)
-            lb = get(node.local_bounds.lower_bounds, j, lb_global)
-            ub = get(node.local_bounds.upper_bounds, j, ub_global)
-            @assert lb >= lb_global
-            @assert ub <= ub_global
-            if lb ≈ ub
-                # variable already fixed
-                continue
-            end
-            gj = grad[j]
-            safety_tolerance = 2.0
-            rhs =
-                tree.incumbent - tree.root.problem.f(x) +
-                safety_tolerance * dual_gap +
-                sqrt(eps(tree.incumbent))
-            if ≈(x[j], lb, atol=tree.options.atol, rtol=tree.options.rtol)
-                if !isapprox(gj, 0, atol=1e-5)
-                    num_potential_tightenings += 1
-                end
-                if gj > 0
-                    Mlb = 0
-                    bound_tightened = true
-                    @debug "starting tightening ub $(rhs)"
-                    while 0.99 * (Mlb * gj + μ / 2 * Mlb^2) <= rhs
-                        Mlb += 1
-                        if lb + Mlb - 1 == ub
-                            bound_tightened = false
-                            break
-                        end
-                    end
-                    if bound_tightened
-                        new_bound = lb + Mlb - 1
-                        @debug "found UB tightening $ub -> $new_bound"
-                        node.local_bounds[j, :lessthan] = new_bound
-                        num_tightenings += 1
-                        if haskey(tree.root.problem.integer_variable_bounds, (j, :lessthan))
-                            @assert node.local_bounds[j, :lessthan] <=
-                                    tree.root.problem.integer_variable_bounds[j, :lessthan]
-                        end
-                    end
-                end
-            elseif ≈(x[j], ub, atol=tree.options.atol, rtol=tree.options.rtol)
-                if !isapprox(gj, 0, atol=1e-5)
-                    num_potential_tightenings += 1
-                end
-                if gj < 0
-                    Mub = 0
-                    bound_tightened = true
-                    @debug "starting tightening lb $(rhs)"
-                    while -0.99 * (Mub * gj + μ / 2 * Mub^2) <= rhs
-                        Mub += 1
-                        if ub - Mub + 1 == lb
-                            bound_tightened = false
-                            break
-                        end
-                    end
-                    if bound_tightened
-                        new_bound = ub - Mub + 1
-                        @debug "found LB tightening $lb -> $new_bound"
-                        node.local_bounds[j, :greaterthan] = new_bound
-                        num_tightenings += 1
-                        if haskey(tree.root.problem.integer_variable_bounds, (j, :greaterthan))
-                            @assert node.local_bounds[j, :greaterthan] >=
-                                    tree.root.problem.integer_variable_bounds[j, :greaterthan]
-                        end
-                    end
-                end
-            end
-        end
-        @debug "# tightenings $num_tightenings"
-        node.local_tightenings = num_tightenings
-        node.local_potential_tightenings = num_potential_tightenings
-    end
-end
-
-"""
-Save the gradient of the root solution (i.e. the relaxed solution) and the 
-corresponding lower and upper bounds.
-"""
-function store_data_global_tightening(tree, node, x, dual_gap)
-    if tree.root.options[:global_dual_tightening] && node.std.id == 1
-        @debug "storing root node info for tightening"
-        grad = similar(x)
-        tree.root.problem.g(grad, x)
-        safety_tolerance = 2.0
-        tree.root.global_tightening_rhs[] = -tree.root.problem.f(x) + safety_tolerance * dual_gap
-        for j in tree.root.problem.integer_variables
-            if haskey(tree.root.problem.integer_variable_bounds.upper_bounds, j)
-                ub = tree.root.problem.integer_variable_bounds[j, :lessthan]
-                if ≈(x[j], ub, atol=tree.options.atol, rtol=tree.options.rtol) && grad[j] < 0
-                    tree.root.global_tightening_root_info.upper_bounds[j] = (grad[j], ub)
-                end
-            end
-            if haskey(tree.root.problem.integer_variable_bounds.lower_bounds, j)
-                lb = tree.root.problem.integer_variable_bounds[j, :greaterthan]
-                if ≈(x[j], lb, atol=tree.options.atol, rtol=tree.options.rtol) && grad[j] > 0
-                    tree.root.global_tightening_root_info.lower_bounds[j] = (grad[j], lb)
-                end
-            end
-        end
-    end
-end
-
-"""
-Use the gradient of the root node to tighten the global bounds.
-"""
-function global_tightening(tree, node)
-    # new incumbent: check global fixings
-    if tree.root.options[:global_dual_tightening] && tree.root.updated_incumbent[]
-        num_tightenings = 0
-        rhs = tree.incumbent + tree.root.global_tightening_rhs[]
-        @assert isfinite(rhs)
-        for (j, (gj, lb)) in tree.root.global_tightening_root_info.lower_bounds
-            ub = get(tree.root.problem.integer_variable_bounds.upper_bounds, j, Inf)
-            ub_new = get(tree.root.global_tightening_root_info.upper_bounds, j, Inf)
-            ub = min(ub, ub_new)
-            Mlb = 0
-            bound_tightened = true
-            lb = lb
-            while Mlb * gj <= rhs
-                Mlb += 1
-                if lb + Mlb - 1 == ub
-                    bound_tightened = false
-                    break
-                end
-            end
-            if bound_tightened
-                new_bound = lb + Mlb - 1
-                @debug "found global UB tightening $ub -> $new_bound"
-                if haskey(tree.root.global_tightenings.upper_bounds, j)
-                    if tree.root.global_tightenings.upper_bounds[j] != new_bound
-                        num_tightenings += 1
-                    end
-                else
-                    num_tightenings += 1
-                end
-                tree.root.global_tightenings.upper_bounds[j] = new_bound
-            end
-        end
-        for (j, (gj, ub)) in tree.root.global_tightening_root_info.upper_bounds
-            lb = get(tree.root.problem.integer_variable_bounds.lower_bounds, j, -Inf)
-            lb_new = get(tree.root.global_tightening_root_info.lower_bounds, j, -Inf)
-            lb = max(lb, lb_new)
-            Mub = 0
-            bound_tightened = true
-            ub = ub
-            while -Mub * gj <= rhs
-                Mub += 1
-                if ub - Mub + 1 == lb
-                    bound_tightened = false
-                    break
-                end
-            end
-            if bound_tightened
-                new_bound = ub - Mub + 1
-                @debug "found global LB tightening $lb -> $new_bound"
-                if haskey(tree.root.global_tightenings.lower_bounds, j)
-                    if tree.root.global_tightenings.lower_bounds[j] != new_bound
-                        num_tightenings += 1
-                    end
-                else
-                    num_tightenings += 1
-                end
-                tree.root.global_tightenings.lower_bounds[j] = new_bound
-            end
-        end
-        node.global_tightenings = num_tightenings
-    end
-end
-
-"""
-Tighten the lower bound using strong convexity of the objective.
-"""
-function tightening_strong_convexity(tree, x, lower_bound)
-    μ = tree.root.options[:strong_convexity]
-    if μ > 0
-        @debug "Using strong convexity $μ"
-        strong_convexity_bound = lower_bound
-        num_fractional = 0
-        for j in tree.root.problem.integer_variables
-            if x[j] > floor(x[j]) + 1e-6 && x[j] < ceil(x[j]) - 1e-6
-                num_fractional += 1
-                new_left_increment = μ / 2 * (x[j] - floor(x[j]))^2
-                new_right_increment = μ / 2 * (ceil(x[j]) - x[j])^2
-                new_increment = min(new_left_increment, new_right_increment)
-                strong_convexity_bound += new_increment
-            end
-        end
-        @debug "Strong convexity: $lower_bound -> $strong_convexity_bound"
-        @assert num_fractional == 0 || strong_convexity_bound > lower_bound
-        lower_bound = strong_convexity_bound
-    end
-    return lower_bound
-end
-
-
 """
 Returns the solution vector of the relaxed problem at the node
 """
diff --git a/src/tightenings.jl b/src/tightenings.jl
new file mode 100644
index 000000000..2cde12a4f
--- /dev/null
+++ b/src/tightenings.jl
@@ -0,0 +1,300 @@
+"""
+Tightening of the bounds at node level. Children node inherit the updated bounds.
+"""
+function dual_tightening(tree, node, x, dual_gap)
+    if tree.root.options[:dual_tightening] && isfinite(tree.incumbent)
+        grad = similar(x)
+        tree.root.problem.g(grad, x)
+        num_tightenings = 0
+        num_potential_tightenings = 0
+        μ = tree.root.options[:strong_convexity]
+        for j in tree.root.problem.integer_variables
+            lb_global = get(tree.root.problem.integer_variable_bounds, (j, :greaterthan), -Inf)
+            ub_global = get(tree.root.problem.integer_variable_bounds, (j, :lessthan), Inf)
+            lb = get(node.local_bounds.lower_bounds, j, lb_global)
+            ub = get(node.local_bounds.upper_bounds, j, ub_global)
+            @assert lb >= lb_global
+            @assert ub <= ub_global
+            if lb ≈ ub
+                # variable already fixed
+                continue
+            end
+            gj = grad[j]
+            safety_tolerance = 2.0
+            rhs =
+                tree.incumbent - tree.root.problem.f(x) +
+                safety_tolerance * dual_gap +
+                sqrt(eps(tree.incumbent))
+            if ≈(x[j], lb, atol=tree.options.atol, rtol=tree.options.rtol)
+                if !isapprox(gj, 0, atol=1e-5)
+                    num_potential_tightenings += 1
+                end
+                if gj > 0
+                    Mlb = 0
+                    bound_tightened = true
+                    @debug "starting tightening ub $(rhs)"
+                    while 0.99 * (Mlb * gj + μ / 2 * Mlb^2) <= rhs
+                        Mlb += 1
+                        if lb + Mlb - 1 == ub
+                            bound_tightened = false
+                            break
+                        end
+                    end
+                    if bound_tightened
+                        new_bound = lb + Mlb - 1
+                        @debug "found UB tightening $ub -> $new_bound"
+                        node.local_bounds[j, :lessthan] = new_bound
+                        num_tightenings += 1
+                        if haskey(tree.root.problem.integer_variable_bounds, (j, :lessthan))
+                            @assert node.local_bounds[j, :lessthan] <=
+                                    tree.root.problem.integer_variable_bounds[j, :lessthan]
+                        end
+                    end
+                end
+            elseif ≈(x[j], ub, atol=tree.options.atol, rtol=tree.options.rtol)
+                if !isapprox(gj, 0, atol=1e-5)
+                    num_potential_tightenings += 1
+                end
+                if gj < 0
+                    Mub = 0
+                    bound_tightened = true
+                    @debug "starting tightening lb $(rhs)"
+                    while -0.99 * (Mub * gj + μ / 2 * Mub^2) <= rhs
+                        Mub += 1
+                        if ub - Mub + 1 == lb
+                            bound_tightened = false
+                            break
+                        end
+                    end
+                    if bound_tightened
+                        new_bound = ub - Mub + 1
+                        @debug "found LB tightening $lb -> $new_bound"
+                        node.local_bounds[j, :greaterthan] = new_bound
+                        num_tightenings += 1
+                        if haskey(tree.root.problem.integer_variable_bounds, (j, :greaterthan))
+                            @assert node.local_bounds[j, :greaterthan] >=
+                                    tree.root.problem.integer_variable_bounds[j, :greaterthan]
+                        end
+                    end
+                end
+            end
+        end
+        @debug "# tightenings $num_tightenings"
+        node.local_tightenings = num_tightenings
+        node.local_potential_tightenings = num_potential_tightenings
+    end
+end
+
+"""
+Save the gradient of the root solution (i.e. the relaxed solution) and the 
+corresponding lower and upper bounds.
+"""
+function store_data_global_tightening(tree, node, x, dual_gap)
+    if tree.root.options[:global_dual_tightening] && node.std.id == 1
+        @debug "storing root node info for tightening"
+        grad = similar(x)
+        tree.root.problem.g(grad, x)
+        safety_tolerance = 2.0
+        tree.root.global_tightening_rhs[] = -tree.root.problem.f(x) + safety_tolerance * dual_gap
+        for j in tree.root.problem.integer_variables
+            if haskey(tree.root.problem.integer_variable_bounds.upper_bounds, j)
+                ub = tree.root.problem.integer_variable_bounds[j, :lessthan]
+                if ≈(x[j], ub, atol=tree.options.atol, rtol=tree.options.rtol) && grad[j] < 0
+                    tree.root.global_tightening_root_info.upper_bounds[j] = (grad[j], ub)
+                end
+            end
+            if haskey(tree.root.problem.integer_variable_bounds.lower_bounds, j)
+                lb = tree.root.problem.integer_variable_bounds[j, :greaterthan]
+                if ≈(x[j], lb, atol=tree.options.atol, rtol=tree.options.rtol) && grad[j] > 0
+                    tree.root.global_tightening_root_info.lower_bounds[j] = (grad[j], lb)
+                end
+            end
+        end
+    end
+end
+
+"""
+Use the gradient of the root node to tighten the global bounds.
+"""
+function global_tightening(tree, node)
+    # new incumbent: check global fixings
+    if tree.root.options[:global_dual_tightening] && tree.root.updated_incumbent[]
+        num_tightenings = 0
+        rhs = tree.incumbent + tree.root.global_tightening_rhs[]
+        @assert isfinite(rhs)
+        for (j, (gj, lb)) in tree.root.global_tightening_root_info.lower_bounds
+            ub = get(tree.root.problem.integer_variable_bounds.upper_bounds, j, Inf)
+            ub_new = get(tree.root.global_tightening_root_info.upper_bounds, j, Inf)
+            ub = min(ub, ub_new)
+            Mlb = 0
+            bound_tightened = true
+            lb = lb
+            while Mlb * gj <= rhs
+                Mlb += 1
+                if lb + Mlb - 1 == ub
+                    bound_tightened = false
+                    break
+                end
+            end
+            if bound_tightened
+                new_bound = lb + Mlb - 1
+                @debug "found global UB tightening $ub -> $new_bound"
+                if haskey(tree.root.global_tightenings.upper_bounds, j)
+                    if tree.root.global_tightenings.upper_bounds[j] != new_bound
+                        num_tightenings += 1
+                    end
+                else
+                    num_tightenings += 1
+                end
+                tree.root.global_tightenings.upper_bounds[j] = new_bound
+            end
+        end
+        for (j, (gj, ub)) in tree.root.global_tightening_root_info.upper_bounds
+            lb = get(tree.root.problem.integer_variable_bounds.lower_bounds, j, -Inf)
+            lb_new = get(tree.root.global_tightening_root_info.lower_bounds, j, -Inf)
+            lb = max(lb, lb_new)
+            Mub = 0
+            bound_tightened = true
+            ub = ub
+            while -Mub * gj <= rhs
+                Mub += 1
+                if ub - Mub + 1 == lb
+                    bound_tightened = false
+                    break
+                end
+            end
+            if bound_tightened
+                new_bound = ub - Mub + 1
+                @debug "found global LB tightening $lb -> $new_bound"
+                if haskey(tree.root.global_tightenings.lower_bounds, j)
+                    if tree.root.global_tightenings.lower_bounds[j] != new_bound
+                        num_tightenings += 1
+                    end
+                else
+                    num_tightenings += 1
+                end
+                tree.root.global_tightenings.lower_bounds[j] = new_bound
+            end
+        end
+        node.global_tightenings = num_tightenings
+    end
+end
+
+"""
+Tighten the lower bound using strong convexity and/or sharpness of the objective.
+"""
+function tightening_lowerbound(tree, node, x, lower_bound)
+    μ = tree.root.options[:strong_convexity]
+    M = tree.root.options[:sharpness_constant]
+    θ = tree.root.options[:sharpness_exponent]
+
+    if μ > 0 || (M > 0 && θ != Inf)
+        @debug "Tightening lower bound using strong convexity $μ and/or sharpness ($θ, $M)"
+        num_fractional = 0
+        bound_improvement = 0.0
+        for j in tree.root.problem.integer_variables
+            if x[j] > floor(x[j]) + 1e-6 && x[j] < ceil(x[j]) - 1e-6
+                num_fractional += 1
+                new_left_increment = (x[j] - floor(x[j]))^2
+                new_right_increment = (ceil(x[j]) - x[j])^2
+                new_increment = min(new_left_increment, new_right_increment)
+                bound_improvement += new_increment
+            end
+        end
+        strong_convexity_bound = lower_bound
+        sharpness_bound = -Inf
+
+        # strong convexity
+        if μ > 0 
+            @debug "Using strong convexity $μ"
+            strong_convexity_bound += μ / 2 * bound_improvement
+            @debug "Strong convexity: $lower_bound -> $strong_convexity_bound"
+            @assert num_fractional == 0 || strong_convexity_bound > lower_bound 
+        end
+
+        # sharpness
+        if M > 0 && θ != Inf
+
+            @debug "Using sharpness θ=$θ and M=$M"
+            fx = tree.root.problem.f(x)
+
+            sharpness_bound = M^(- 1 / θ) * 1/ 2 * (sqrt(bound_improvement) - M / 2 * node.dual_gap^θ)^(1 / θ) + fx - node.dual_gap
+            @debug "Sharpness: $lower_bound -> $sharpness_bound"
+            @assert num_fractional == 0 || sharpness_bound >= lower_bound "$(num_fractional) == 0 || $(sharpness_bound) > $(lower_bound)"
+        end
+
+        lower_bound = max(strong_convexity_bound, sharpness_bound)
+    end
+    return lower_bound
+end
+
+"""
+Use strong convexity and/or sharpness to potentially remove one of the children nodes.
+If both sharpness and strong convexity parameters are provided, strong convexity is preferred.
+"""
+function prune_children(tree, node, lower_bound_base, x, vidx)
+    prune_left = false
+    prune_right = false
+
+    μ = tree.root.options[:strong_convexity]
+    M = tree.root.options[:sharpness_constant]
+    θ = tree.root.options[:sharpness_exponent]
+
+    if μ > 0 || (M > 0 && θ != Inf)
+        bound_improvement = 0.0
+        for j in tree.root.problem.integer_variables
+            if vidx == j
+                continue
+            end
+            bound_improvement += min((x[j] - floor(x[j]))^2, (ceil(x[j]) - x[j])^2)
+        end
+        new_bound_left = bound_improvement + (x[vidx] - floor(x[vidx]))^2
+        new_bound_right = bound_improvement + (ceil(x[vidx]) - x[vidx])^2
+
+        # strong convexity
+        if μ > 0
+            new_bound_left = lower_bound_base + μ / 2 * new_bound_left
+            new_bound_right = lower_bound_base + μ / 2 * new_bound_right
+
+            if new_bound_left > tree.incumbent
+                @debug "prune left, from $(node.lb) -> $new_bound_left, ub $(tree.incumbent), lb $(node.lb)"
+                prune_left = true
+            end
+            if new_bound_right > tree.incumbent
+                @debug "prune right, from $(node.lb) -> $new_bound_right, ub $(tree.incumbent), lb $(node.lb)"
+                prune_right = true
+            end
+            @assert !(
+                (new_bound_left > tree.incumbent + tree.root.options[:dual_gap]) &&
+                (new_bound_right > tree.incumbent + tree.root.options[:dual_gap])
+            ) 
+        # sharpness
+        elseif M > 0 && θ != Inf
+            fx = tree.root.problem.f(x)
+
+            new_bound_left = M^(- 1 / θ) * 1 / 2 * (sqrt(new_bound_left) - M / 2 * node.dual_gap^θ)^(1 / θ) + fx - node.dual_gap
+            new_bound_right = M^(- 1 / θ) * 1 / 2 * (sqrt(new_bound_right) - M / 2 * node.dual_gap^θ)^(1 / θ) + fx - node.dual_gap
+
+            if new_bound_left > tree.incumbent
+                @debug "prune left, from $(node.lb) -> $new_bound_left, ub $(tree.incumbent), lb $(node.lb)"
+                prune_left = true
+            end
+            if new_bound_right > tree.incumbent
+                @debug "prune right, from $(node.lb) -> $new_bound_right, ub $(tree.incumbent), lb $(node.lb)"
+                prune_right = true
+            end
+            @assert !(
+                (new_bound_left > tree.incumbent + tree.root.options[:dual_gap]) &&
+                (new_bound_right > tree.incumbent + tree.root.options[:dual_gap])
+            ) 
+        end
+    end
+   
+    # If both nodes are pruned, when one of them has to be equal to the incumbent.
+    # Thus, we have proof of optimality by strong convexity.
+    if prune_left && prune_right
+        tree.lb = min(new_bound_left, new_bound_right)
+    end
+
+    return prune_left, prune_right
+end
diff --git a/test/interface_test.jl b/test/interface_test.jl
index 572f0b557..cbd812c13 100644
--- a/test/interface_test.jl
+++ b/test/interface_test.jl
@@ -117,6 +117,63 @@ end
     end
 end
 
+@testset "Normbox - strong convexity and sharpness" begin
+    function f(x)
+        return 0.5 * sum((x[i] - diffi[i])^2 for i in eachindex(x))
+    end
+    function grad!(storage, x)
+        @. storage = x - diffi
+    end
+
+    @testset "Strong convexity" begin
+        int_vars = collect(1:n)
+        lbs = zeros(n)
+        ubs = ones(n)
+
+        sblmo = Boscia.CubeSimpleBLMO(lbs, ubs, int_vars)
+        μ = 1.0
+
+        x, _, result =
+            Boscia.solve(f, grad!, sblmo, lbs[int_vars], ubs[int_vars], int_vars, n, strong_convexity=μ)
+
+        @test x == round.(diffi)
+        @test isapprox(f(x), f(result[:raw_solution]), atol=1e-6, rtol=1e-3)
+    end
+
+    @testset "Sharpness" begin
+        int_vars = collect(1:n)
+        lbs = zeros(n)
+        ubs = ones(n)
+
+        sblmo = Boscia.CubeSimpleBLMO(lbs, ubs, int_vars)
+        θ = 1/2
+        M = 2.0
+
+        x, _, result =
+            Boscia.solve(f, grad!, sblmo, lbs[int_vars], ubs[int_vars], int_vars, n, sharpness_constant=M, sharpness_exponent=θ)
+
+        @test x == round.(diffi)
+        @test isapprox(f(x), f(result[:raw_solution]), atol=1e-6, rtol=1e-3)
+    end
+
+    @testset "Strong convexity and sharpness" begin
+        int_vars = collect(1:n)
+        lbs = zeros(n)
+        ubs = ones(n)
+
+        sblmo = Boscia.CubeSimpleBLMO(lbs, ubs, int_vars)
+        μ = 1.0
+        θ = 1/2
+        M = 2.0
+
+        x, _, result =
+            Boscia.solve(f, grad!, sblmo, lbs[int_vars], ubs[int_vars], int_vars, n, strong_convexity=μ, sharpness_constant=M, sharpness_exponent=θ)
+
+        @test x == round.(diffi)
+        @test isapprox(f(x), f(result[:raw_solution]), atol=1e-6, rtol=1e-3)
+    end
+end
+
 @testset "Start with Active Set" begin
 
     function f(x)