Improve BnB

interfaces/daqp-julia/test/core_tests.jl

@@ -156,6 +156,48 @@ end
     x,_,_,info = DAQPBase.quadprog(H,f,A,bu,bl,sense)
     @test norm(x-[0;1;1]) < tol
 
+    # A binary constraint can be integer feasible at a zero-dual endpoint
+    # without belonging to the active set. Do not branch in that case.
+    ndeg = 8
+    Hdeg = Matrix{Float64}(I, ndeg, ndeg)
+    fdeg = zeros(ndeg)
+    bdeg_u = ones(ndeg)
+    bdeg_l = zeros(ndeg)
+    sdeg = fill(Cint(DAQPBase.BINARY), ndeg)
+    xdeg, _, efdeg, ideg = DAQPBase.quadprog(
+        Hdeg, fdeg, zeros(0, ndeg), bdeg_u, bdeg_l, sdeg)
+    @test efdeg == DAQPBase.OPTIMAL
+    @test norm(xdeg, Inf) < tol
+    @test ideg.nodes == 1
+
+    # Cover the same zero-dual case for general binary constraints.
+    xgen, _, efgen, igen = DAQPBase.quadprog(
+        Hdeg, fdeg, Hdeg, bdeg_u, bdeg_l, sdeg)
+    @test efgen == DAQPBase.OPTIMAL
+    @test norm(xgen, Inf) < tol
+    @test igen.nodes == 1
+
+    # BnB cleanup may shorten, but must never lengthen, the valid prefix of
+    # the cached triangular solve. Lengthening it can reuse stale xldl/zldl
+    # entries and produce a CSP that violates its own active equalities.
+    dreuse = DAQPBase.Model()
+    Hreuse = Matrix{Float64}(I, 2, 2)
+    Areuse = ones(1, 2)
+    bureuse = [1.0, 1.0, 1.0]
+    blreuse = [0.0, 0.0, 1.0]
+    sreuse = Cint[DAQPBase.BINARY, DAQPBase.BINARY, DAQPBase.EQUALITY]
+    DAQPBase.setup(
+        dreuse, Hreuse, zeros(2), Areuse, bureuse, blreuse, sreuse)
+    wsreuse = unsafe_load(Ptr{DAQPBase.Workspace}(dreuse.work))
+    @test wsreuse.n_active == 1
+    reuse_field = findfirst(==(:reuse_ind), fieldnames(DAQPBase.Workspace))
+    reuse_offset = fieldoffset(DAQPBase.Workspace, reuse_field)
+    reuse_ptr = Ptr{Cint}(Ptr{UInt8}(dreuse.work) + reuse_offset)
+    unsafe_store!(reuse_ptr, 0)
+    ccall((:daqp_node_cleanup_workspace, DAQPBase.libdaqp), Cvoid,
+          (Cint, Ptr{DAQPBase.Workspace}), wsreuse.n_active, dreuse.work)
+    @test unsafe_load(Ptr{DAQPBase.Workspace}(dreuse.work)).reuse_ind == 0
+
 end
 
 @testset "Model interface" begin
@@ -467,6 +509,29 @@ end
     x,fval,exitflag,info = quadprog(H,f,A,bupper,blower,sense; settings=s)
     @test exitflag == DAQPBase.OPTIMAL
     @test norm(xref-x) < tol
