using LinearAlgebra, Statistics, StatsBase, Distributions, JuMP, DelimitedFiles
using Gurobi, Random, Ipopt, Mosek, MosekTools, Convex, SCS, CSV, DataFrames

function calculate_analytic_center_qp_exp(α,β,C,d)
    n = size(C)[2]
    m1 = Model(Ipopt.Optimizer)
    @variable(m1, x[1:n])
    @variable(m1, t)

    @constraint(m1, C*x .<= d)
    @NLconstraint(m1, [i in 1:n], exp(x[i]) <= β)
    @NLconstraint(m1, t == sum(exp(-x[i]) for i in 1:n))
    @NLconstraint(m1, log(t) <= α)

    @NLobjective(m1, Max, sum(log(d[i] - sum(C[i,j]*x[j] for j in 1:n)) for i in 1:size(C)[1]) +
                          sum(log(β - exp(x[i])) for i in 1:n) +
                          log(α - log(t)))
    optimize!(m1)
    return JuMP.value.(x)
end


function generate_hyperplane_mip_qp_exp(X,α,β,C,d,δ1,use_ac)
    M = 1e5
    n1, L1 = size(X)
    if mod(L1,2) == 0
        q = Int64(L1/2)
    else
        q = Int64((L1+1)/2)
    end
    m1 = Model(Gurobi.Optimizer)
    @variable(m1, f[1:n1])
    @variable(m1, l)
    @variable(m1, z[1:L1], Bin)
    @constraint(m1, sum(f[i] for i in 1:n1) >= δ1)
    if use_ac
        xac = calculate_analytic_center_qp_exp(α,β,C,d)
        @constraint(m1, f'*xac == l)
    end
    @constraint(m1, [i in 1:L1], f'*X[:,i] <= l + M*(1-z[i]))
    @constraint(m1, [i in 1:L1], f'*X[:,i] >= l - M*z[i])
    @constraint(m1, sum(z[i] for i in 1:L1) == q)
    @objective(m1, Min, f'*f)
    optimize!(m1)
    if termination_status(m1) == MOI.OPTIMAL
        return false, JuMP.value.(f), JuMP.value.(l)
    else
        return true, zeros(n), 0.0, zeros(n)
    end
end



function calculate_ub_qp_exp(Y,P,q,r,α,β,C,d)
    best_ub = 1e6
    if length(Y) > size(q)[1]
        for i in 1:size(Y)[2]
            x = Y[:,i]
            cur_ub = x'*P*x + q'*x + r
            # cur_ub = qp_exp_ub_ipopt(P,q,r,α,β,C,d,x)
            if cur_ub < best_ub
                best_ub = cur_ub
            end
        end
    else
        x = Y
        cur_ub = x'*P*x + q'*x + r
        # cur_ub = qp_exp_ub_ipopt(P,q,r,α,β,C,d,x)
        if cur_ub < best_ub
            best_ub = cur_ub
        end
    end
    return best_ub
end


function solve_nonconv_quad_exp_new(P,q,r,α,β,C,d,use_lmi,use_new)
    n = size(q)[1]

    m1 = Model(Mosek.Optimizer)
    @variable(m1, x[1:n])
    @variable(m1, z[1:n])
    @variable(m1, t[1:n])

    @variable(m1, X[1:n,1:n], Symmetric)
    @variable(m1, Z[1:n,1:n], Symmetric)
    @variable(m1, T[1:n,1:n], Symmetric)
    @variable(m1, V[1:n,1:n])
    @variable(m1, Q[1:n,1:n])
    @variable(m1, W[1:n,1:n])

    @constraint(m1, [i in 1:n], X[i,i] >= 0)
    @constraint(m1, [i in 1:n], Z[i,i] >= 0)
    @constraint(m1, [i in 1:n], T[i,i] >= 0)

    @constraint(m1, sum(z[i] for i in 1:n) <= 1)
    @constraint(m1, [i in 1:n], [- x[i] - α, 1, z[i]] in MOI.ExponentialCone())
    @constraint(m1, sum(t[i] for i in 1:n) <= β)
    @constraint(m1, [i in 1:n], [x[i], 1, t[i]] in MOI.ExponentialCone())

    @constraint(m1, [i in 1:n], [- x[i] - α + sum(V[i,j] for j in 1:n) +
                                 α*sum(z[j] for j in 1:n),
                                 1 - sum(z[j] for j in 1:n),
                                 z[i] - sum(Z[i,j] for j in 1:n)] in MOI.ExponentialCone())

    @constraint(m1, β - sum(t[j] for j in 1:n) - β*sum(z[j] for j in 1:n) + sum(Q[i,j] for i in 1:n for j in 1:n) >= 0)

    @constraint(m1, [i in 1:n], [x[i] - sum(V[i,j] for j in 1:n),
                                 1 - sum(z[j] for j in 1:n),
                                 t[i] - sum(Q[j,i] for j in 1:n)] in MOI.ExponentialCone())

    @constraint(m1, sum(Z[i,j] for i in 1:n for j in 1:n) - 2*sum(z[i] for i in 1:n) + 1 >= 0)

    @constraint(m1, [i in 1:n], [-β*x[i] - α*β + sum(W[i,j] for j in 1:n) + α*sum(t[j] for j in 1:n),
                                 β - sum(t[j] for j in 1:n),
                                 β*z[i] - sum(Q[i,j] for j in 1:n)] in MOI.ExponentialCone())

    @constraint(m1, [i in 1:n], [β*x[i] - sum(W[i,j] for j in 1:n),
                                 β - sum(t[j] for j in 1:n),
                                 β*t[i] - sum(T[i,j] for j in 1:n)] in MOI.ExponentialCone())
    @constraint(m1, sum(T[i,j] for i in 1:n for j in 1:n) - 2*sum(t[j] for j in 1:n) + β^2 >= 0)
    @variable(m1, Q_aux_1[1:n,1:n])
    @variable(m1, V_aux_1[1:n,1:n])
    @variable(m1, z_aux_1[1:n,1:n])
    @constraint(m1, Q_aux_1 .== Q)
    @constraint(m1, V_aux_1 .== V)
    @constraint(m1, [i in 1:n], z_aux_1[:,i] .== z)
    @constraint(m1, [i in 1:n, j in 1:n], [V_aux_1[i,j], z_aux_1[j,i], Q_aux_1[j,i]] in MOI.ExponentialCone())

    @constraint(m1, [i in 1:n, j in 1:n], [- V[i,j] - α*z[j], z[j], Z[i,j]] in MOI.ExponentialCone())
    @constraint(m1, [i in 1:n, j in i:n], [- x[i] - x[j] - 2*α, 1, Z[i,j]] in MOI.ExponentialCone())
    @constraint(m1, [i in 1:n, j in 1:n], [- x[i] - α + x[j], 1, Q[i,j]] in MOI.ExponentialCone())
    @constraint(m1, [i in 1:n, j in 1:n], [- W[i,j] - α*t[j], t[j], Q[i,j]] in MOI.ExponentialCone())
    @variable(m1, W_aux_1[1:n,1:n])
    @variable(m1, T_aux_1[1:n,1:n])
    @variable(m1, t_aux_1[1:n,1:n])
    @constraint(m1, W_aux_1 .== W)
    @constraint(m1, T_aux_1 .== T)
    @constraint(m1, [i in 1:n], t_aux_1[:,i] .== t)
    @constraint(m1, [i in 1:n, j in 1:n], [W_aux_1[i,j], t_aux_1[j,i], T_aux_1[i,j]] in MOI.ExponentialCone())
    @constraint(m1, [i in 1:n, j in i:n], [x[i] + x[j], 1, T[i,j]] in MOI.ExponentialCone())

    # Exponential Cone Decomposition
    @constraint(m1, [i in 1:n], - x[i] - α + 1 <= z[i])
    @constraint(m1, [i in 1:n], x[i] + 1 <= t[i])

    @constraint(m1, [i in 1:n], z[i] - sum(Z[i,j] for j in 1:n) + x[i] + α - 1 -
                                sum(V[i,j] for j in 1:n) - (α-1)*sum(z[j] for j in 1:n) >= 0)

    @constraint(m1, [i in 1:n, j in 1:n], [- V[i,j] - X[i,j] - α*x[i] + x[i] - α*(z[j]+x[j]+α-1),
                                            z[j] + x[j] + α - 1,
                                            Z[i,j] + V[j,i] + (α-1)*z[i]] in MOI.ExponentialCone())

    @constraint(m1, [i in 1:n, j in 1:n], [V[i,j] + X[i,j] + (α-1)*x[i],
                                           z[j] + x[j] + α - 1,
                                           Q[j,i] + W[j,i] + (α-1)*t[i] ] in MOI.ExponentialCone())

    @constraint(m1, [i in 1:n, j in 1:n], (α-1)*(z[j] + x[j] + α - 1) + Z[i,j] + V[j,i] + (α-1)*z[i] +
                                           V[i,j] + X[i,j] + (α-1)*x[i] >= 0)

    @constraint(m1, [i in 1:n], t[i] - sum(Q[j,i] for j in 1:n) - x[i] +
                                sum(V[i,j] for j in 1:n) - 1 + sum(z[j] for j in 1:n) >= 0)
    @constraint(m1, [i in 1:n, j in 1:n], [- W[i,j] - α*t[j] + X[i,j] + α*(x[j] + 1) + x[i],
                                             t[j] - x[j] - 1,
                                             Q[i,j] - V[j,i] - z[i]] in MOI.ExponentialCone())
    @constraint(m1, [i in 1:n, j in 1:n], [W[i,j] - X[i,j] - x[i],
                                           t[j] - x[j] - 1,
                                           T[i,j]-Q[j,i]- t[i]] in MOI.ExponentialCone())

    @constraint(m1, [i in 1:n, j in 1:n], X[i,j]+T[i,j]-W[i,j]-W[j,i]+x[i]+x[j]-t[i]-t[j]+1 >= 0)
    @constraint(m1, [i in 1:n, j in 1:n], Q[j,i]+W[j,i]-X[i,j]-V[i,j]+α*t[i]-t[i]-α*x[i]+x[i]-z[j]-x[j]-α+1 >= 0)
    @constraint(m1, [j in 1:n], β*z[j] - sum(Q[j,i] for i in 1:n) + β*x[j] - sum(W[j,i] for i in 1:n) -
                                α*sum(t[i] for i in 1:n) + α*β - β + sum(t[i] for i in 1:n) >= 0)

    @constraint(m1, [j in 1:n], β*t[j] - sum(T[i,j] for i in 1:n) - β*x[j] + sum(W[j,i] for i in 1:n) -
                                β + sum(t[i] for i in 1:n) >= 0)
    if length(C) > 0
        L = size(C)[1]
        @constraint(m1, C*x .<= d)
        @constraint(m1, d*x'*C' .+ C*x*d' .<= d*d' .+ C*X*C')
        @constraint(m1, C*x .- C*sum(V[:,j] for j in 1:n) .<= d*(1-sum(z[j] for j in 1:n)))
        @constraint(m1, β*C*x .- C*(sum(W[:,j] for j in 1:n)) .<= (β-sum(t[j] for j in 1:n))*d)
        @constraint(m1, [i in 1:n, j in 1:L], [- d[j]*x[i] - d[j]*α + C[j,:]'*X[:,i] + α*C[j,:]'*x,
                                                d[j]-C[j,:]'*x, d[j]*z[i] - C[j,:]'*V[:,i]] in MOI.ExponentialCone())
        @constraint(m1, [i in 1:n, j in 1:L], [d[j]*x[i]-C[j,:]'*X[:,i], d[j]-C[j,:]'*x, d[j]*t[i] - C[j,:]'*W[:,i]] in MOI.ExponentialCone())

    end

    if use_lmi
        @constraint(m1, [X V W x; V' Z Q z; W' Q' T t; x' z' t' 1] in PSDCone())
    end
    @objective(m1, Min, tr(P*X) + q'*x + r)
    optimize!(m1)
    if termination_status(m1) == MOI.OPTIMAL #|| termination_status(m1) == MOI.SLOW_PROGRESS
        return true, JuMP.value.(x), JuMP.value.(z), JuMP.value.(X), JuMP.value.(V), objective_value(m1)
    else
        return false, zeros(n), zeros(n), zeros(n,n), zeros(n,n), 1e6
    end
end

function is_feasible_qp_exp_new(y,α,β)
    n = size(y)[1]
    cnt_1 = (log(sum(exp(-y[i]) for i in 1:n)) <= α)*1
    cnt_2 = (sum(exp(y[i]) for i in 1:n) <= β)*1
    return cnt_1 + cnt_2 == 2
end

function calculate_cand_vectors_qp_exp_new(x,z,X,V,α,β)
    Y = [x]
    for j in 1:size(x)[1]
        if x[j] == 0
            push!(Y,x)
        else
            if is_feasible_qp_exp_new(X[:,j]/x[j],α,β)
                push!(Y,X[:,j]/x[j])
            end
        end
    end

    for j in 1:size(z)[1]
        if z[j] == 0
            push!(Y,x)
        else
            if is_feasible_qp_exp_new(V[:,j]/z[j],α,β)
                push!(Y,V[:,j]/z[j])
            end
        end
    end

    return Y
end

function qp_exp_ub_ipopt(P,q,r,α,β,C,d,x0)
    n = size(q)[1]
    m1 = Model(Ipopt.Optimizer)
    @variable(m1, x[1:n])
    @variable(m1, t)

    for j in 1:n
        JuMP.set_start_value(x[j], x0[j])
    end
    if length(C) > 0
        @constraint(m1, C*x .<= d)
    end

    @NLconstraint(m1, [i in 1:n], exp(x[i]) <= β)
    @NLconstraint(m1, t == sum(exp(-x[i]) for i in 1:n))
    @NLconstraint(m1, log(t) <= α)

    @NLobjective(m1, Min, sum(P[i,j]*x[i]*x[j] for i in 1:n for j in 1:n) +
                          sum(q[j]*x[j] for j in 1:n) +r)
    optimize!(m1)
    return objective_value(m1)
end

function calculate_different_vectors(X,x)
    Y = x
    for i in 1:length(X)
        if norm(X[i] .- x) > 1e-3
            Y = hcat(Y,X[i])
        end
    end
    return Y
end


function rpt_branch_and_split_qp_exp_new(P,q,r,α,β,δ,δ1,use_lmi,use_new)
    C, d = [], []
    gen_hyper = 0

    # Root Node
    res_root, xlb, zlb, Xlb, Vlb, lb_root = solve_nonconv_quad_exp_new(P,q,r,α,β,C,d,use_lmi,use_new)
    ub_cand_list = calculate_cand_vectors_qp_exp_new(xlb,zlb,Xlb,Vlb,α,β)
    Y_split = calculate_different_vectors(ub_cand_list,xlb)
    ub_root = calculate_ub_qp_exp(Y_split,P,q,r,α,β,C,d)

    UB, LB, opt_sol, opt_val = ub_root, lb_root, xlb, ub_root
    # unvisited_nodes = Queue{Vector{}}()
    nodes_list = []
    ub_list, lb_list, hyp_list_f, hyp_list_l = [], [], [], []
    push!(ub_list, UB)
    push!(lb_list, LB)

    t0_1 = time_ns()
    total_time = 0.0

    while UB - LB > δ && total_time < 3600 && length(Y_split) > size(P)[1]
        # Split
        use_ac = false
        if gen_hyper < 1
            use_ac = false
        end
        res, fopt, lopt = generate_hyperplane_mip_qp_exp(Y_split,α,β,C,d,δ1,use_ac)
        if res
            L1 = size(Y_split)[2]
            q1 = Int64(round(L1/2))
            res_1, fopt, lopt = generate_hyperplane_mip_qp_exp(Y_split[:,1:q1],α,β,C,d,δ1,use_ac)
        end
        fr, lr, fl, ll  = fopt, lopt, -fopt, -lopt
        if gen_hyper < 1
            Cr, dr, Cl, dl = Matrix(fr'), [lr], Matrix(fl'), [ll]
        else
            Cr, dr, Cl, dl = vcat(C,fr'), vcat(d,lr), vcat(C,fl'), vcat(d,ll)
        end
        gen_hyper += 1
        push!(hyp_list_f, fopt)
        push!(hyp_list_l, lopt)


        # Right child
        res_right, xlb_right, zlb_right, Xlb_right, Vlb_right, lb_right = solve_nonconv_quad_exp_new(P,q,r,α,β,Cr,dr,use_lmi,use_new)
        if res_right
            cand_list_right = calculate_cand_vectors_qp_exp_new(xlb_right,zlb_right,Xlb_right,Vlb_right,α,β)
            Y_right = calculate_different_vectors(cand_list_right,xlb_right)
            ub_right = calculate_ub_qp_exp(Y_right,P,q,r,α,β,Cr,dr)
            if lb_right < UB
                push!(nodes_list,[ub_right,lb_right,xlb_right,Cr,dr,Y_right])
            end
        end
        # Left child
        res_left, xlb_left, zlb_left, Xlb_left, Vlb_left, lb_left = solve_nonconv_quad_exp_new(P,q,r,α,β,Cl,dl,use_lmi,use_new)
        if res_left
            cand_list_left = calculate_cand_vectors_qp_exp_new(xlb_left,zlb_left,Xlb_left,Vlb_left,α,β)
            Y_left = calculate_different_vectors(cand_list_left,xlb_left)
            ub_left = calculate_ub_qp_exp(Y_left,P,q,r,α,β,Cl,dl)
            if lb_left < UB
                push!(nodes_list,[ub_left,lb_left,xlb_left,Cl,dl,Y_left])
            end
        end

        if isempty(nodes_list)
            if res_right
                push!(ub_list, ub_right)
                push!(lb_list, lb_right)
            end
            if res_left
                push!(ub_list, ub_left)
                push!(lb_list, lb_left)
            end
            break
        else
            ind = argmin([nodes_list[i][2] for i in 1:length(nodes_list)])
            cur_node = nodes_list[ind]
            deleteat!(nodes_list, ind)
            cur_ub, cur_lb = cur_node[1], cur_node[2]
            C, d, Y_split = cur_node[4], cur_node[5], cur_node[6]
            LB = cur_lb
            if cur_ub < UB
                UB = cur_ub
                opt_sol, opt_val = cur_node[3], cur_node[1]
            end
            push!(ub_list, cur_ub)
            push!(lb_list, cur_lb)
        end
        t0_2 = time_ns()
        total_time = (t0_2-t0_1)*10^(-9)
    end
    return opt_sol, opt_val, gen_hyper, ub_list, lb_list, hyp_list_f, hyp_list_l
end

function matrix_string_to_float(A)
    m = size(A)[2]
    for j in 1:m
        if typeof(A[1,j]) != Float64
            A[:,j] = parse.(Float64, A[:,j])
        end
    end
    return A
end


δ = 1e-4
δ1 = 1
use_new = true



# Instance 1
n, L = 5, 10
id_mat = zeros(n,n) + I
P = -0.5*id_mat
q = -5*ones(n)
r = -12.5*n
α = 2
β = 20
C, d = [], []
use_lmi = false
res_root, xlb, zlb, Xlb, Vlb, lb = solve_nonconv_quad_exp_new(P,q,r,α,β,C,d,use_lmi,use_new)
ub_cand_list = calculate_cand_vectors_qp_exp_new(xlb,zlb,Xlb,Vlb,α,β)
Y_list = calculate_different_vectors(ub_cand_list,xlb)
ub = calculate_ub_qp_exp(Y_list,P,q,r,α,β,C,d)

println("------------------------------------------------")
println("Lower Bound:")
println(lb)
println("Upper Bound:")
println(ub)

use_lmi = false
t1 = time_ns()
x_opt, obj_opt, gen_hyper, ub_list, lb_list, hyp_list_f, hyp_list_l = rpt_branch_and_split_qp_exp_new(P,q,r,α,β,δ,δ1,use_lmi,use_new)
t2 = time_ns()
total_time = (t2-t1)*10^(-9)

use_lmi = true
t1 = time_ns()
x_opt, obj_opt_sdp, gen_hyper_sdp, ub_list, lb_list, hyp_list_f, hyp_list_l = rpt_branch_and_split_qp_exp_new(P,q,r,α,β,δ,δ1,use_lmi,use_new)
t2 = time_ns()
total_time_sdp = (t2-t1)*10^(-9)

println("---------------------------------------------------")
println(obj_opt)
println(obj_opt_sdp)

println("---------------------------------------------------")
println(total_time)
println(total_time_sdp)


# writedlm("QP_EXP_P1_P.csv", P, ',')
# writedlm("QP_EXP_P1_q.csv", q, ',')


# Instance 2
n, L = 5, 10
id_mat = zeros(n,n) + I
P = -0.5*id_mat
q = -7*ones(n)
r = -24.5*n
α = 2
β = 20
C, d = [], []
use_lmi = false
res_root, xlb, zlb, Xlb, Vlb, lb = solve_nonconv_quad_exp_new(P,q,r,α,β,C,d,use_lmi,use_new)
ub_cand_list = calculate_cand_vectors_qp_exp_new(xlb,zlb,Xlb,Vlb,α,β)
Y_list = calculate_different_vectors(ub_cand_list,xlb)
ub = calculate_ub_qp_exp(Y_list,P,q,r,α,β,C,d)

println("------------------------------------------------")
println("Lower Bound:")
println(lb)
println("Upper Bound:")
println(ub)

use_lmi = false
t1 = time_ns()
x_opt, obj_opt, gen_hyper, ub_list, lb_list, hyp_list_f, hyp_list_l = rpt_branch_and_split_qp_exp_new(P,q,r,α,β,δ,δ1,use_lmi,use_new)
t2 = time_ns()
total_time = (t2-t1)*10^(-9)

use_lmi = true
t1 = time_ns()
x_opt, obj_opt_sdp, gen_hyper_sdp, ub_list, lb_list, hyp_list_f, hyp_list_l = rpt_branch_and_split_qp_exp_new(P,q,r,α,β,δ,δ1,use_lmi,use_new)
t2 = time_ns()
total_time_sdp = (t2-t1)*10^(-9)

println("---------------------------------------------------")
println(obj_opt)
println(obj_opt_sdp)

println("---------------------------------------------------")
println(total_time)
println(total_time_sdp)


writedlm("QP_EXP_P2_P.csv", P, ',')
writedlm("QP_EXP_P2_q.csv", q, ',')

######################   Run Averages   ########################

# Instance 3
n = 10
P_mat = CSV.read("QP_EXP_P3_Pmat.csv", DataFrame, header = false)
q_mat = CSV.read("QP_EXP_P3_qmat.csv", DataFrame, header = false)
P_mat = convert(Matrix, P_mat[:,:])
q_mat = convert(Matrix, q_mat[:,:])
r = 0
α, β = 5, 15
use_lmi = true
times, obj_vals, tree_depths = [], [], []
upper_bound_list, lower_bound_list, hyp_list_f, hyp_list_l = [], [], [], []
lbs_root, ubs_root = [], []
i, cnt = 1, 1
while i < 10*n
    P = P_mat[:,i:i+n-1]
    q = q_mat[:,cnt]
    C, d = [], []

    t1 = time_ns()
    x_opt, obj_opt, gen_hyper, ub_list, lb_list, hyp_list_f, hyp_list_l = rpt_branch_and_split_qp_exp_new(P,q,r,α,β,δ,δ1,use_lmi,use_new)
    t2 = time_ns()
    total_time = (t2-t1)*10^(-9)
    push!(times,total_time)
    push!(obj_vals,obj_opt)
    push!(tree_depths,gen_hyper)
    push!(lower_bound_list, lb_list)
    push!(upper_bound_list, ub_list)
    i += n
    cnt += 1
    println("------------------------------------------")
    println(cnt)
end
println("----------------------------")
println("Optimal objective:")
println(mean(obj_vals))
println("----------------------------")
println("Total time:")
println(mean(times))
println("----------------------------")
println("Total depth:")
println(mean(tree_depths))


# Instance 4
n = 20
P_mat = CSV.read("QP_EXP_P4_Pmat.csv", DataFrame, header = false)
q_mat = CSV.read("QP_EXP_P4_qmat.csv", DataFrame, header = false)
P_mat = convert(Matrix, P_mat[:,:])
q_mat = convert(Matrix, q_mat[:,:])
α, β = 5, 50
r = 0
use_lmi = false
times, obj_vals, tree_depths = [], [], []
upper_bound_list, lower_bound_list, hyp_list_f, hyp_list_l = [], [], [], []
lbs_root, ubs_root = [], []
i, cnt = 1, 1
while i < 10*n
    P = P_mat[:,i:i+n-1]
    q = q_mat[:,cnt]
    C, d = [], []

    t1 = time_ns()
    x_opt, obj_opt, gen_hyper, ub_list, lb_list, hyp_list_f, hyp_list_l = rpt_branch_and_split_qp_exp_new(P,q,r,α,β,δ,δ1,use_lmi,use_new)
    t2 = time_ns()
    total_time = (t2-t1)*10^(-9)
    push!(times,total_time)
    push!(obj_vals,obj_opt)
    push!(tree_depths,gen_hyper)
    push!(lower_bound_list, lb_list)
    push!(upper_bound_list, ub_list)
    i += n
    cnt += 1
    println("------------------------------------------")
    println(cnt)
end
println("----------------------------")
println("Optimal objective:")
println(mean(obj_vals))
println("----------------------------")
println("Total time:")
println(mean(times))
println("----------------------------")
println("Total depth:")
println(mean(tree_depths))


# Instance 5
n = 50
P_mat = CSV.read("QP_EXP_P5_Pmat.csv", DataFrame, header = false)
q_mat = CSV.read("QP_EXP_P5_qmat.csv", DataFrame, header = false)
P_mat = convert(Matrix, P_mat[:,:])
q_mat = convert(Matrix, q_mat[:,:])
r = 0
α, β = 5, 100
use_lmi = false
times, obj_vals, tree_depths = [], [], []
upper_bound_list, lower_bound_list, hyp_list_f, hyp_list_l = [], [], [], []
lbs_root, ubs_root = [], []
i, cnt = 1, 1
while i < 10*n
    P = P_mat[:,i:i+n-1]
    q = q_mat[:,cnt]
    C, d = [], []
    res_root, xlb, zlb, Xlb, Vlb, lb_root_cur = solve_nonconv_quad_exp_new(P,q,r,α,β,C,d,use_lmi,use_new)
    ub_cand_list = calculate_cand_vectors_qp_exp_new(xlb,zlb,Xlb,Vlb,α,β)
    Y_list = calculate_different_vectors(ub_cand_list,xlb)
    ub_root_cur = calculate_ub_qp_exp(Y_list,P,q,r,α,β,C,d)
    push!(lbs_root, lb_root_cur)
    push!(ubs_root, ub_root_cur)

    t1 = time_ns()
    x_opt, obj_opt, gen_hyper, ub_list, lb_list, hyp_list_f, hyp_list_l = rpt_branch_and_split_qp_exp_new(P,q,r,α,β,δ,δ1,use_lmi,use_new)
    t2 = time_ns()
    total_time = (t2-t1)*10^(-9)
    push!(times,total_time)
    push!(obj_vals,obj_opt)
    push!(tree_depths,gen_hyper)
    push!(lower_bound_list, lb_list)
    push!(upper_bound_list, ub_list)
    i += n
    cnt += 1
    println("------------------------------------------")
    println(cnt)
end
println("----------------------------")
println("Optimal objective:")
println(mean(obj_vals))
println("----------------------------")
println("Total time:")
println(mean(times))
println("----------------------------")
println("Total depth:")
println(mean(tree_depths))

# Instance 6
n = 100
P_mat = CSV.read("QP_EXP_P6_Pmat.csv", DataFrame, header = false)
q_mat = CSV.read("QP_EXP_P6_qmat.csv", DataFrame, header = false)
P_mat = convert(Matrix, P_mat[:,:])
q_mat = convert(Matrix, q_mat[:,:])
P_mat = matrix_string_to_float(P_mat)
P_mat = Float64.(P_mat)
r = 0
α, β = 5, 500
use_lmi = false
times, obj_vals, tree_depths = [], [], []
upper_bound_list, lower_bound_list, hyp_list_f, hyp_list_l = [], [], [], []
lbs_root, ubs_root = [], []
i, cnt = 1, 1
while i < 10*n
    P = P_mat[:,i:i+n-1]
    q = q_mat[:,cnt]
    C, d = [], []
    res_root, xlb, zlb, Xlb, Vlb, lb_root_cur = solve_nonconv_quad_exp_new(P,q,r,α,β,C,d,use_lmi,use_new)
    ub_cand_list = calculate_cand_vectors_qp_exp_new(xlb,zlb,Xlb,Vlb,α,β)
    Y_list = calculate_different_vectors(ub_cand_list,xlb)
    ub_root_cur = calculate_ub_qp_exp(Y_list,P,q,r,α,β,C,d)
    push!(lbs_root, lb_root_cur)
    push!(ubs_root, ub_root_cur)

    t1 = time_ns()
    x_opt, obj_opt, gen_hyper, ub_list, lb_list, hyp_list_f, hyp_list_l = rpt_branch_and_split_qp_exp_new(P,q,r,α,β,δ,δ1,use_lmi,use_new)
    t2 = time_ns()
    total_time = (t2-t1)*10^(-9)
    push!(times,total_time)
    push!(obj_vals,obj_opt)
    push!(tree_depths,gen_hyper)
    push!(lower_bound_list, lb_list)
    push!(upper_bound_list, ub_list)
    i += n
    cnt += 1
    println("------------------------------------------")
    println(cnt)
end
println("----------------------------")
println("Optimal objective:")
println(mean(obj_vals))
println("----------------------------")
println("Total time:")
println(mean(times))
println("----------------------------")
println("Total depth:")
println(mean(tree_depths))
