Add functionality for counting outer iterations

CMakeLists.txt

@@ -1,5 +1,5 @@
 cmake_minimum_required(VERSION 3.10)
-project(daqp VERSION 0.9.0)
+project(daqp VERSION 0.9.1)
 
 set(CMAKE_POSITION_INDEPENDENT_CODE ON)
 

codegen/codegen.c

@@ -55,7 +55,7 @@ void render_daqp_workspace(DAQPWorkspace* work, const char *fname, const char *d
     }
 
     // Write Hierarchical
-    if(work->nh > 1 ){
+    if(DAQP_IS_HIERARCHICAL(work)){
         fprintf(fh, "#define DAQP_HIERARCHICAL\n");
         fprintf(fh, "extern int %sbreak_points[%d];\n", prefix, work->nh);
         char varname[256];
@@ -88,7 +88,7 @@ void write_daqp_workspace_h(FILE *f, DAQPWorkspace* work, const char* prefix){
     const int ms = work->ms;
     int ntot = n;
     // Account for soft constraints
-    if(work->nh > 1){
+    if(DAQP_IS_HIERARCHICAL(work)){
         int ns = 0, start=0;
         for(i = 0; i < work->nh; i++){
             ns = (ns  > work->break_points[i]-start) ? ns : work->break_points[i]-start;
@@ -216,7 +216,7 @@ void write_daqp_workspace_src(FILE* f, DAQPWorkspace* work, const char* prefix){
     else
         fprintf(f, "&%sbnb_work,\n", prefix);
     // Hierarhical
-    if(work->nh > 1)
+    if(DAQP_IS_HIERARCHICAL(work))
         fprintf(f, "%d,%sbreak_points,\n", work->nh, prefix);
     else
         fprintf(f, "0, NULL,\n");

docs/docs/settings.md

@@ -49,6 +49,7 @@ Table of contents
 |-5|Nonconvex problem|
 |-6|Initial working set overdetermined|
 |-7|Time limit reached|
+|-8|Unsupported problem|
 
 ## Constraint classification
 The type of a constraint is classified through an integer value (called sense), where the bits in this integer encode different properties: 

include/constants.h

@@ -48,6 +48,7 @@ extern "C" {
 #define DAQP_EXIT_NONCONVEX -5
 #define DAQP_EXIT_OVERDETERMINED_INITIAL -6
 #define DAQP_EXIT_TIMELIMIT -7
+#define DAQP_EXIT_UNSUPPORTED -8
 
 // UPDATE LDP MASKS
 #define DAQP_UPDATE_Rinv 1

include/types.h

@@ -263,6 +263,9 @@ typedef struct{
     c_float *Mu;
 }DAQPWorkspace;
 
+#define DAQP_IS_HIERARCHICAL(work) \
+    ((work)->break_points != NULL && (work)->nh > 1)
+
 # ifdef __cplusplus
 }
 # endif // ifdef __cplusplus

interfaces/daqp-julia/src/api.jl

@@ -370,7 +370,7 @@ function update(daqp::DAQPBase.Model, H,f,A,bupper,blower,sense=nothing,break_po
         update_mask |= DAQP_UPDATE_sense
     end
 
-    if(!isnothing(break_points) && length(break_points)== work.nh)
+    if(!isnothing(break_points) && length(break_points)== qp.nh)
         daqp.qpj.break_points .= break_points
         update_mask |= DAQP_UPDATE_hierarchy
     end

interfaces/daqp-julia/src/constants.jl

@@ -33,6 +33,7 @@ const  ITERLIMIT      = -4
 const  NONCONVEX      = -5
 const  OVERDETERMINED = -6
 const  TIMELIMIT      = -7
+const  UNSUPPORTED    = -8
 
 const flag2status= Dict{Int,Symbol}(3 => :Constrained_Point,
                                     2 => :Soft_Optimal,
@@ -43,4 +44,5 @@ const flag2status= Dict{Int,Symbol}(3 => :Constrained_Point,
                                    -4 => :Iteration_Limit,
                                    -5 => :Nonconvex_Problem,
                                    -6 => :Initial_Overdetermined,
-                                   -7 => :Time_Limit)
+                                   -7 => :Time_Limit,
+                                   -8 => :Unsupported_Problem)

interfaces/daqp-julia/test/core_tests.jl

@@ -362,6 +362,53 @@ end
 
 end
 
+# Only daqp_hiqp can handle a hierarchy
+@testset "Hierarchical dispatch guard" begin
+    A = [1.0 1 1; 1 -1 0; 3 1 -1]
+    bu = [ones(3);1;0.5;20]
+    bl = [-ones(3);-1e30;0.5;10]
+    sense = zeros(Cint,6)
+    Hsing = [1.0 0 0; 0 1 0; 0 0 0]
+    Hpd   = [2.0 0 0; 0 2 0; 0 0 2]
+    bps   = [3;4;5;6]
+
+    # Rejected up front, before any solve is attempted
+    d = DAQPBase.Model()
+    exitflag,_ = DAQPBase.setup(d,Hsing,zeros(3),A,bu,bl,sense;break_points = bps)
+    @test exitflag == DAQPBase.UNSUPPORTED
+    @test !d.has_model
+
+    # A hierarchy with a positive definite Hessian is unaffected
+    d = DAQPBase.Model()
+    exitflag,_ = DAQPBase.setup(d,Hpd,zeros(3),A,bu,bl,sense;break_points = bps)
+    @test exitflag >= 0
+    x,fval,exitflag,info = solve(d)
+    @test exitflag > 0
+    @test info.nodes == 1
+    @test unsafe_load(d.work).nh == 4
+
+    # An update that turns the Hessian singular must report it too. Nothing
+    # downstream re-checks, so the caller is expected to act on the flag
+    # rather than solve a workspace that was rejected.
+    @test DAQPBase.update(d,Hsing,nothing,nothing,nothing,nothing) == DAQPBase.UNSUPPORTED
+
+    # The rejection must leave nh/break_points intact -- they outlive the
+    # solve (later solves, codegen), so nothing may overwrite them on the way.
+    @test unsafe_load(d.work).nh == 4
+    srcdir = tempname()
+    DAQPBase.codegen(d,dir=srcdir)
+    @test occursin("DAQP_HIERARCHICAL", read(joinpath(srcdir,"daqp_workspace.h"),String))
+    @test occursin("break_points[4]", read(joinpath(srcdir,"daqp_workspace.c"),String))
+    rm(srcdir,recursive=true)
+
+    # A hierarchy paired with a non-symmetric AVI is rejected the same way
+    d = DAQPBase.Model()
+    Hns = [2.0 1 0; -1 2 0; 0 0 2]
+    exitflag,_ = DAQPBase.setup(d,Hns,zeros(3),A,bu,bl,sense;break_points = bps, is_avi=true)
+    @test exitflag == DAQPBase.UNSUPPORTED
+    @test !d.has_model
+end
+
 @testset "Trivial infeasible" begin
     H = [6.837677669279314 1.3993262799977795 1.9781574256330445 0.7988389688453156;
          1.3993262799977795 4.91607513347457 0.8347008717503388 0.964319980996552;
@@ -420,8 +467,9 @@ end
     # Test avi interface by passing settings
     xref,H,f,A,b = generate_test_avi(n,m)
     s = settings(DAQPBase.Model(), Dict(:iter_limit => 1))
-    _,_,exitflag,_ = DAQPBase.avi(H,f,A,b; settings=s)
+    _,_,exitflag,info = DAQPBase.avi(H,f,A,b; settings=s)
     @test exitflag == DAQPBase.ITERLIMIT
+    @test info.nodes == 1
 
     # Test that update does not cause a segfault for AVIs
     n = 10; m = 50

interfaces/daqp-matlab/daqp.m

@@ -195,7 +195,7 @@ function soften_constraints(this,ids)
             % TODO: update workspace
             this
         end
-        function update(this,H,f,A,bupper,blower,sense)
+        function exitflag = update(this,H,f,A,bupper,blower,sense)
             update_mask = int32(0);
             if(size(H,1)==this.n && size(H,2) == this.n)
                 this.H = H;
@@ -218,7 +218,7 @@ function update(this,H,f,A,bupper,blower,sense)
                 this.sense = int32(sense);
                 update_mask = update_mask+16;
             end
-            daqpmex('update', this.work_ptr,...
+            exitflag = daqpmex('update', this.work_ptr,...
                 this.H,this.f,this.A,this.bupper,this.blower,this.sense,...
                 update_mask);
         end

interfaces/daqp-matlab/daqpmex.c

@@ -215,7 +215,8 @@ void mexFunction( int nlhs, mxArray *plhs[],
 	  work->qp->sense= (int *)mxGetPr(prhs[7]);
 	  // Update LDP with new QP data
 	  const int update_mask = (int)mxGetScalar(prhs[8]);
-	  daqp_update_ldp(update_mask,work,work->qp);
+	  const int error_flag = daqp_update_ldp(update_mask,work,work->qp);
+	  plhs[0] = mxCreateDoubleScalar(error_flag);
 	}
     else if (!strcmp("codegen", cmd)) {
         char fname[64];

interfaces/daqp-python/daqp.pyx

@@ -156,7 +156,7 @@ def solve(double[:, :] H, double[:] f, double[:, :] A,
            * setup_time : Time for settings up the problem
            * solve_time : Time for solving the problem
            * iterations : Number of performed iterations
-           * nodes      : Explored nodes in branch-and-bound tree
+           * nodes      : Explored BnB nodes, or outer iterations for proximal/AVI solves
            * lam        : Optimal dual solution
     """
 

interfaces/daqp-python/setup.py

@@ -46,7 +46,7 @@
         include_dirs=[str(csrc_dir / 'include')])
 
 setup(name='daqp',
-        version='0.9.0',
+        version='0.9.1',
         description='DAQP: A dual active-set QP solver',
         url='http://github.com/darnstrom/daqp',
         author='Daniel Arnström',

interfaces/daqp-python/test/example_test.py

@@ -493,8 +493,46 @@ def test_positive_eps_forces_full_proximal(self):
         self.assertEqual(flag_ref, 1)
         self.assertEqual(flag_force, 1)
         self.assertGreater(info_force['iterations'], 1)
+        self.assertGreater(info_force['nodes'], 1)
+        self.assertLessEqual(info_force['nodes'], info_force['iterations'])
         np.testing.assert_allclose(x_force, x_ref, atol=1e-6)
 
+    def test_single_break_point_stays_non_hierarchical(self):
+        """A lone break point must not turn into a hierarchy across solves.
+
+        Fewer than two levels is no hierarchy, so installing the array must
+        drop it: break_points != NULL is what marks nh as describing a
+        hierarchy, and the proximal loop reports its outer-iteration count
+        through nh once it does not. Keeping a one-entry array here would let
+        that count be read back as a hierarchy on the next solve.
+        """
+        H = np.array([[4.0, 1.0], [1.0, 3.0]], dtype=c_double)
+        f = np.array([1.0, 2.0], dtype=c_double)
+        A = np.array([[1.0, 1.0]], dtype=c_double)
+        bupper = np.array([5.0, 5.0, 5.0], dtype=c_double)
+        blower = np.array([-5.0, -5.0, -5.0], dtype=c_double)
+        sense = np.array([0, 0, 0], dtype=c_int)
+        break_points = np.array([3], dtype=np.intc)
+
+        d = daqp.Model()
+        d.settings = {'eps_prox': 1e-2}
+        setup_flag, _ = d.setup(
+            H, f, A, bupper, blower, sense, break_points=break_points)
+        self.assertGreaterEqual(setup_flag, 0)
+        self.assertEqual(d.update(break_points=break_points), 0)
+
+        x1, _, flag1, info1 = d.solve()
+        self.assertEqual(flag1, 1)
+        self.assertGreater(info1['iterations'], 1)
+        # nh is free to count outer iterations here, so it must report them
+        self.assertGreater(info1['nodes'], 1)
+        self.assertLessEqual(info1['nodes'], info1['iterations'])
+
+        # The second solve must not see a hierarchy fabricated by the first.
+        x2, _, flag2, _ = d.solve()
+        self.assertEqual(flag2, 1)
+        np.testing.assert_allclose(x1, x2, atol=1e-9)
+
     def test_consistent_redundant_equalities_are_ignored(self):
         """Dependent consistent equalities do not make setup fail."""
         H = np.eye(2, dtype=c_double)

src/api.c

@@ -6,6 +6,7 @@
 
 // Solve problem from a given workspace and measure setup and solve time
 void daqp_solve(DAQPResult *res, DAQPWorkspace *work){
+    if(work->break_points == NULL) work->nh = 1;
     // Put back an elimination that a previous solve retrieved
     if((res->exitflag = daqp_eq_reinstall(work)) < 0) return;
 #ifdef PROFILING
@@ -19,7 +20,7 @@ void daqp_solve(DAQPResult *res, DAQPWorkspace *work){
             if(work->avi == NULL || work->avi->is_symmetric){
                 if(work->bnb != NULL)
                     res->exitflag = daqp_bnb(work);
-                else if(work->nh > 1)
+                else if(DAQP_IS_HIERARCHICAL(work))
                     res->exitflag = daqp_hiqp(work,res->lam);
                 else
                     res->exitflag = daqp_ldp(work);
@@ -329,7 +330,7 @@ void allocate_daqp_workspace(DAQPWorkspace *work, int n, int ns){
 #endif
 
     work->bnb = NULL;
-    work->nh = 0;
+    work->nh = 1;
     work->break_points = NULL;
     work->avi = NULL;
     work->eq = NULL;
@@ -459,7 +460,7 @@ void daqp_extract_result(DAQPResult* res, DAQPWorkspace* work){
     // Extract dual solution for ordinary QPs. Hierarchical QPs populate the
     // output duals in daqp_hiqp(), where they represent the soft-level
     // penalties used by the public interfaces.
-    if(res->lam != NULL && work->nh < 2){
+    if(res->lam != NULL && !DAQP_IS_HIERARCHICAL(work)){
         for(i=0;i<work->m;i++)
             res->lam[i] = 0;
         for(i=0;i<work->n_active;i++)
@@ -482,7 +483,12 @@ void daqp_extract_result(DAQPResult* res, DAQPWorkspace* work){
     // info
     res->soft_slack = work->soft_slack;
     res->iter = work->iterations;
-    res->nodes = (work->bnb == NULL) ? 1 : work->bnb->nodecount;
+    if(work->bnb != NULL)
+        res->nodes = work->bnb->nodecount;
+    else if(DAQP_IS_HIERARCHICAL(work))
+        res->nodes = 1;
+    else
+        res->nodes = work->nh;
 
     // Expand a reduced equality-eliminated result and restore the full LDP.
     daqp_eq_retrieve(res,work);

src/avi.c

@@ -15,12 +15,16 @@ int daqp_solve_avi(DAQPWorkspace *work) {
     int retry_requested = 0;
     c_float minimum_newton_residual = DAQP_INF;
 
+    // nh act as a counter for outer iterations
+    work->nh = 0;
+
     // Initial avi iterate
     for(i=0; i < work->n; i++) work->avi->x[i] = work->x[i];
 
     // Start the iterations
     // TODO iter_limit should be the for tot_iter...
     for (k = 0; k < work->settings->iter_limit; k++) {
+        work->nh++;
         // Compute xtemp = H*x + f - (Hsym + I)x
         for(i=0, disp=0; i < work->n; i++){
             sum = sum2 = 0.0;
@@ -103,13 +107,15 @@ int daqp_solve_avi(DAQPWorkspace *work) {
     }
     if(retry_requested){
         int original_limit = work->settings->iter_limit;
+        int previous_outer_iterations = work->nh;
         int retry_flag = daqp_retry_avi_with_reduced_rho(work);
         if(retry_flag < 0) return retry_flag;
         if(retry_flag > 0 && k+1 < original_limit){
             work->settings->iter_limit = original_limit-(k+1);
             retry_flag = daqp_solve_avi(work);
             work->settings->iter_limit = original_limit;
             work->iterations += tot_iter;
+            work->nh += previous_outer_iterations;
             return retry_flag;
         }
         work->iterations = tot_iter;

src/daqp_prox.c

@@ -29,6 +29,9 @@ int daqp_prox(DAQPWorkspace *work){
     c_float eta = work->settings->eta_prox;
     c_float eps;
 
+    // nh act as a counter for outer iterations
+    work->nh = 0;
+
     // The outer iteration works in the full space; the elimination is applied
     // to each inner least-distance problem below
     daqp_eq_restore(work);
@@ -122,6 +125,7 @@ int daqp_prox(DAQPWorkspace *work){
          * Solve the (regularised) least-distance problem
          * ----------------------------------------------------------------*/
         work->u = work->x;
+        work->nh++;
         exitflag = daqp_ldp(work);
 
         total_iter += work->iterations;

src/eq_elim.c

@@ -126,7 +126,7 @@ static void apply_Q(const DAQPEqElim* eq, c_float* y){
 
 static int is_eq_elim_eligible(const DAQPWorkspace* work){
     if((work->avi != NULL && !work->avi->is_symmetric) ||
-            work->bnb != NULL || work->nh > 1) return 0;
+            work->bnb != NULL || DAQP_IS_HIERARCHICAL(work)) return 0;
     if(work->qp == NULL || work->qp->A == NULL) return 0;
     /*
      * A singular Hessian is handled by the proximal method, which solves a
@@ -705,7 +705,8 @@ void daqp_eq_retrieve(DAQPResult* res, DAQPWorkspace* work){
     int i, j;
     if(eq == NULL || eq->neq == 0) return;
     if(eq->installed) expand_solution(work,work->u,eq->tmp);
-    if((eq->installed || eq->expanded) && res->lam != NULL && work->nh < 2){
+    if((eq->installed || eq->expanded) && res->lam != NULL &&
+            !DAQP_IS_HIERARCHICAL(work)){
         // Scatter the multipliers onto the constraints of the original problem
         const int m_r = eq->m_r;
         for(i = eq->m; i > m_r; ) res->lam[--i] = 0;
@@ -728,7 +729,7 @@ void daqp_eq_retrieve(DAQPResult* res, DAQPWorkspace* work){
         res->fval *= 0.5;
     }
     eq->expanded = 0;
-    if(res->lam != NULL && work->nh < 2){
+    if(res->lam != NULL && !DAQP_IS_HIERARCHICAL(work)){
         compute_lam_eq(work,res->lam);
         for(i = 0; i < eq->neq; i++) res->lam[eq->eq_ids[i]] = eq->lam_eq[i];
         for(i = 0; i < eq->ndrop; i++) res->lam[eq->drop_ids[i]] = 0;

src/hierarchical.c

@@ -9,7 +9,7 @@ int daqp_hiqp(DAQPWorkspace *work, c_float *lambda){
     int exitflag=0;
 
     // If only one hiearchy -> just solve normal LDP
-    if( work->nh < 2) return daqp_ldp(work);
+    if(!DAQP_IS_HIERARCHICAL(work)) return daqp_ldp(work);
 
     // Reset lambda for output
     if(lambda != NULL) for(i=0;i<work->m;i++) lambda[i]=0;