+
+    # Node relaxations in BnB can each finish before the inner solver's
+    # periodic timer check. The tree-level check must still enforce the limit.
+    rng_bnb = Xoshiro(1)
+    nt, nbt, mt = 30, 14, 6
+    Qt = randn(rng_bnb, nt, nt)
+    Ht = Matrix(Qt' * Qt / nt + 0.2I)
+    target = 0.15 .+ 0.7rand(rng_bnb, nt)
+    ft = -Ht * target
+    At = zeros(mt, nt)
+    for row in 1:mt
+        At[row, 1:nbt] .= 0.2 .+ rand(rng_bnb, nbt)
+        At[row, nbt+1:end] .= 0.1randn(rng_bnb, nt - nbt)
+    end
+    center = At * target
+    width = 0.15 .+ 0.15rand(rng_bnb, mt)
+    but = vcat(ones(nbt), fill(2.0, nt - nbt), center + width)
+    blt = vcat(zeros(nbt), fill(-2.0, nt - nbt), center - width)
+    st = vcat(fill(Cint(DAQPBase.BINARY), nbt), zeros(Cint, nt - nbt + mt))
+    s = settings(DAQPBase.Model(), Dict(:time_limit => 1e-9))
+    _, _, exitflag, info = quadprog(Ht, ft, At, but, blt, st; settings=s)
+    @test exitflag == DAQPBase.TIMELIMIT
+    @test info.nodes <= 32
 end
 
 @testset "Semi-proximal method" begin

src/bnb.c

@@ -1,4 +1,24 @@
 #include "bnb.h"
+#ifdef PROFILING
+#include "utils.h"
+#endif
+
+static c_float daqp_binary_diff(const int id, DAQPWorkspace* work){
+    int j, disp;
+    c_float diff = 0.5*(work->dupper[id]+work->dlower[id]);
+
+    if(id < work->ms){//Simple bound
+        if(work->Rinv==NULL) diff -= work->u[id]; //Hessian is identity
+        else{
+            for(j=id,disp=id+DAQP_R_OFFSET(id,work->n);j<work->n;j++)
+                diff -= work->Rinv[disp++]*work->u[j];
+        }
+    }
+    else{//General bound; daqp_add_infeasible already computed M*u
+        diff -= work->Mu[id-work->ms];
+    }
+    return diff;
+}
 
 int daqp_bnb(DAQPWorkspace* work){
     int branch_id, exitflag;
@@ -27,6 +47,17 @@ int daqp_bnb(DAQPWorkspace* work){
 
         node = work->bnb->tree+(--work->bnb->n_nodes);
         exitflag = daqp_process_node(node,work); // Solve relaxation
+#ifdef PROFILING
+        // Individual relaxations are often too short to reach the timer check
+        // in daqp_ldp, so also enforce the limit across the BnB tree.
+        if(work->timer != NULL && (work->bnb->nodecount&31)==0){
+            toc((DAQPtimer*)work->timer);
+            if(get_time((DAQPtimer*)work->timer) > work->settings->time_limit){
+                exitflag = DAQP_EXIT_TIMELIMIT;
+                break;
+            }
+        }
+#endif
         // Cut conditions
         if(exitflag==DAQP_EXIT_INFEASIBLE) continue; // Dominance cut
         if(exitflag<0) break; // Inner solver failed => exit loop
@@ -54,7 +85,7 @@ int daqp_bnb(DAQPWorkspace* work){
         work->settings->fval_bound = fval_bound0;
         // Let work->u point to the best feasible solution
         swp_ptr=work->u; work->u= work->xold; work->xold=swp_ptr;
-        return DAQP_EXIT_OPTIMAL;
+        return exitflag < DAQP_EXIT_INFEASIBLE ? exitflag : DAQP_EXIT_OPTIMAL;
     }
 }
 
@@ -74,9 +105,9 @@ int daqp_process_node(DAQPNode* node, DAQPWorkspace* work){
         else{
             daqp_add_upper_lower(node->bin_id,work);
             work->sense[DAQP_REMOVE_LOWER_FLAG(node->bin_id)] |= DAQP_IMMUTABLE; //Equality
-            if(work->sing_ind != DAQP_EMPTY_IND) // Need to cold start to not miss integer feasible
+            if(work->sing_ind != DAQP_EMPTY_IND){ // Need to cold start to not miss integer feasible
                 daqp_setup_cold_bnb(node,work);
-
+            }
         }
         // Add binary constraint
     }
@@ -85,6 +116,9 @@ int daqp_process_node(DAQPNode* node, DAQPWorkspace* work){
     work->bnb->itercount += work->iterations;
 
     if(exitflag == DAQP_EXIT_CYCLE){// Try to repair (cold start)
+        // A cycle can be caused by stale cached forward-substitution data.
+        // Force the retained fixed prefix to be recomputed during repair.
+        work->reuse_ind=0;
         daqp_setup_cold_bnb(node,work);
         exitflag = daqp_ldp(work);
         work->bnb->itercount += work->iterations;
@@ -94,31 +128,28 @@ int daqp_process_node(DAQPNode* node, DAQPWorkspace* work){
 }
 
 int daqp_get_branch_id(DAQPWorkspace* work){
-    int i, j, disp;
-    int branch_id = DAQP_EMPTY_IND;
-    c_float diff;
+    int i;
+    int id = DAQP_EMPTY_IND;
+    c_float diff, dist, tol;
 
     for(i=0; i < work->bnb->nb; i++){
-        branch_id = work->bnb->bin_ids[i];
+        id = work->bnb->bin_ids[i];
         // Skip fixed binary constraints
-        if(DAQP_IS_ACTIVE(branch_id)) continue;
+        if(DAQP_IS_ACTIVE(id)) continue;
 
         // Compute signed distance from midpoint between bounds
-        diff = 0.5*(work->dupper[branch_id]+work->dlower[branch_id]);
-        if(branch_id < work->ms){//Simple bound
-            if(work->Rinv==NULL) diff -= work->u[branch_id]; //Hessian is identity
-            else{
-                for(j=branch_id,disp=branch_id+DAQP_R_OFFSET(branch_id,work->n);j<work->n;j++)
-                    diff -= work->Rinv[disp++]*work->u[j];
-            }
-        }
-        else{//General bound
-            for(j=0,disp=work->n*(branch_id-work->ms);j<work->n;j++)
-                diff -= work->M[disp++]*work->u[j];
-        }
+        diff = daqp_binary_diff(id,work);
+
+        // A zero-dual binary constraint can lie at an endpoint without being
+        // active. It is already integer feasible and does not need branching.
+
+        dist = 0.5*(work->dupper[id]-work->dlower[id])-(diff < 0 ? -diff : diff);
+        tol = work->settings->primal_tol;
+        if(work->scaling != NULL) tol *= work->scaling[id];
+        if(dist <= tol) continue;
 
-        // Branch on the first infeasible binary variable
-        return diff < 0 ? branch_id : DAQP_ADD_LOWER_FLAG(branch_id);
+        // Explore the endpoint nearest to the relaxation first.
+        return diff < 0 ? id : DAQP_ADD_LOWER_FLAG(id);
     }
 
     return DAQP_EMPTY_IND;
@@ -150,7 +181,9 @@ void daqp_node_cleanup_workspace(int n_clean, DAQPWorkspace* work){
     // Reset workspace
     work->sing_ind=DAQP_EMPTY_IND;
     work->n_active=n_clean;
-    work->reuse_ind=n_clean;
+    // Truncation cannot make a previously invalid cached prefix reusable.
+    if(work->reuse_ind > n_clean)
+        work->reuse_ind=n_clean;
 }