/* ========================================================================== */
/* === Source/Mongoose_QPBoundary.cpp ======================================= */
/* ========================================================================== */

/* -----------------------------------------------------------------------------
 * Mongoose Graph Partitioning Library  Copyright (C) 2017-2018,
 * Scott P. Kolodziej, Nuri S. Yeralan, Timothy A. Davis, William W. Hager
 * Mongoose is licensed under Version 3 of the GNU General Public License.
 * Mongoose is also available under other licenses; contact authors for details.
 * -------------------------------------------------------------------------- */

/*
   Move all components of x to boundary of the feasible region

    0 <= x <= 1, a'x = b, lo <= b <= hi

   while decreasing the cost function. The algorithm has the following parts

   1. For each i in the free set, see if x_i can be feasibly pushed to either
      boundary while decreasing the cost.

   2. For each i in the bound set, see if x_i can be feasibly flipped to
      opposite boundary while decreasing the cost.

   3. For each i in the free list with a_{ij} = 0 and with j free,
      move either x_i or x_j to the boundary while decreasing
      the cost. The adjustments has the form x_i = s/a_i and x_j = -s/a_j
      where s is a scalar factor. These adjustments must decrease cost.

   4. For the remaining i in the free list, take pair x_i and x_j and
      apply adjustments of the same form as in #2 above to push at least one
      component to boundary. The quadratic terms can only decrease the
      cost function. We choose the sign of s such that g_i x_i + g_j x_j <= 0.
      Hence, this adjustment cannot increase the cost.
 */

/* ========================================================================== */

#include "Mongoose_QPBoundary.hpp"
#include "Mongoose_Debug.hpp"
#include "Mongoose_Internal.hpp"

#define EMPTY (-1)

