From 28e7e93302519f32688d688668139c16dd6a541a Mon Sep 17 00:00:00 2001 From: dhendryc <92737336+dhendryc@users.noreply.github.com> Date: Fri, 10 May 2024 15:28:59 +0200 Subject: [PATCH] Changes from SEA camera ready version (#2) --- ODWB/src/opt_design_boscia.jl | 117 +++++++++------ ODWB/src/opt_design_custom_BB.jl | 249 ++++++++++++++++++------------- ODWB/src/opt_design_scip.jl | 28 ++-- ODWB/src/utilities.jl | 164 ++++++++++++-------- ODWB/test/test_derivatives.jl | 21 ++- 5 files changed, 351 insertions(+), 228 deletions(-) diff --git a/ODWB/src/opt_design_boscia.jl b/ODWB/src/opt_design_boscia.jl index 3045601..069d529 100644 --- a/ODWB/src/opt_design_boscia.jl +++ b/ODWB/src/opt_design_boscia.jl @@ -41,12 +41,31 @@ # m - number of possible experiments # A = [v_1^T,.., v_m^T], so the rows of A correspond to the different experiments -function solve_opt(seed, m, n, time_limit, criterion, corr; full_callback=true, write = true, verbose = true, use_scip=false, do_strong_branching=false, use_shadow_set=false, lazy_tolerance=2.0, use_heuristics=false, use_tightening=false, long_runs=false, options_run=false) +function solve_opt(seed, m, n, time_limit, criterion, corr; full_callback=true, p=0, write = true, verbose = true, use_scip=false, do_strong_branching=false, use_shadow_set=false, lazy_tolerance=2.0, use_heuristics=false, use_tightening=false, long_runs=false, options_run=false, fw_verbose=false, specific_seed=false) if criterion == "AF" || criterion == "DF" - A, C, N, ub = build_data(seed, m, n, true, corr) + A, C, N, ub, _ = build_data(seed, m, n, true, corr; scaling_C=long_runs && criterion != "AF" && criterion != "DF") else - A, _, N, ub = build_data(seed, m, n, false, corr) + A, _, N, ub, _ = build_data(seed, m, n, false, corr; scaling_C=long_runs) + end + + # parameter tunning + if !options_run + use_heuristics = true + if !(criterion in ["D","DF"]) + use_shadow_set = true + elseif !(criterion in ["A","AF"]) + lazy_tolerance = 1.5 + end + end + + if long_runs + use_heuristics = true + if criterion in ["A","AF"] + use_shadow_set = true + elseif criterion in ["D","DF"] + lazy_tolerance = 1.5 + end end if use_scip @@ -74,39 +93,46 @@ function solve_opt(seed, m, n, time_limit, criterion, corr; full_callback=true, result = 0.0 domain_oracle = build_domain_oracle(A, n) + println("build function") if criterion == "A" - f, grad! = build_a_criterion(A, false, μ=1e-4, build_safe=true) + f, grad! = build_a_criterion(A, false, μ=1e-4, build_safe=true, long_run=long_runs) + #f, grad! = build_general_trace(A, -1, false, build_safe=true) # Log(Tr(X^{-1})) elseif criterion == "AF" - f, grad! = build_a_criterion(A, true, C=C) + f, grad! = build_a_criterion(A, true, C=C, long_run=long_runs) + #f, grad! = build_general_trace(A, -1, true, C=C) # Log(Tr(X^{-1})) elseif criterion =="D" - f, grad! = build_d_criterion(A, false, μ=1e-4, build_safe=true) + f, grad! = build_d_criterion(A, false, μ=1e-4, build_safe=true, long_run=long_runs) elseif criterion == "DF" - f, grad! = build_d_criterion(A, true, C=C) + f, grad! = build_d_criterion(A, true, C=C, long_run=long_runs) else error("Invalid criterion!") end - if criterion == "AF" || criterion == "DF" + if criterion in ["AF","DF"] direction = collect(1.0:m) - x0 = compute_extreme_point(lmo, direction) - active_set= FrankWolfe.ActiveSet([(1.0, x0)]) + x0 = compute_extreme_point(lmo, direction) + active_set= FrankWolfe.ActiveSet([(1.0, x0)]) + # set same incumbent as for Co-BnB + z = greedy_incumbent_fusion(A,m,n,N,ub) + # Precompile - x, _, result = Boscia.solve(f, grad!, lmo; verbose=false, time_limit=10, active_set=active_set, branching_strategy=branching_strategy, use_shadow_set=use_shadow_set, dual_tightening=use_tightening, global_dual_tightening=use_tightening, lazy_tolerance=lazy_tolerance, custom_heuristics=[heu]) + x, _, result = Boscia.solve(f, grad!, lmo; verbose=false, time_limit=10, active_set=active_set, branching_strategy=branching_strategy, use_shadow_set=use_shadow_set, dual_tightening=use_tightening, global_dual_tightening=use_tightening, lazy_tolerance=lazy_tolerance, custom_heuristics=[heu], start_solution=z) - x, _, result = Boscia.solve(f, grad!, lmo; verbose=verbose, time_limit=time_limit, active_set=active_set, branching_strategy=branching_strategy, use_shadow_set=use_shadow_set, dual_tightening=use_tightening, global_dual_tightening=use_tightening, lazy_tolerance=lazy_tolerance, custom_heuristics=[heu] ) + # Actual Run + x, _, result = Boscia.solve(f, grad!, lmo; verbose=verbose, time_limit=time_limit, active_set=active_set, branching_strategy=branching_strategy, use_shadow_set=use_shadow_set, dual_tightening=use_tightening, global_dual_tightening=use_tightening, lazy_tolerance=lazy_tolerance, custom_heuristics=[heu], fw_verbose=fw_verbose, start_solution=z) else - # Precompile _, active_set, S = build_start_point2(A, m, n, N, ub) z = greedy_incumbent(A, m, n, N, ub) - x, _, result = Boscia.solve(f, grad!, lmo; verbose=false, time_limit=10, active_set=active_set, domain_oracle=domain_oracle, start_solution=z, dual_tightening=use_tightening, global_dual_tightening=use_tightening, lazy_tolerance=lazy_tolerance, branching_strategy=branching_strategy, use_shadow_set=use_shadow_set, custom_heuristics=[heu]) #line_search=StepSizeRule, + + # Precompile + x, _, result = Boscia.solve(f, grad!, lmo; verbose=false, time_limit=10, active_set=active_set, domain_oracle=domain_oracle, start_solution=z, dual_tightening=use_tightening, global_dual_tightening=use_tightening, lazy_tolerance=lazy_tolerance, branching_strategy=branching_strategy, use_shadow_set=use_shadow_set, custom_heuristics=[heu]) - # real run - # Find a good start point _, active_set, S = build_start_point2(A, m, n, N, ub) - # initial upper bound z = greedy_incumbent(A, m, n, N, ub) - x, _, result = Boscia.solve(f, grad!, lmo; verbose=verbose, time_limit=time_limit, active_set=active_set, domain_oracle=domain_oracle, start_solution=z, dual_tightening=use_tightening, global_dual_tightening=use_tightening, lazy_tolerance=lazy_tolerance, branching_strategy=branching_strategy, use_shadow_set=use_shadow_set, custom_heuristics=[heu]) #line_search=StepSizeRule, + + # Actual run + x, _, result = Boscia.solve(f, grad!, lmo; verbose=verbose, time_limit=time_limit, active_set=active_set, domain_oracle=domain_oracle, start_solution=z, dual_tightening=use_tightening, global_dual_tightening=use_tightening, lazy_tolerance=lazy_tolerance, branching_strategy=branching_strategy, use_shadow_set=use_shadow_set, custom_heuristics=[heu], fw_verbose=fw_verbose) end total_time_in_sec=result[:total_time_in_sec] @@ -114,6 +140,9 @@ function solve_opt(seed, m, n, time_limit, criterion, corr; full_callback=true, if occursin("Optimal", result[:status]) status = "OPTIMAL" end + if occursin("Time", result[:status]) + status = "TIME_LIMIT" + end if full_callback lb_list = result[:list_lb] ub_list = result[:list_ub] @@ -124,7 +153,23 @@ function solve_opt(seed, m, n, time_limit, criterion, corr; full_callback=true, end if write - if options_run + if long_runs + # CSV file for the full callback + type = corr ? "correlated" : "independent" + df_full_cb = DataFrame(seed=seed, numberOfExperiments=m, numberOfParameters=n, N=N, time=time_list, lowerBound=lb_list, upperBound =ub_list, termination=status, LMOcalls =list_lmo_calls, localTighteings=list_local_tightening, globalTightenings=list_global_tightening) + file_name_full_cb = "/home/htc/dhendryc/research_projects/MasterThesis/optDesign/csv/full_runs_boscia/long_runs/boscia_"* criterion * "_optimality_" * type *"_" * string(m) * "-" * string(n) * "_" * string(seed) *".csv" + CSV.write(file_name_full_cb, df_full_cb, append=false) + + # CSV file for the results of all instances. + scaled_solution = result[:primal_objective]*m + df = DataFrame(seed=seed, numberOfExperiments=m, numberOfParameters=n, N=N, time=total_time_in_sec, solution=result[:primal_objective], scaled_solution=scaled_solution, dual_gap = result[:dual_gap], rel_dual_gap=result[:rel_dual_gap], ncalls=result[:lmo_calls], num_nodes=result[:number_nodes],termination=status) + file_name = "/home/htc/dhendryc/research_projects/MasterThesis/optDesign/csv/Boscia/long_runs/boscia_" * criterion * "_" * string(m) * "_" * type *"_optimality.csv" + if !isfile(file_name) + CSV.write(file_name, df, append=true, writeheader=true, delim=";") + else + CSV.write(file_name, df, append=true, delim=";") + end + elseif options_run # CSV file for the full callback type = corr ? "correlated" : "independent" option = if use_scip @@ -144,29 +189,13 @@ function solve_opt(seed, m, n, time_limit, criterion, corr; full_callback=true, end df_full_cb = DataFrame(seed=seed, numberOfExperiments=m, numberOfParameters=n, N=N, time=time_list, lowerBound=lb_list, upperBound =ub_list, termination=status, LMOcalls =list_lmo_calls, localTighteings=list_local_tightening, globalTightenings=list_global_tightening) - file_name_full_cb = joinpath(@__DIR__, "../csv/full_runs_boscia/" * option * "/boscia_"* criterion * "_optimality_" * type *"_" * string(m) * "-" * string(n) * "_" * string(seed) *".csv") - CSV.write(file_name_full_cb, df_full_cb, append=false) + file_name_full_cb = "/home/htc/dhendryc/research_projects/MasterThesis/optDesign/csv/full_runs_boscia/" * option * "/boscia_"* criterion * "_optimality_" * type *"_" * string(m) * "-" * string(n) * "_" * string(seed) *".csv" + CSV.write(file_name_full_cb, df_full_cb, append=false) # CSV file for the results of all instances. scaled_solution = result[:primal_objective]*m df = DataFrame(seed=seed, numberOfExperiments=m, numberOfParameters=n, N=N, time=total_time_in_sec, solution=result[:primal_objective], scaled_solution=scaled_solution, dual_gap = result[:dual_gap], rel_dual_gap=result[:rel_dual_gap], ncalls=result[:lmo_calls], num_nodes=result[:number_nodes],termination=status) - file_name = joinpath(@__DIR__, "../csv/Boscia/" * option * "/boscia_" * criterion * "_" * string(m) * "_" * type *"_optimality.csv") - if !isfile(file_name) - CSV.write(file_name, df, append=true, writeheader=true, delim=";") - else - CSV.write(file_name, df, append=true, delim=";") - end - elseif long_runs - # CSV file for the full callback - type = corr ? "correlated" : "independent" - df_full_cb = DataFrame(seed=seed, numberOfExperiments=m, numberOfParameters=n, N=N, time=time_list, lowerBound=lb_list, upperBound =ub_list, termination=status, LMOcalls =list_lmo_calls, localTighteings=list_local_tightening, globalTightenings=list_global_tightening) - file_name_full_cb = joinpath(@__DIR__, "../csv/full_runs_boscia/long_runs/boscia_"* criterion * "_optimality_" * type *"_" * string(m) * "-" * string(n) * "_" * string(seed) *".csv") - CSV.write(file_name_full_cb, df_full_cb, append=false) - - # CSV file for the results of all instances. - scaled_solution = result[:primal_objective]*m - df = DataFrame(seed=seed, numberOfExperiments=m, numberOfParameters=n, N=N, time=total_time_in_sec, solution=result[:primal_objective], scaled_solution=scaled_solution, dual_gap = result[:dual_gap], rel_dual_gap=result[:rel_dual_gap], ncalls=result[:lmo_calls], num_nodes=result[:number_nodes],termination=status) - file_name = joinpath(@__DIR__, "../csv/Boscia/long_runs/boscia_" * criterion * "_" * string(m) * "_" * type *"_optimality.csv") + file_name = "/home/htc/dhendryc/research_projects/MasterThesis/optDesign/csv/Boscia/" * option * "/boscia_" * criterion * "_" * string(m) * "_" * type *"_optimality.csv" if !isfile(file_name) CSV.write(file_name, df, append=true, writeheader=true, delim=";") else @@ -176,13 +205,17 @@ function solve_opt(seed, m, n, time_limit, criterion, corr; full_callback=true, # CSV file for the full callback type = corr ? "correlated" : "independent" df_full_cb = DataFrame(seed=seed, numberOfExperiments=m, numberOfParameters=n, N=N, time=time_list, lowerBound=lb_list, upperBound =ub_list, termination=status, LMOcalls =list_lmo_calls, localTighteings=list_local_tightening, globalTightenings=list_global_tightening) - file_name_full_cb = joinpath(@__DIR__, "../csv/full_runs_boscia/boscia_"* criterion * "_optimality_" * type *"_" * string(m) * "-" * string(n) * "_" * string(seed) *".csv") - CSV.write(file_name_full_cb, df_full_cb, append=false) + file_name_full_cb = "/home/htc/dhendryc/research_projects/MasterThesis/optDesign/csv/full_runs_boscia/boscia_"* criterion * "_optimality_" * type *"_" * string(m) * "-" * string(n) * "_" * string(seed) *".csv" + CSV.write(file_name_full_cb, df_full_cb, append=false) # CSV file for the results of all instances. scaled_solution = result[:primal_objective]*m df = DataFrame(seed=seed, numberOfExperiments=m, numberOfParameters=n, N=N, time=total_time_in_sec, solution=result[:primal_objective], scaled_solution=scaled_solution, dual_gap = result[:dual_gap], rel_dual_gap=result[:rel_dual_gap], ncalls=result[:lmo_calls], num_nodes=result[:number_nodes],termination=status) - file_name = joinpath(@__DIR__, "../csv/Boscia/boscia_" * criterion * "_" * string(m) * "_" * type *"_optimality.csv") + if specific_seed + file_name = "/home/htc/dhendryc/research_projects/MasterThesis/optDesign/csv/Boscia/boscia_" * criterion * "_" * string(m) * "_" * string(n) * "_" * type * "_" * string(seed) * "_optimality.csv" + else + file_name = "/home/htc/dhendryc/research_projects/MasterThesis/optDesign/csv/Boscia/boscia_" * criterion * "_" * string(m) * "_" * type *"_optimality.csv" + end if !isfile(file_name) CSV.write(file_name, df, append=true, writeheader=true, delim=";") else @@ -203,10 +236,10 @@ function solve_opt(seed, m, n, time_limit, criterion, corr; full_callback=true, end @show result[:primal_objective] @show result[:dual_gap] + @show result[:primal_objective] - result[:dual_gap] @show x @show sum(x) @show N - # @assert isapprox(sum(x), N) end return x, result diff --git a/ODWB/src/opt_design_custom_BB.jl b/ODWB/src/opt_design_custom_BB.jl index 6040425..3bf5152 100644 --- a/ODWB/src/opt_design_custom_BB.jl +++ b/ODWB/src/opt_design_custom_BB.jl @@ -1,8 +1,5 @@ """ Custom Branch-and-Bound for the Exact Optimal Design Problem - -From "A Branch-and-Bound Algorithm for the Exact Optimal Experimental Design Problem" -by Selin Damla Ahipasaoglu """ """ @@ -13,6 +10,7 @@ mutable struct CustomBBNode <: Bonobo.AbstractNode lower_bounds::Vector{Float64} upper_bounds::Vector{Float64} solution::Vector{Float64} + w::Vector{Float64} time::Float64 iteration_count::Int64 end @@ -45,32 +43,40 @@ Create the information for the childern nodes. """ function Bonobo.get_branching_nodes_info(tree::Bonobo.BnBTree{<:CustomBBNode}, node, vidx) x = Bonobo.get_relaxed_values(tree, node) - @assert node.lower_bounds[vidx] < node.upper_bounds[vidx] + w = node.w + N = tree.root.N + @assert node.lower_bounds[vidx] < node.upper_bounds[vidx] "Branching lb: $(node.lower_bounds[vidx]) ub: $(node.upper_bounds[vidx]) x[vidx]=$(x[vidx])" frac_val = x[vidx] # Left child keeps the lower bounds and gets a new upper bound at vidx. left_bounds = copy(node.upper_bounds) - left_bounds[vidx] = floor(frac_val) + left_bounds[vidx] = floor(frac_val)/N left_solution = copy(node.solution) left_solution[vidx] = floor(frac_val) + left_w = copy(node.w) + left_w[vidx] = floor(frac_val)/N node_info_left = (lower_bounds=node.lower_bounds, upper_bounds = left_bounds, solution = left_solution, + w = left_w, time = 0.0, iteration_count = 0 ) # Right child keeps the upper bounds and gets a new lower bound at vidx. right_bounds = copy(node.lower_bounds) - right_bounds[vidx] = ceil(frac_val) + right_bounds[vidx] = ceil(frac_val)/N right_solution = copy(node.solution) right_solution[vidx] = ceil(frac_val) + right_w = copy(node.w) + right_w[vidx] = ceil(frac_val)/N node_info_right = (lower_bounds=right_bounds, upper_bounds = node.upper_bounds, solution = right_solution, + w = right_w, time = 0.0, iteration_count = 0 ) @@ -87,35 +93,35 @@ function Bonobo.evaluate_node!(tree, node) iter_count = 0 # get feasible start point - x = copy(node.solution) - x = project_onto_feasible_region(x, node, tree.root.N) + N = tree.root.N + w = copy(node.w) + w = alternating_projection(w, node, 1) # If the node is not reachable, i.e. we can never satisfy the knapsack add_constraint - if x === nothing + if w === nothing return NaN, NaN end + x = N*w - @assert check_feasibilty(x, node, tree) - - gradient = similar(x) - tree.root.grad(gradient, x) - linesearch_workspace = FrankWolfe.build_linesearch_workspace(FrankWolfe.Adaptive(), x, gradient) + gradient = similar(w) + tree.root.grad(gradient, w) + linesearch_workspace = FrankWolfe.build_linesearch_workspace(FrankWolfe.Adaptive(), w, gradient) - jdx = find_max_gradient(x, -gradient, node.upper_bounds) - kdx = find_min_gradient(x, -gradient, node.lower_bounds) + jdx = find_max_gradient(w, -gradient, node.upper_bounds, tree.root) + kdx = find_min_gradient(w, -gradient, node.lower_bounds, tree.root) - while gradient[jdx]/gradient[kdx] - 1 > 1e-5 && iter_count < 10000 + while gradient[jdx]/gradient[kdx] - 1 > 1e-6 && iter_count < 10000 # find theta - d = zeros(length(x)) + d = zeros(length(w)) d[jdx] = 1.0 d[kdx] = -1.0 - theta_max = min(node.upper_bounds[jdx] - x[jdx], x[kdx] - node.lower_bounds[kdx]) + theta_max = min(node.upper_bounds[jdx] - w[jdx], w[kdx] - node.lower_bounds[kdx]) theta = tree.root.linesearch( node, tree.root.f, tree.root.grad, gradient, - x, + w, d, theta_max, linesearch_workspace, @@ -125,33 +131,35 @@ function Bonobo.evaluate_node!(tree, node) ) # weight shift - x[jdx] += theta - x[kdx] -= theta + w[jdx] += theta + w[kdx] -= theta # compute the new jdx and kdx - tree.root.grad(gradient, x) - jdx = find_max_gradient(x, -gradient, node.upper_bounds) - kdx = find_min_gradient(x, -gradient, node.lower_bounds) + tree.root.grad(gradient, w) + jdx = find_max_gradient(w, -gradient, node.upper_bounds, tree.root) + kdx = find_min_gradient(w, -gradient, node.lower_bounds, tree.root) iter_count += 1 - if !check_feasibilty(x, node, tree) - @show node.upper_bounds - x - @show x - node.lower_bounds - @show sum(x) + if !check_feasibilty(w, node) + @show node.upper_bounds - w + @show w - node.lower_bounds + @show sum(w) end - @assert check_feasibilty(x, node, tree) + @assert check_feasibilty(w, node) end time = float(Dates.value(Dates.now() - time_ref)) node.time = time node.iteration_count = iter_count - obj_value = tree.root.f(x) + x = N*w + node.w = w node.solution = x + obj_value = tree.root.f(x) - if isapprox(sum(isapprox.(x, round.(x); atol=1e-6, rtol=5e-2)), tree.root.nvars) && check_feasibilty(round.(x), node, tree) - #println("Integer solution found.") + if isapprox(sum(isapprox.(x, round.(x); atol=1e-6, rtol=5e-2)), tree.root.nvars) && check_feasibilty(round.(x),node.lower_bounds*N, node.upper_bounds*N, N) + println("Integer solution found.") node.solution = round.(x) primal = tree.root.f(node.solution) return obj_value, primal @@ -161,34 +169,31 @@ function Bonobo.evaluate_node!(tree, node) end """ -Project point on the feasible region which is the intersection of the hyperbox -and the probability simplex scaled by integer N. + Alternating projection + +See: https://en.wikipedia.org/wiki/Projections_onto_convex_sets """ -function project_onto_feasible_region(x, node, N) - # the difference in weight you need to distribute - s_prod = sum(x) - N - # we are already on the probability simplex - if s_prod == 0 - return x - end - # If the difference is negative, we might crash the upper bounds. - # Else, we could crash the lower bounds. - bounds = s_prod > 0 ? node.lower_bounds : node.upper_bounds - b = (x - s_prod*ones(length(x))-bounds)*sign(s_prod) - set = findall(x -> x > 0, b) - - # If the set is empty, we reached an infeasible node! - if isempty(set) +function alternating_projection(x, node, N) + if sum(node.lower_bounds) > N || sum(node.upper_bounds) < N || !all(node.lower_bounds <= node.upper_bounds) return nothing end - dir = zeros(length(x)) - dir[set] .= 1 - x = x - s_prod/sum(dir) * dir + function project_onto_cube(x, lb, ub) + return clamp.(x, lb, ub) + end + function project_on_prob_simplex(x, N) + n = length(x) + return x - (sum(x)-N)/n * ones(n) + end + + iter = 0 - @assert isapprox(sum(x), N; atol=1e-4, rtol=1e-2) - @assert sum(node.upper_bounds - x .>= 0) == length(x) - @assert sum(x- node.lower_bounds .>= 0) == length(x) + while !check_feasibilty(x, node.lower_bounds, node.upper_bounds, N) && iter < 100000 + x = project_onto_cube(x, node.lower_bounds, node.upper_bounds) + x = project_on_prob_simplex(x, N) + iter += 1 + end + @debug "Number of iteration in the alternating projection: $(iter)" return x end @@ -196,17 +201,36 @@ end """ Check feasibility of a given point. """ -function check_feasibilty(x, lb, ub, N) +function check_feasibilty(x, lb, ub, N; tol=0.0) m = length(x) + if all(x.>=lb .- tol) && all(x.<=ub .+ tol) && isapprox(sum(x), N) + return true + else + @debug "Sum of w $(sum(x))" + @debug "Lower bounds: $(lb)" + @debug "Iterate: $(x)" + @debug "Upper bounds: $(ub)" + @debug "Difference iterate and lower bounds: $(x-lb)" + @debug "Difference iterate and upper bounds: $(ub-x)" + return false + end return sum(x.>=lb) == m && sum(x.<=ub) == m && isapprox(sum(x), N) end -check_feasibilty(x, node, tree) = check_feasibilty(x, node.lower_bounds, node.upper_bounds, tree.root.N) +check_feasibilty(x, node; N=1, tol=1e-9) = check_feasibilty(x, node.lower_bounds, node.upper_bounds, N, tol=tol) + +""" +Returns the list of indices corresponding to the integral variables. +""" +function Bonobo.get_branching_indices(root::ConstrainedBoxMIProblem) + return root.integer_vars +end """ Find maximum gradient entry where the upper bound has yet been met. """ -function find_max_gradient(x, gradient, ub) - bool_vec = x .< ub +function find_max_gradient(x, gradient, ub, root; tol=0.0) + int_var = Bonobo.get_branching_indices(root) + bool_vec = x[int_var] .< ub[int_var] .+ tol idx_set = findall(x -> x == 1, bool_vec) return findfirst(x -> x == maximum(gradient[idx_set]), gradient) end @@ -214,50 +238,52 @@ end """ Find minimum gradient entry where the lower bound has yet been met. """ -function find_min_gradient(x, gradient, lb) - bool_vec = x .> lb +function find_min_gradient(x, gradient, lb, root; tol=0.0) + int_var = Bonobo.get_branching_indices(root) + bool_vec = x[int_var] .> lb[int_var] .- tol idx_set = findall(x -> x == 1, bool_vec) return findfirst(x -> x == minimum(gradient[idx_set]), gradient) end -""" -Returns the list of indices corresponding to the integral variables. -""" -function Bonobo.get_branching_indices(root::ConstrainedBoxMIProblem) - return root.integer_vars -end - - """ Set up the problem and solve it with the custom BB algorithm. """ -function solve_opt_custom(seed, m, n, time_limit, criterion, corr; write = true, verbose= true) +function solve_opt_custom(seed, m, n, time_limit, criterion, corr; p=0, write = true, verbose= true, long_run=false, print_iter=100, specific_seed=false) # build data - A, C, N, ub = build_data(seed, m, n, criterion in ["AF","DF"], corr) + A, C, N, ub, C_hat = build_data(seed, m, n, criterion in ["AF","DF", "GTIF"], corr, scaling_C=long_run) # build function - p = 1 if criterion == "A" || criterion == "AF" p = -1 elseif criterion == "D" || criterion == "DF" p = 0 end @assert p <= 0 - if criterion in ["A","D"] - C = nothing + if criterion in ["A","D","GTI"] + C_hat = nothing end - f, grad!, linesearch = build_matrix_means_objective(A, p,C=C) + build_safe = criterion in ["A","D","GTI"] + f, grad!, linesearch = build_matrix_means_objective(A, p, C_hat=C_hat, build_safe = build_safe) # create problem and tree - if criterion in ["A","D"] + if criterion in ["A","D", "GTI"] x0 = greedy_incumbent(A, m, n, N, ub) + m_hat = m + ub_hat = ub + N_hat = N + lb_hat = zeros(m) else - x0= greedy_incumbent_fusion(A,m,n,N,ub) + x0 = vcat(greedy_incumbent_fusion(A,m,n,N,ub), fill(1, 2n)) + m_hat = m + 2n + ub_hat = vcat(ub, ones(2n)) + N_hat = N + 2n + lb_hat = vcat(zeros(m), ones(2n)) end - root = ConstrainedBoxMIProblem(f, grad!, linesearch, m, collect(1:m), N, time_limit) + root = ConstrainedBoxMIProblem(f, grad!, linesearch, m_hat, collect(1:m), N_hat, time_limit) nodeExample = CustomBBNode( Bonobo.BnBNodeInfo(1, 0.0, 0.0), - zeros(m), - ub, + lb_hat, + ub_hat, + fill(-1.0, length(x0)), fill(-1.0, length(x0)), 0.0, 0 @@ -267,6 +293,7 @@ function solve_opt_custom(seed, m, n, time_limit, criterion, corr; write = true, Value = Vector{Float64} tree = Bonobo.initialize(; traverse_strategy = Bonobo.BFS(), + branch_strategy = Bonobo.MOST_INFEASIBLE(), Node = Node, Solution = Bonobo.DefaultSolution{Node, Value}, root = root, @@ -274,15 +301,16 @@ function solve_opt_custom(seed, m, n, time_limit, criterion, corr; write = true, ) Bonobo.set_root!(tree, ( - lower_bounds = zeros(length(x0)), - upper_bounds = ub, + lower_bounds = lb_hat/N_hat, + upper_bounds = ub_hat/N_hat, solution = x0, + w = x0/N_hat, time=0.0, iteration_count=0 )) time_ref = Dates.now() - callback = build_callback(tree, time_ref, verbose) + callback = build_callback(tree, time_ref, verbose, print_iter) # dummy solution - In case the process stops because of the time limit and no solution as been found yet. dummy_solution = Bonobo.DefaultSolution(Inf, fill(-1.0, length(x0)), nodeExample) @@ -306,43 +334,56 @@ function solve_opt_custom(seed, m, n, time_limit, criterion, corr; write = true, if tree.root.Stage == Boscia.OPT_TREE_EMPTY @assert sum(sol_x .!= dummy_solution.solution) == length(sol_x) end - @show sol_x + @show sol_x[1:m] solution = f(sol_x) @show solution + @show tree.num_nodes # check feasibility - feasible = isfeasible(seed, m, n, criterion, sol_x, corr, N=N) + feasible = isfeasible(seed, m, n, criterion, sol_x[1:m], corr, N=N) if !feasible - @show m, sum(ub - sol_x .>= 0) - @show m, sum(sol_x .>= 0) - @show N, sum(sol_x) + @show m, sum(ub - sol_x[1:m] .>= 0) + @show m, sum(sol_x[1:m] .>= 0) + @show N, sum(sol_x[1:m]) end @assert feasible # Calculate objective with respect to the criteria used in Boscia if criterion in ["A","AF"] - f_check, _ = build_a_criterion(A, criterion == "AF", C=C) + f_check, _ = build_a_criterion(A, criterion == "AF", C=C, build_safe = criterion=="A") else - f_check, _ = build_d_criterion(A, criterion == "DF", C=C) + f_check, _ = build_d_criterion(A, criterion == "DF", C=C, build_safe = criterion=="D") end - solution_scaled = f_check(sol_x) - #solution_scaled = criterion == "A" ? exp(solution) : solution - n*log(n) + solution_scaled = f_check(sol_x[1:m]) @show solution_scaled # Check the solving stage @assert tree.root.Stage != Boscia.SOLVING + @show tree.root.Stage status = if tree.root.Stage in [Boscia.OPT_GAP_REACHED, Boscia.OPT_TREE_EMPTY] "optimal" else - "Time limit reached" + "TIME_LIMIT" end + dual_gap = tree.incumbent - tree.lb + @show dual_gap + if write + if criterion in ["GTI","GTIF"] + criterion = criterion * "_" * string(Int64(p*100)) + end # write data into file time_per_nodes = time/tree.num_nodes type = corr ? "correlated" : "independent" - df = DataFrame(seed=seed, numberOfExperiments=m, numberOfParameters=n, N=N, time=time, time_per_nodes=time_per_nodes, solution=solution, solution_scaled=solution_scaled, number_nodes=tree.num_nodes, termination=status) - file_name = joinpath(@__DIR__, "../csv/CustomBB/customBB_" * criterion * "_" * string(m) * "_" * type * "_optimality.csv") + df = DataFrame(seed=seed, numberOfExperiments=m, numberOfParameters=n, N=N, time=time, time_per_nodes=time_per_nodes, solution=solution, solution_scaled=solution_scaled, dual_gap=dual_gap, number_nodes=tree.num_nodes, termination=status) + if long_run + file_name = "/home/htc/dhendryc/research_projects/MasterThesis/optDesign/csv/CustomBB/long_runs/custombb_" * criterion * "_" * string(m) * "_" * type * "_optimality.csv" + elseif specific_seed + file_name = "/home/htc/dhendryc/research_projects/MasterThesis/optDesign/csv/CustomBB/custombb_" * criterion * "_" * string(m) * "_" * string(n) * "_" * type * "_" * string(seed) * "_optimality.csv" + else + file_name = "/home/htc/dhendryc/research_projects/MasterThesis/optDesign/csv/CustomBB/custombb_" * criterion * "_" * string(m) * "_" * type * "_optimality.csv" + end if !isfile(file_name) CSV.write(file_name, df, append=true, writeheader=true) else @@ -352,24 +393,28 @@ function solve_opt_custom(seed, m, n, time_limit, criterion, corr; write = true, return sol_x end -function build_callback(tree, time_ref, verbose) +function build_callback(tree, time_ref, verbose, print_iter) + iter = 1 return function callback(tree, node; node_infeasible = false, worse_than_incumbent = false) time = float(Dates.value(Dates.now() - time_ref)) - if (node.id == 1 || mod(tree.num_nodes, 100) == 0 ) && verbose - println("ID: $(node.id) num_nodes: $(tree.num_nodes) LB: $(tree.lb) Incumbent: $(tree.incumbent) Total Time in s: $(time/1000.0) Time node in s: $(node.time/1000.0) Iterations: $(node.iteration_count)") + if (node.id == 1 || mod(iter, print_iter) == 0 || Bonobo.terminated(tree)) && verbose + println("Iter: $(iter) ID: $(node.id) num_nodes: $(tree.num_nodes) LB: $(tree.lb) Incumbent: $(tree.incumbent) Total Time in s: $(time/1000.0) Time node in s: $(node.time/1000.0) Iterations: $(node.iteration_count)") if node_infeasible println("Node not feasible!") elseif worse_than_incumbent println("Node cut because it is worse than the incumbent") end end + iter += 1 # break if time is met if tree.root.time_limit < Inf if time / 1000.0 ≥ tree.root.time_limit - @assert tree.root.Stage == Boscia.SOLVING - tree.root.Stage = Boscia.TIME_LIMIT_REACHED + if tree.root.Stage == Boscia.SOLVING + @assert tree.root.Stage == Boscia.SOLVING + tree.root.Stage = Boscia.TIME_LIMIT_REACHED + end end end end @@ -389,6 +434,7 @@ function Bonobo.terminated(tree::Bonobo.BnBTree{<:CustomBBNode}) absgap = tree.incumbent - tree.lb if absgap ≤ tree.root.abs_gap tree.root.Stage = Boscia.OPT_GAP_REACHED + @show absgap return true end @@ -402,6 +448,7 @@ function Bonobo.terminated(tree::Bonobo.BnBTree{<:CustomBBNode}) end if dual_gap ≤ tree.root.rel_gap tree.root.Stage = Boscia.OPT_GAP_REACHED + @show dual_gap return true end diff --git a/ODWB/src/opt_design_scip.jl b/ODWB/src/opt_design_scip.jl index 7a00ab4..18269bd 100644 --- a/ODWB/src/opt_design_scip.jl +++ b/ODWB/src/opt_design_scip.jl @@ -1,14 +1,13 @@ -### Optimal design with SCIP -function solve_opt_scip(seed, m, n, time_limit, criterion, corr; write=true, verbose=true) - if criterion == "A" - error("SCIP OA does not work with A-opt") - elseif criterion == "D" - error("SCIP OA does not work with D-opt") +# Optimal design with SCIP +function solve_opt_scip(seed, m, n, time_limit, criterion, corr; p=0, write=true, verbose=true, long_run=false, specific_seed = false) + if criterion in ["A","D"] + error("SCIP OA only works with the Fusion Problems") elseif criterion == "AF" - A, C, N, ub =build_data(seed, m, n, true, corr) + A, C, N, ub, _ = build_data(seed, m, n, true, corr) f, grad! = build_a_criterion(A, true, C=C) + #f, grad! = build_general_trace(A, -1, true, C=C) # Log(Tr(X^{-1})) elseif criterion == "DF" - A, C, N, ub =build_data(seed, m, n, true, corr) + A, C, N, ub, _ = build_data(seed, m, n, true, corr) f, grad! = build_d_criterion(A, true, C=C) else error("Invalid criterion!") @@ -34,7 +33,13 @@ function solve_opt_scip(seed, m, n, time_limit, criterion, corr; write=true, ver if write df = DataFrame(seed=seed, numberOfExperiments=m, numberOfParameters=n, N=N,time=time_scip, solution=solution_scip, dual=dual_objective, termination=termination_scip, calls=ncalls_scip) - file_name = joinpath(@__DIR__, "../csv/SCIP/scip_" * criterion * "_" * string(m) * "_" * type * "_optimality" * ".csv") + if long_run + file_name = "/home/htc/dhendryc/research_projects/MasterThesis/optDesign/csv/SCIP/long_run/scip_" * criterion * "_" * string(m) * "_" * type * "_optimality" * ".csv" + elseif specific_seed + file_name = file_name = "/home/htc/dhendryc/research_projects/MasterThesis/optDesign/csv/SCIP/scip_" * criterion * "_" * string(m) * "_" * string(n) * "_" * type * "_optimality_" * string(seed) * ".csv" + else + file_name = "/home/htc/dhendryc/research_projects/MasterThesis/optDesign/csv/SCIP/scip_" * criterion * "_" * string(m) * "_" * type * "_optimality_" * ".csv" + end if !isfile(file_name) CSV.write(file_name, df, append=true, writeheader=true) else @@ -62,8 +67,3 @@ function build_scip_optimizer(m, N, ub, limit, f, grad!, verbose) return lmo, epigraph_ch, x, lmo_check end - - - - - diff --git a/ODWB/src/utilities.jl b/ODWB/src/utilities.jl index c8cffd7..3e83c00 100644 --- a/ODWB/src/utilities.jl +++ b/ODWB/src/utilities.jl @@ -7,7 +7,8 @@ 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(seed, m, n, fusion, corr) +function build_data(seed, m, n, fusion, corr; scaling_C=false) + @show scaling_C # set up Random.seed!(seed) if corr @@ -22,8 +23,9 @@ function build_data(seed, m, n, fusion, corr) A = rand(m,n) @assert rank(A) == n # check that A has the desired rank! end - C = rand(2n, n) - C = transpose(C)*C + C_hat = rand(2n, n) + C = scaling_C ? 1/2n*transpose(C_hat)*C_hat : transpose(C_hat)*C_hat + @assert rank(C) == n if fusion @@ -35,14 +37,14 @@ function build_data(seed, m, n, fusion, corr) ub = rand(1.0:u, m) end - return A, C, N, ub + return A, C, N, ub, C_hat end """ Build LMO for the problems. Used in Boscia and SCIP. """ function build_lmo(o, m, N, ub) - MOI.set(o, MOI.Silent(), true) + MOI.set(o, MOI.Silent(), false) MOI.empty!(o) x = MOI.add_variables(o, m) for i in 1:m @@ -72,9 +74,12 @@ end """ Build function for the A-criterion. """ -function build_a_criterion(A, fusion; μ=0.0, C=nothing, build_safe=false) +function build_a_criterion(A, fusion; μ=1e-4, C=nothing, build_safe=false, long_run=false) m, n = size(A) - a=m + a = m + if !fusion && m in [100,120] + μ = 1e-3 + end domain_oracle = build_domain_oracle(A, n) if fusion && C === nothing @@ -90,15 +95,14 @@ function build_a_criterion(A, fusion; μ=0.0, C=nothing, build_safe=false) end function grad_a!(storage, x) - #x = BigFloat.(x) # Setting can be useful for numerical tricky problems X = fusion ? C + transpose(A)*diagm(x)*A : 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 #float.(storage) # in case of x .= BigFloat(x) - end + return storage + end function f_a_safe(x) if !domain_oracle(x) @@ -114,14 +118,13 @@ function build_a_criterion(A, fusion; μ=0.0, C=nothing, build_safe=false) if !domain_oracle(x) return fill(Inf, length(x)) end - #x = BigFloat.(x) # Setting can be useful for numerical tricky problems X = fusion ? C + transpose(A)*diagm(x)*A : 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 #float.(storage) # in case of x .= BigFloat(x) + return storage end if build_safe @@ -134,9 +137,10 @@ end """ Build function for the D-criterion. """ -function build_d_criterion(A, fusion; μ =1e-5, C=nothing, build_safe=false) +function build_d_criterion(A, fusion; μ =0.0, C=nothing, build_safe=false, long_run=false) m, n = size(A) - a=m + a = m + γ = long_run ? m : 1 domain_oracle = build_domain_oracle(A, n) if fusion && C === nothing @@ -145,18 +149,17 @@ function build_d_criterion(A, fusion; μ =1e-5, C=nothing, build_safe=false) function f_d(x) X = fusion ? C + transpose(A)*diagm(x)*A : transpose(A)*diagm(x)*A + Matrix(μ *I, n, n) - X = Symmetric(X) - return -log(det(X))/a + X = Symmetric(X * 1/γ) + return float(-log(det(X))/a) end function grad_d!(storage, x) X = fusion ? C + transpose(A)*diagm(x)*A : transpose(A)*diagm(x)*A + Matrix(μ *I, n, n) - X= Symmetric(X) + X= Symmetric(X * 1/γ) F = cholesky(X) for i in 1:length(x) - storage[i] = 1/a * LinearAlgebra.tr(-(F \ A[i,:] )*transpose(A[i,:])) + storage[i] = 1/a * LinearAlgebra.tr(-(F \ A[i,:] )*transpose(A[i,:])) * 1/γ end - # https://stackoverflow.com/questions/46417005/exclude-elements-of-array-based-on-index-julia return storage end @@ -165,8 +168,8 @@ function build_d_criterion(A, fusion; μ =1e-5, C=nothing, build_safe=false) return Inf end X = fusion ? C + transpose(A)*diagm(x)*A : transpose(A)*diagm(x)*A + Matrix(μ *I, n, n) - X = Symmetric(X) - return -log(det(X))/a + X = Symmetric(X*1/γ) + return float(-log(det(X))/a) end function grad_d_safe!(storage, x) @@ -174,10 +177,10 @@ function build_d_criterion(A, fusion; μ =1e-5, C=nothing, build_safe=false) return fill(Inf, length(x)) end X = fusion ? C + transpose(A)*diagm(x)*A : transpose(A)*diagm(x)*A + Matrix(μ *I, n, n) - X= Symmetric(X) + X= Symmetric(X*1/γ) F = cholesky(X) for i in 1:length(x) - storage[i] = 1/a * LinearAlgebra.tr(-(F \ A[i,:] )*transpose(A[i,:])) + storage[i] = 1/a * LinearAlgebra.tr(-(F \ A[i,:] )*transpose(A[i,:])) * 1/γ end # https://stackoverflow.com/questions/46417005/exclude-elements-of-array-based-on-index-julia return storage @@ -192,69 +195,75 @@ end function build_general_trace(A, p, fusion; C=nothing, μ=0.0, build_safe=false) m, n = size(A) - a=m + a = 1 domain_oracle = build_domain_oracle(A, n) if fusion && C === nothing @error("For the fusion problem, please provide a matrix C.") end - function f_a(x) + function f_gti(x) X = fusion ? C + transpose(A)*diagm(x)*A : transpose(A)*diagm(x)*A + Matrix(μ *I, n, n) - X = Symmetric(X^p) - return LinearAlgebra.tr(X)/a + X= Symmetric(X) + if p == 0 + return -log(det(X)) + end + return -1/p * log(LinearAlgebra.tr(Symmetric(X^p))) # 1/n * end - function grad_a!(storage, x) - #x = BigFloat.(x) # Setting can be useful for numerical tricky problems + function grad_gti!(storage, x) X = fusion ? C + transpose(A)*diagm(x)*A : transpose(A)*diagm(x)*A + Matrix(μ *I, n, n) - X = Symmetric(X^(p-1)) - #F = cholesky(X) - for i in 1:length(x) - storage[i] = (p/a) * A[i,:]' * X * A[i,:] + X=Symmetric(X) + a = p == 0 ? -1 : -1/(LinearAlgebra.tr(Symmetric(X^p))) + X =Symmetric(X^(p-1)) + for i in 1:m + storage[i] = a* A[i,:]' * X * A[i,:] end - return storage #float.(storage) # in case of x .= BigFloat(x) + return storage end - function f_a_safe(x) + function f_gti_safe(x) if !domain_oracle(x) return Inf end X = fusion ? C + transpose(A)*diagm(x)*A : transpose(A)*diagm(x)*A + Matrix(μ *I, n, n) - X = Symmetric(X^p) - # X_inv = LinearAlgebra.inv(X) - return LinearAlgebra.tr(X)/a + X= Symmetric(X) + if p == 0 + return -log(det(1/n*X)) + end + return -1/p * log(LinearAlgebra.tr(Symmetric(X^p))) # 1/n * end - function grad_a_safe!(storage, x) + function grad_gti_safe!(storage, x) if !domain_oracle(x) return fill(Inf, length(x)) end - #x = BigFloat.(x) # Setting can be useful for numerical tricky problems X = fusion ? C + transpose(A)*diagm(x)*A : transpose(A)*diagm(x)*A + Matrix(μ *I, n, n) - X = Symmetric(X^(p-1)) - # F = cholesky(X) - for i in 1:length(x) - storage[i] = (p/a) * A[i,:]' * X * A[i,:] + X=Symmetric(X) + a = p == 0 ? -1 : -1/(LinearAlgebra.tr(Symmetric(X^p))) + X =Symmetric(X^(p-1)) + for i in 1:m + storage[i] = a* A[i,:]' * X * A[i,:] end - return storage #float.(storage) # in case of x .= BigFloat(x) + return storage end if build_safe - return f_a_safe, grad_a_safe! + return f_gti_safe, grad_gti_safe! end - return f_a, grad_a! + return f_gti, grad_gti! end """ Build general matrix means objective: log(ϕ(X)) """ -function build_matrix_means_objective(A,p;C=nothing) +function build_matrix_means_objective(A,p;μ=0.0, C_hat=nothing, build_safe=false) m,n = size(A) + domain_oracle = build_domain_oracle(A, n) function inf_matrix(x) - X = C===nothing ? A' * diagm(x) * A : C + A' * diagm(x) * A + X = C_hat===nothing ? A' * diagm(x) * A + Matrix(μ *I, n, n) : C_hat' * diagm(x[m+1:m+2n]) * C_hat + A' * diagm(x[1:m]) * A return X end @@ -278,6 +287,33 @@ function build_matrix_means_objective(A,p;C=nothing) return storage end + function f_safe(x) + if !domain_oracle(x) + return Inf + end + X = inf_matrix(x) + X=Symmetric(X) + if p == 0 + @assert det(1/n*X) > 0 "Determinat is $(det(1/n*X)) and x is $(x)" + return -log(det(1/n*X)) + end + return -1/p * log(LinearAlgebra.tr(Symmetric(X^p))) # 1/n * + end + + function grad_safe!(storage, x) + if !domain_oracle(x) + return fill(Inf, length(x)) + end + X = inf_matrix(x) + X=Symmetric(X) + a = p == 0 ? -1 : -1/(LinearAlgebra.tr(Symmetric(X^p))) + X =Symmetric(X^(p-1)) + for i in 1:m + storage[i] = a* A[i,:]' * X * A[i,:] + end + return storage + end + function linesearch(node, f, grad!, gradient, x, d, theta_max, linesearch_workspace, iter_count, jdx, kdx) theta = 0.0 X = inf_matrix(x) @@ -287,13 +323,11 @@ function build_matrix_means_objective(A,p;C=nothing) w_k = A[kdx,:]' * X_inv * A[kdx,:] # D Optimality if p == 0 - if isapprox(w_jk^2, w_k*w_j) - theta = min(node.upper_bounds[jdx] - x[jdx], x[kdx] - node.lower_bounds[kdx]) - elseif w_jk^2 < w_k*w_j + if w_jk^2 < w_k*w_j theta_bar = (w_j-w_k) / (2(w_j*w_k - w_jk^2)) theta = min(node.upper_bounds[jdx] - x[jdx], x[kdx] - node.lower_bounds[kdx], theta_bar) else - error("W_jk^2 : $(w_jk^2) is greater than w_j*w_k: $(w_j*w_k)") + theta = min(node.upper_bounds[jdx] - x[jdx], x[kdx] - node.lower_bounds[kdx]) end # A Optimality elseif p == -1 @@ -325,7 +359,7 @@ function build_matrix_means_objective(A,p;C=nothing) grad!, gradient, x, - d, + -d, theta_max, linesearch_workspace, FrankWolfe.InplaceEmphasis(), @@ -334,6 +368,9 @@ function build_matrix_means_objective(A,p;C=nothing) return theta end + if build_safe + return f_safe, grad_safe!, linesearch + end return f, grad!, linesearch end @@ -419,10 +456,10 @@ function linearly_independent_rows(A, m ,n) end """ -Build start point used in Boscia in case of A-opt and D-opt. +Build start point used in FrankWolfe and Boscia in case of A-opt and D-opt. 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_point2(A, m, n, N, ub) # Get n linearly independent rows of A S = linearly_independent_rows(A,m,n) @@ -432,7 +469,6 @@ function build_start_point2(A, m, n, N, ub) E = [] V = Vector{Float64}[] - while !isempty(setdiff(S, findall(x-> !(iszero(x)),x))) v = zeros(m) while sum(v) < N @@ -511,23 +547,21 @@ end Check if a point is linear feasible with respect to the original model """ function isfeasible(seed, m, n, criterion, x, corr; N=0) - if criterion == "A" || criterion == "D" - A, _, N, ub = build_data(seed, m, n, false, corr) - elseif criterion == "AF"|| criterion == "DF" - A, C, N, ub = build_data(seed, m, n, true, corr) + if criterion in ["A","D","GTI"] + A, _, N, ub, _ = build_data(seed, m, n, false, corr) + elseif criterion in ["AF","DF","GTIF"] + A, C, N, ub, _ = build_data(seed, m, n, true, corr) end if sum(x) < N - 1e-4 return false elseif sum(x.>=0-1e-4) != m return false - elseif sum(ub-x.>= 0-1e-4) != m + elseif sum(ub - x.>= 0-1e-4) != m return false - elseif criterion in ["D", "A"] && m - sum(iszero.(x)) < n + elseif criterion in ["D", "A", "GTI"] && m - sum(iszero.(x)) < n return false else return true end end - - diff --git a/ODWB/test/test_derivatives.jl b/ODWB/test/test_derivatives.jl index e09042b..c52d3d4 100644 --- a/ODWB/test/test_derivatives.jl +++ b/ODWB/test/test_derivatives.jl @@ -6,6 +6,10 @@ using Random using ODWB using Test +seed = rand(UInt64) +@show seed +Random.seed!(seed) + """ Check if the gradient using finite differences matches the grad! provided. Copied from FrankWolfe package: https://github.com/ZIB-IOL/FrankWolfe.jl/blob/master/examples/plot_utils.jl @@ -29,7 +33,7 @@ end n = Int(floor(dim/4)) @show dim, n gradient = rand(dim) - A, _, _, _ = ODWB.build_data(seed,dim, n, false, false) + A, _, _, _, _ = ODWB.build_data(seed,dim, n, false, false) f, grad! = ODWB.build_a_criterion(A, false, μ=1e-2) @test check_gradients(grad!, f, gradient) @@ -41,8 +45,8 @@ end n = Int(floor(dim/4)) @show dim, n gradient = rand(dim) - A, _, _, _ = ODWB.build_data(seed,dim, n, false, false) - f, grad! = ODWB.build_d_criterion(A, false) + A, _, _, _, _ = ODWB.build_data(seed,dim, n, false, false) + f, grad! = ODWB.build_d_criterion(A, false, build_safe=true) @test check_gradients(grad!, f, gradient) end @@ -53,7 +57,7 @@ end n = Int(floor(dim/4)) @show dim, n gradient = rand(dim) - A, C, _, _ = ODWB.build_data(seed,dim, n, true, false) + A, C, _, _, _ = ODWB.build_data(seed,dim, n, true, false) f, grad! = ODWB.build_a_criterion(A, true, C=C) @test check_gradients(grad!, f, gradient) @@ -65,7 +69,7 @@ end n = Int(floor(dim/4)) @show dim, n gradient = rand(dim) - A, C, _, _ = ODWB.build_data(seed,dim, n, true, false) + A, C, _, _, _ = ODWB.build_data(seed,dim, n, true, false) f, grad! = ODWB.build_d_criterion(A,true, C=C) @test check_gradients(grad!, f, gradient) @@ -79,7 +83,12 @@ end n = Int(floor(dim/4)) @show dim, n gradient = rand(dim) - A, _, _, _ = ODWB.build_data(seed,dim, n, true, false) + A, _, _, _, _ = ODWB.build_data(seed,dim, n, true, false) + f, grad! = ODWB.build_matrix_means_objective(A,p) + + @test check_gradients(grad!, f, gradient) + + A, _, _, _, _ = ODWB.build_data(seed,dim, n, false, false) f, grad! = ODWB.build_matrix_means_objective(A,p) @test check_gradients(grad!, f, gradient)