namespace Mongoose
{

void QPBoundary(EdgeCutProblem *graph, const EdgeCut_Options *options, QPDelta *QP)
{
    (void)options; // Unused variable
    /* ---------------------------------------------------------------------- */
    /* Step 0. read in the needed arrays                                      */
    /* ---------------------------------------------------------------------- */

    /* input and output */

    //--- FreeSet
    Int nFreeSet        = QP->nFreeSet;
    Int *FreeSet_list   = QP->FreeSet_list; /* list for free indices */
    Int *FreeSet_status = QP->FreeSet_status;
    /* FreeSet_status [i] = +1, -1, or 0
       if x_i = 1, 0, or 0 < x_i < 1 */
    //---

    PR(("Mongoose_QPBoundary nFreeSet %ld\n", nFreeSet));

    if (nFreeSet == 0)
    {
        // quick return if FreeSet is empty
        return;
    }

    double *x    = QP->x;        /* current estimate of solution */
    double *grad = QP->gradient; /* gradient at current x */
    Int ib       = QP->ib;       /* ib = +1, -1, or 0 ,
         if b = hi, lo, or lo < b < hi, respectively.  Note there are cases
         where roundoff occurs, and ib can be zero even though b == lo or
         b == hi.  The value of be can even be < lo or > hi, but only by a tiny
         amount of roundoff error.  This is OK. */

    double b = QP->b; /* current value for a'x */

    /* problem specification for the graph G */
    Int n      = graph->n; /* problem dimension */
    double *Ex = graph->x; /* numerical values for edge weights */
    Int *Ei    = graph->i; /* adjacent vertices for each vertex */
    Int *Ep    = graph->p; /* points into Ex or Ei */
    double *a  = graph->w; /* a'x = b, lo <= b <= hi */

    double lo = QP->lo;
    double hi = QP->hi;

    /* work array */
    double *D = QP->D; /* diagonal of quadratic */

    PR(("\n----- QPBoundary start: [\n"));
    DEBUG(QPcheckCom(graph, options, QP, 1, QP->nFreeSet, QP->b)); // check b

    /* ---------------------------------------------------------------------- */
    /* Step 1. if lo < b < hi, then for each free k,                          */
    /*         see if x_k can be pushed to 0 or 1                             */
    /* ---------------------------------------------------------------------- */

    DEBUG(FreeSet_dump("QPBoundary start", n, FreeSet_list, nFreeSet,
                       FreeSet_status, 0, x));

    PR(("Boundary 1 start: ib %ld lo %g b %g hi %g b-lo %g hi-b %g\n", ib, lo,
        b, hi, b - lo, hi - b));

    Int kfree2 = 0;
    for (Int kfree = 0; kfree < nFreeSet; kfree++)
    {
        // Once b becomes bounded, the remainder of the FreeSet is unchanged,
        // and no further changes are made to x.  However, this loop must still
        // continue, so as to compact the FreeSet from deletions made by earlier
        // iterations.

        // get the next k from the FreeSet
        Int k = FreeSet_list[kfree];

        PR(("Step 1: k %ld  x[k] %g  ib %ld b %g\n", k, x[k], ib, b));

        // only modify x[k] if ib == 0 (which means lo < b < hi)
        if (ib == 0)
        {
            double delta_xk;
            double ak = (a) ? a[k] : 1;
            if (grad[k] > 0.0)
            {
                // decrease x [k]
                delta_xk = (b - lo) / ak; // note that delta_xk > 0
                if (delta_xk < x[k])
                {
                    // x [k] decreases by delta_xk but does not hit zero
                    // b hits the lower bound, lo
                    ib = -1;
                    b  = lo;
                    x[k] -= delta_xk;
                    //--- keep k in the FreeSet
                    FreeSet_list[kfree2++] = k;
                }
                else
                {
                    // x [k] hits lower bound of zero
                    // b does not hit lo; still between lower and upper bound
                    delta_xk          = x[k];
                    x[k]              = 0.;
                    FreeSet_status[k] = -1;
                    b -= delta_xk * ak;
                    //--- remove k from the FreeSet by not incrementing kfree2
                }
            }
            else
            {
                // increase x [k]
                delta_xk = (b - hi) / ak; // note that delta_xk < 0
                if (delta_xk < x[k] - 1.)
                {
                    // x [k] hits upper bound of one
                    // b does not reach hi; still between lower and upper bound
                    delta_xk          = x[k] - 1.;
                    x[k]              = 1.;
                    FreeSet_status[k] = +1;
                    b -= delta_xk * ak;
                    //--- remove k from the FreeSet by not incrementing kfree2
                }
                else
                {
                    // x [k] increases by -delta_xk but does not hit one
                    // b hits the upper bound, hi.
                    ib = +1;
                    b  = hi;
                    x[k] -= delta_xk;
                    //--- keep k in the FreeSet
                    FreeSet_list[kfree2++] = k;
                }
            }
            // x [k] has dropped by delta_xk, so update the gradient
            for (Int p = Ep[k]; p < Ep[k + 1]; p++)
            {
                grad[Ei[p]] += delta_xk * ((Ex) ? Ex[p] : 1);
            }
            grad[k] += delta_xk * D[k];
        }
        else
        {
            // b is at lo or hi and thus x [k] is not changed.
            // Once this happens, the remainder of this loop does this next
            // step only, and no further changes are made to x and the FreeSet.
            //--- keep k in the FreeSet
            FreeSet_list[kfree2++] = k;
        }
    }

    // update the size of the FreeSet, after pruning
    nFreeSet = kfree2;

    /* ---------------------------------------------------------------------- */
    /* Step 2. Examine flips of x_k from 0 to 1 or from 1 to 0 */
    /* ---------------------------------------------------------------------- */

    PR(("Boundary step 2:\n"));

    for (Int k = 0; k < n; k++)
    {
        Int FreeSet_status_k = FreeSet_status[k];
        if (FreeSet_status_k == 0)
        {
            // k is in FreeSet so it cannot be simply flipped 0->1 or 1->0
            continue;
        }

        // k not in FreeSet, so no changes here to FreeSet

        double ak = (a) ? a[k] : 1;
        if (FreeSet_status_k > 0) /* try changing x_k from 1 to 0 */
        {
            if (b - ak >= lo)
            {
                if (0.5 * D[k] + grad[k] >= 0) /* flip lowers cost */
                {
                    b -= ak;
                    ib                = (b <= lo ? -1 : 0);
                    x[k]              = 0.0;
                    FreeSet_status[k] = -1;
                }
            }
        }
        else /* try changing x_k from 0 to 1 */
        {
            if (b + ak <= hi)
            {
                if (grad[k] - 0.5 * D[k] <= 0) /* flip lowers cost */
                {
                    b += ak;
                    ib                = (b >= hi ? 1 : 0);
                    x[k]              = 1.0;
                    FreeSet_status[k] = +1;
                }
            }
        }

        if (FreeSet_status_k != FreeSet_status[k])
        {
            if (FreeSet_status_k == 1) /* x [k] was 1, now it is 0 */
            {
                for (Int p = Ep[k]; p < Ep[k + 1]; p++)
                {
                    grad[Ei[p]] += (Ex) ? Ex[p] : 1;
                }
                grad[k] += D[k];
            }
            else /* x [k] was 0, now it is 1 */
            {
                for (Int p = Ep[k]; p < Ep[k + 1]; p++)
                {
                    grad[Ei[p]] -= (Ex) ? Ex[p] : 1;
                }
                grad[k] -= D[k];
            }
        }
        // DEBUG (QPcheckCom (graph, options, QP, 1, nFreeSet, b)) ;         //
        // check b
    }

    /* ---------------------------------------------------------------------- */
    // quick return if FreeSet is now empty
    /* ---------------------------------------------------------------------- */

    if (nFreeSet == 0)
    {
        PR(("Boundary quick: ib %ld lo %g b %g hi %g b-lo %g hi-b %g\n", ib, lo,
            b, hi, b - lo, hi - b));
        QP->nFreeSet = nFreeSet;
        QP->b        = b;
        QP->ib       = ib;
        PR(("------- QPBoundary end ]\n"));
        return;
    }

    /* ---------------------------------------------------------------------- */
    /* Step 3. Search for a_{ij} = 0 in the free index set */
    /* ---------------------------------------------------------------------- */

    // look for where both i and j are in the FreeSet,
    // but i and j are not adjacent in the graph G.

    DEBUG(FreeSet_dump("step 3", n, FreeSet_list, nFreeSet, FreeSet_status, 0,
                       x));

    // for each j in FreeSet, except for the last one
    for (Int jfree = 0; jfree < nFreeSet - 1; jfree++)
    {

        // get j from the FreeSet
        Int j = FreeSet_list[jfree];
        if (j == EMPTY)
        {
            // j has already been deleted, skip it
            continue;
        }

        /* -------------------------------------------------------------- */
        /* find i and j both free and where a_{ij} = 0 */
        /* -------------------------------------------------------------- */

        // mark all vertices i adjacent to j in the FreeSet
        for (Int p = Ep[j]; p < Ep[j + 1]; p++)
        {
            Int i = Ei[p];
            ASSERT(i != j); // graph has no self edges
            graph->mark(i);
        }
        graph->mark(j);

        // for each i that follows after j in the FreeSet
        for (Int ifree = jfree + 1; ifree < nFreeSet; ifree++)
        {

            // get i from the FreeSet
            Int i = FreeSet_list[ifree];
            if (i == EMPTY)
            {
                // i has already been deleted it; skip it
                continue;
            }

            if (!graph->isMarked(i))
            {
                // vertex i is not adjacent to j in the graph G
                double aj = (a) ? a[j] : 1;
                double ai = (a) ? a[i] : 1;
                double xi = x[i];
                double xj = x[j];

                /* cost change if x_j increases dx_j = s/a_j, dx_i = s/a_i */
                double s;
                Int bind1, bind2;
                if (aj * (1. - xj) < ai * xi) // x_j hits upper bound
                {
                    s     = aj * (1. - xj);
                    bind1 = 1;
                }
                else /* x_i hits lower bound */
                {
                    s     = ai * xi;
                    bind1 = 0;
                }
                double dxj = s / aj;
                double dxi = -s / ai;
                double c1  = (grad[j] - .5 * D[j] * dxj) * dxj
                            + (grad[i] - .5 * D[i] * dxi) * dxi;

                /* cost change if x_j decreases dx_j = s/a_j, dx_i = s/a_i */
                if (aj * xj < ai * (1. - xi)) // x_j hits lower bound
                {
                    s     = -aj * xj;
                    bind2 = -1;
                }
                else /* x_i hits upper bound */
                {
                    s     = -ai * (1. - xi);
                    bind2 = 0;
                }
                dxj       = s / aj;
                dxi       = -s / ai;
                double c2 = (grad[j] - 0.5 * D[j] * dxj) * dxj
                            + (grad[i] - 0.5 * D[i] * dxi) * dxi;

                Int new_FreeSet_status;
                if (c1 < c2) /* increase x_j */
                {
                    if (bind1 == 1)
                    {
                        // j is bound (not i) and x_j becomes 1
                        dxj  = 1. - xj;
                        dxi  = -aj * dxj / ai;
                        x[j] = 1.;
                        x[i] += dxi;
                        new_FreeSet_status = +1; /* j is bound at 1 */
                    }
                    else // bind1 is zero
                    {
                        // i is bound (not j) and x_i becomes 0
                        dxi  = -xi;
                        dxj  = -ai * dxi / aj;
                        x[i] = 0.;
                        x[j] += dxj;
                        new_FreeSet_status = -1; /* i is bound at 0 */
                    }
                }
                else
                {
                    if (bind2 == -1)
                    {
                        // j is bound (not i) and x_j becomes 0
                        bind1 = 1;
                        x[j]  = 0.;
                        x[i] += dxi;
                        new_FreeSet_status = -1; /* j is bound at 0 */
                    }
                    else /* x_i = 1 */
                    {
                        // i is bound (not j) and x_i becomes 1
                        bind1 = 0;
                        x[i]  = 1;
                        x[j] += dxj;
                        new_FreeSet_status = +1; /* i is bound at 1 */
                    }
                }

                for (Int p = Ep[j]; p < Ep[j + 1]; p++)
                {
                    grad[Ei[p]] -= ((Ex) ? Ex[p] : 1) * dxj;
                }
                for (Int p = Ep[i]; p < Ep[i + 1]; p++)
                {
                    grad[Ei[p]] -= ((Ex) ? Ex[p] : 1) * dxi;
                }
                grad[j] -= D[j] * dxj;
                grad[i] -= D[i] * dxi;

                // Remove either i or j from the FreeSet.  Note that it
                // is possible for both x[i] and x[j] to reach their bounds
                // at the same time.  Only one is removed from the FreeSet;
                // the other will be removed later.

                if (bind1)
                {
                    // remove j from the FreeSet by setting its place to EMPTY
                    PR(("(b1):remove j = %ld from the FreeSet\n", j));
                    ASSERT(j == FreeSet_list[jfree]);
                    ASSERT(FreeSet_status[j] == 0);
                    FreeSet_list[jfree] = EMPTY;
                    FreeSet_status[j]   = new_FreeSet_status;
                    ASSERT(FreeSet_status[j] != 0);
                    //---
                    // no longer consider j, so skip all of remainder of i loop
                    break;
                }
                else
                {
                    // remove i from the FreeSet by setting its place to EMPTY
                    PR(("(b2):remove i = %ld from the FreeSet\n", i));
                    ASSERT(i == FreeSet_list[ifree]);
                    ASSERT(FreeSet_status[i] == 0);
                    FreeSet_list[ifree] = EMPTY;
                    FreeSet_status[i]   = new_FreeSet_status;
                    ASSERT(FreeSet_status[i] != 0);
                    //---
                    // keep j, and consider it with the next i
                    continue;
                }
            }
        }

        // clear the marks from all the vertices
        graph->clearMarkArray();
    }

    // remove deleted vertices from the FreeSet
    kfree2 = 0;
    for (Int kfree = 0; kfree < nFreeSet; kfree++)
    {
        Int k = FreeSet_list[kfree];
        if (k != EMPTY)
        {
            // keep k in the FreeSet
            FreeSet_list[kfree2++] = k;
            ASSERT(0 <= k && k < n);
            ASSERT(FreeSet_status[k] == 0);
        }
    }
    nFreeSet = kfree2;

    DEBUG(FreeSet_dump("step 3 done", n, FreeSet_list, nFreeSet, FreeSet_status,
                       1, x));

    DEBUG(QPcheckCom(graph, options, QP, 1, nFreeSet, b)); // check b

#ifndef NDEBUG
    // the vertices in the FreeSet now form a single clique.  Check this.
    // this test is for debug mode only
    ASSERT(nFreeSet >= 1); // we can have 1 or more vertices still in FreeSet
    for (Int kfree = 0; kfree < nFreeSet; kfree++)
    {
        // j must be adjacent to all other vertices in the FreeSet
        Int j               = FreeSet_list[kfree];
        Int nfree_neighbors = 0;
        for (Int p = Ep[j]; p < Ep[j + 1]; p++)
        {
            Int i = Ei[p];
            ASSERT(i != j);
            if (FreeSet_status[i] == 0)
                nfree_neighbors++;
        }
        ASSERT(nfree_neighbors == nFreeSet - 1);
    }
#endif

    /* ---------------------------------------------------------------------- */
    /* Step 4. dxj = s/aj, dxi = -s/ai, choose s with g_j dxj + g_i dxi <= 0 */
    /* ---------------------------------------------------------------------- */

    DEBUG(FreeSet_dump("step 4 starts", n, FreeSet_list, nFreeSet,
                       FreeSet_status, 0, x));

    // consider pairs of vertices in the FreeSet, until only one is left
    while (nFreeSet > 1)
    {
        /* free variables: 0 < x_j < 1 */
        /* choose s so that first derivative terms decrease */

        // i and j are the last two vertex in the FreeSet_list, as in:
        // FreeSet_list = [ .... i j ]
        // at the end of this iteration, one will be deleted, thus becoming
        // FreeSet_list = [ .... j ]
        // or
        // FreeSet_list = [ .... i ]
        Int j = FreeSet_list[nFreeSet - 1];
        ASSERT(FreeSet_status[j] == 0);
        Int i = FreeSet_list[nFreeSet - 2];
        ASSERT(FreeSet_status[i] == 0);

        double ai = (a) ? a[i] : 1;
        double aj = (a) ? a[j] : 1;
        double xi = x[i];
        double xj = x[j];

        Int new_FreeSet_status;
        Int bind1;
        double dxj, dxi, s = grad[j] / aj - grad[i] / ai;

        if (s < 0.) /* increase x_j */
        {
            if (aj * (1. - xj) < ai * xi) /* x_j hits upper bound */
            {
                dxj  = 1. - xj;
                dxi  = -aj * dxj / ai;
                x[j] = 1.;
                x[i] += dxi;
                new_FreeSet_status = +1;
                bind1              = 1; /* x_j is bound at 1 */
            }
            else /* x_i hits lower bound */
            {
                dxi  = -xi;
                dxj  = -ai * dxi / aj;
                x[i] = 0.;
                x[j] += dxj;
                new_FreeSet_status = -1;
                bind1              = 0; /* x_i is bound at 0 */
            }
        }
        else /* decrease x_j */
        {
            if (aj * xj < ai * (1. - xi)) /* x_j hits lower bound */
            {
                dxj  = -xj;
                dxi  = -aj * dxj / ai;
                x[j] = 0;
                x[i] += dxi;
                new_FreeSet_status = -1;
                bind1              = 1; /* x_j is bound */
            }
            else /* x_i hits upper bound */
            {
                dxi  = 1 - xi;
                dxj  = -ai * dxi / aj;
                x[i] = 1;
                x[j] += dxj;
                new_FreeSet_status = +1;
                bind1              = 0; /* x_i is bound */
            }
        }

        for (Int k = Ep[j]; k < Ep[j + 1]; k++)
        {
            grad[Ei[k]] -= ((Ex) ? Ex[k] : 1) * dxj;
        }
        for (Int k = Ep[i]; k < Ep[i + 1]; k++)
        {
            grad[Ei[k]] -= ((Ex) ? Ex[k] : 1) * dxi;
        }
        grad[j] -= D[j] * dxj;
        grad[i] -= D[i] * dxi;

        // ---------------------------------------------------------------------
        // the following 2 cases define the next j in the iteration:
        // ---------------------------------------------------------------------

        // Remove either i or j from the FreeSet.  Note that it is possible for
        // both x[i] and x[j] to reach their bounds at the same time.  Only one
        // is removed from the FreeSet; the other will be removed later.

        if (bind1)
        {
            // j is bound.
            // remove j from the FreeSet, and keep i.  The FreeSet_list was
            // FreeSet_list = [ .... i j ] becomes FreeSet_list = [ .... i ]
            PR(("(b3):remove j = %ld from the FreeSet\n", j));
            ASSERT(FreeSet_status[j] == 0);
            FreeSet_status[j] = new_FreeSet_status;
            ASSERT(FreeSet_status[j] != 0);
        }
        else
        {
            // i is bound.
            // remove i from the FreeSet, and keep j.  The FreeSet_list was
            // FreeSet_list = [ .... i j ] becomes FreeSet_list = [ .... j ]
            PR(("(b4):remove i = %ld from the FreeSet\n", i));
            ASSERT(FreeSet_status[i] == 0);
            FreeSet_status[i] = new_FreeSet_status;
            ASSERT(FreeSet_status[i] != 0);
            // shift j down by one in the list, thus discarding j.
            ASSERT(FreeSet_list[nFreeSet - 2] == i);
            FreeSet_list[nFreeSet - 2] = j;
        }

        // one fewer vertex in the FreeSet (i or j removed)
        nFreeSet--;

        DEBUG(FreeSet_dump("step 4", n, FreeSet_list, nFreeSet, FreeSet_status,
                           0, x));
        DEBUG(QPcheckCom(graph, options, QP, 1, nFreeSet, b)); // check b
    }

    DEBUG(FreeSet_dump("wrapup", n, FreeSet_list, nFreeSet, FreeSet_status, 0,
                       x));

    /* ---------------------------------------------------------------------- */
    /* step 5: at most one free variable remaining */
    /* ---------------------------------------------------------------------- */

    ASSERT(nFreeSet == 0 || nFreeSet == 1);

    PR(("Step 5: ib %ld lo %g b %g hi %g b-lo %g hi-b %g\n", ib, lo, b, hi,
        b - lo, hi - b));

    if (nFreeSet == 1) /* j is free, optimize over x [j] */
    {
        // j is the first and only item in the FreeSet
        Int j = FreeSet_list[0];
        PR(("ONE AND ONLY!! j = %ld x[j] %g\n", j, x[j]));

        Int bind1  = 0;
        double aj  = (a) ? a[j] : 1;
        double dxj = (hi - b) / aj;
        PR(("dxj %g  x[j] %g  (1-x[j]): %g\n", dxj, x[j], 1 - x[j]));
        if (dxj < 1. - x[j])
        {
            bind1 = 1;
        }
        else
        {
            dxj = 1. - x[j];
        }

        Int bind2  = 0;
        double dxi = (lo - b) / aj;
        PR(("dxi %g  x[j] %g  (-x[j]): %g\n", dxi, x[j], -x[j]));
        if (dxi > -x[j])
        {
            bind2 = 1;
        }
        else
        {
            dxi = -x[j];
        }

        double c1 = (grad[j] - 0.5 * D[j] * dxj) * dxj;
        double c2 = (grad[j] - 0.5 * D[j] * dxi) * dxi;
        if (c1 <= c2) /* x [j] += dxj */
        {
            if (bind1)
            {
                PR(("bind1: xj changes from %g", x[j]));
                x[j] += dxj;
                PR((" to %g, b now at hi\n", x[j]));
                ib = +1;
                b  = hi;
            }
            else
            {
                x[j] = 1.;
                b += dxj * aj;
                /// remove j from the FreeSet, which is now empty
                PR(("(b5):remove j = %ld from FreeSet, now empty\n", j));
                ASSERT(FreeSet_status[j] == 0);
                FreeSet_status[j] = 1;
                ASSERT(FreeSet_status[j] != 0);
                nFreeSet--;
                ASSERT(nFreeSet == 0);
            }
        }
        else /* x [j] += dxi */
        {
            dxj = dxi;
            if (bind2)
            {
                PR(("bind2: xj changes from %g", x[j]));
                x[j] += dxj;
                PR((" to %g, b now at lo\n", x[j]));
                ib = -1;
                b  = lo;
            }
            else
            {
                x[j] = 0.;
                b += dxj * aj;
                /// remove j from the FreeSet, which is now empty
                PR(("(b6):remove j = %ld from FreeSet, now empty\n", j));
                ASSERT(FreeSet_status[j] == 0);
                FreeSet_status[j] = -1;
                ASSERT(FreeSet_status[j] != 0);
                nFreeSet--;
                ASSERT(nFreeSet == 0);
            }
        }

        if (dxj != 0.)
        {
            for (Int p = Ep[j]; p < Ep[j + 1]; p++)
            {
                grad[Ei[p]] -= ((Ex) ? Ex[p] : 1) * dxj;
            }
            grad[j] -= D[j] * dxj;
        }
    }

    /* ---------------------------------------------------------------------- */
    // wrapup
    /* ---------------------------------------------------------------------- */

    PR(("QBboundary, done:\n"));
    DEBUG(FreeSet_dump("QPBoundary: done ", n, FreeSet_list, nFreeSet,
                       FreeSet_status, 0, x));
    ASSERT(nFreeSet == 0 || nFreeSet == 1);
    PR(("Boundary done: ib %ld lo %g b %g hi %g b-lo %g hi-b %g\n", ib, lo, b,
        hi, b - lo, hi - b));

    QP->nFreeSet = nFreeSet;
    QP->b        = b;
    QP->ib       = ib;

    // clear the marks from all the vertices
    graph->clearMarkArray();

    DEBUG(QPcheckCom(graph, options, QP, 1, nFreeSet, b)); // check b
    PR(("----- QPBoundary end ]\n"));
}

} // end namespace Mongoose