graph_m4ri.cpp

Go to the documentation of this file.

#include <frechet/reachability/graph_m4ri.h>
#include <concurrency.h>
//  m4ri extension to Boolean Algebra
#if defined(__cplusplus) && ! defined(_MSC_VER)
extern "C" {
#endif
# include <m4ri/m4ri.h>
# include <m4ri/mzd_bool.h>
# include <m4ri/russian_bool.h>
#if defined(__cplusplus) && !defined(_MSC_VER)
}
#endif

#include <iomanip>
#include <boost/cstdint.hpp>
#include <boost/multiprecision/integer.hpp>

using namespace frechet;
using namespace reach;
using namespace app;

Graph::Graph(const GraphModel::ptr amodel)
    : model(amodel),
      origin{Origin::RG,-1,0,-1,-1,nullptr,nullptr,nullptr},
      diagonalElement(-1)
{
    for(Orientation o1=HORIZONTAL; o1 <= VERTICAL; ++o1)
        mask[o1] = IndexRange(o1, 0, model->count2(o1));
    for(Orientation o1=HORIZONTAL; o1 <= VERTICAL; ++o1)
        for(Orientation o2=HORIZONTAL; o2 <= VERTICAL; ++o2)
            mtx[o1][o2] = nullptr;
}

Graph::Graph(const GraphModel::ptr amodel,
      IndexRange range)
      : Graph(amodel)
{
    mask[HORIZONTAL] = range;

    //Q_ASSERT(mask[HORIZONTAL].lower < model->count1(HORIZONTAL));
    Q_ASSERT(mask[HORIZONTAL].len() <= model->count2(HORIZONTAL));
}

Graph::Graph(const GraphModel::ptr amodel,
                            Structure& str)
    : Graph(amodel)
{
    mask[HORIZONTAL].ori = HORIZONTAL;
    mask[HORIZONTAL].lower = model->map_lower(HORIZONTAL,str.bottom().first()->lower());
    mask[HORIZONTAL].upper = model->map_lower(HORIZONTAL,str.bottom().last()->upper());
    //  Note: upper bound is exclusive.
    //  Don't overlap with next RG.
    //Q_ASSERT(mask[HORIZONTAL].lower < model->count1(HORIZONTAL));
    Q_ASSERT(mask[HORIZONTAL].len() <= model->count2(HORIZONTAL));

    allocate(HORIZONTAL,VERTICAL);
    allocate(VERTICAL,HORIZONTAL);
    allocate(VERTICAL,VERTICAL);

    //  insert edges into the graph (EXCEPT H-H egdes!)
    copy(str,VERTICAL);
    copy(str,HORIZONTAL);
    Q_ASSERT(is_upper_triangular(VERTICAL,VERTICAL));
}

Graph::Graph(const Graph& that)
    : Graph(that.model)
{
    copy(that);
}

Graph::Graph(Graph&& that)
    : Graph(that.model)
{
    swap(that);
}

void Graph::allocateAll()  {
    for(Orientation o1=HORIZONTAL; o1 <= VERTICAL; ++o1)
        for(Orientation o2=HORIZONTAL; o2 <= VERTICAL; ++o2)
            allocate(o1,o2);
}

mzd_t * Graph::allocate(Orientation o1, Orientation o2) {
    Q_ASSERT(mtx[o1][o2]==nullptr);
    mtx[o1][o2] = mzd_init(mask[o1].len(),mask[o2].len());

    Q_ASSERT(mtx[o1][o2]->nrows==mask[o1].len());
    Q_ASSERT(mtx[o1][o2]->ncols==mask[o2].len());
    return mtx[o1][o2];
}

frechet::Rect Graph::rect(Orientation o1, Orientation o2) const {
    return Rect(mask[o1].lower,mask[o2].lower,
                mask[o1].upper,mask[o2].upper);
}

Graph::~Graph()
{
    release();
}

void Graph::setOriginPlacement(int di, int dj)
{
    origin.operation=Origin::PLACEMENT;
    origin.blevel = -1;
    origin.succesors = 0;
    origin.i = di;
    origin.j = dj;
    origin.A.reset();
    origin.B.reset();
    origin.P.reset();
}

void Graph::setOriginRG(int i)
{
    origin.operation=Origin::RG;
    origin.blevel = 0;
    origin.succesors = 0;
    origin.i = i;
    origin.j = i+1;
    origin.A.reset();
    origin.B.reset();
    origin.P.reset();
}

void Graph::setOriginMerge2(Graph::ptr A, Graph::ptr B)
{
    origin.operation=Origin::MERGE2;
    origin.blevel = std::max(A->origin.blevel, B->origin.blevel) + 1;
    origin.succesors = 0;
    origin.i = A->origin.i;
    origin.j = A->origin.j + (B->origin.j - B->origin.i);
    origin.A = A;
    origin.B = B;
    origin.P.reset();
    A->origin.succesors++; // count dependencies on A_VV
    B->origin.succesors++; // count dependencies on B_VV
}

void Graph::setOriginMerge3(Graph::ptr A, Graph::ptr B)
{
    origin.operation = Origin::MERGE3;
    Q_ASSERT(A->origin.blevel == 0);
    origin.blevel = B->origin.blevel + 1;
    origin.succesors = 0;
    origin.i = A->origin.i;
    origin.j = A->origin.i;
    origin.A = A;
    origin.B = B;
    origin.P.reset();
    //A->origin.succesors++; don't count dependecies on A_HV and A_VH
    B->origin.succesors++;  // count dependencies on B_VV
}

void Graph::setOriginCombine(Graph::ptr P)
{
    origin.P = P;
}

void Graph::printOrigin(std::ostream& out) const
{
    if (origin.P) 
        out << "COMBINE(";  
    switch(origin.operation)
    {
    case Origin::RG:
        out << "RG("<<origin.i<<","<<origin.j<<")";
        break;
    case Origin::MERGE2:
        Q_ASSERT(origin.A.get() != this);
        Q_ASSERT(origin.B.get() != this);
        //  otherwise we have infinit recursion right here
        out << "MERGE(";
        origin.A->printOrigin(out);
        out << ",";
        origin.B->printOrigin(out);
        out << ")";
        break;
    case Origin::MERGE3:
        Q_ASSERT(origin.A.get() != this);
        Q_ASSERT(origin.B.get() != this);
        //  otherwise we have infinit recursion right here
        out << "MERGE(";
        origin.A->printOrigin(out);
        out << ",";
        origin.B->printOrigin(out);
        out << ",";
        origin.A->printOrigin(out);
        out << ")";
        break;    default:
        Q_ASSERT(false);
    }
    if (origin.P)
        out << ")"; //  COMBINE
}

void Graph::finalize() {
    releaseIfZero();
}

void Graph::release()
{
    for(Orientation o1=HORIZONTAL; o1 <= VERTICAL; ++o1)
        for(Orientation o2=HORIZONTAL; o2 <= VERTICAL; ++o2)
            release(o1,o2);
}

void Graph::releaseIfZero()
{
    for(Orientation o1=HORIZONTAL; o1 <= VERTICAL; ++o1)
        for(Orientation o2=HORIZONTAL; o2 <= VERTICAL; ++o2)
            if (zero(o1,o2))
                release(o1,o2);
}

void Graph::release(Orientation o1, Orientation o2)
{
    if (mtx[o1][o2]) {
        mzd_free(mtx[o1][o2]);
        mtx[o1][o2] = nullptr;
    }
}


Graph& Graph::operator= (const Graph& that)
{
    Q_ASSERT(model == that.model);
    copy(that);
    return *this;
}

Graph& Graph::operator= (Graph&& that)
{
    swap(that);
    return *this;
}

bool Graph::operator== (const Graph& that) const
{
    for(Orientation o1=HORIZONTAL; o1 <= VERTICAL; ++o1)
    {
        if (mask[o1] != that.mask[o1])
            return false;

        for(Orientation o2=HORIZONTAL; o2 <= VERTICAL; ++o2)
        {
            mzd_t* m1 = mtx[o1][o2];
            mzd_t* m2 = that.mtx[o1][o2];
            if (!mzd_equal(m1,m2))
                return false;
        }
    }
    return true;
}

void Graph::copy(const Graph& that)
{
    for(Orientation o1=HORIZONTAL; o1 <= VERTICAL; ++o1)
    {
        this->mask[o1] = that.mask[o1];
        for (Orientation o2 = HORIZONTAL; o2 <= VERTICAL; ++o2)
        {
            mzd_t* m2 = that.mtx[o1][o2];
            mzd_t*& m1 = mtx[o1][o2];

            if (m1 && m2
                    && (m1->nrows == m2->nrows)
                    && (m1->ncols == m2->ncols))
            {
                mzd_copy(m1, m2);
            }
            else
            {
                if (m1)
                    mzd_free(m1);
                if (m2)
                    m1 = mzd_copy(nullptr,m2);
                else
                    m1 = nullptr;
            }
        }
    }

    this->origin = that.origin;
    this->diagonalElement = that.diagonalElement;
}

void Graph::swap(Graph& that)
{
    Q_ASSERT(&model == &that.model);
    std::swap(this->mtx, that.mtx);
    std::swap(this->origin, that.origin);
    std::swap(this->mask, that.mask);
    std::swap(this->diagonalElement, that.diagonalElement);
}

void Graph::copy(const Structure& str, Orientation ori)
{
    /*  Copy from Reachability Structure
     */
    Pointer i1 = str.first(ori).first();
    Pointer i2 = str.second(ori).first();
    for( ; i1; i1=i1->next(), i2=i2->next())
    {
        if (i1->type==NON_ACCESSIBLE) continue;
        //  traverse all reachable intervals
        IndexRange r1 = model->map(i1);

        //  map [l..h] to node indexes
        StructureIterator lh(str,*i1,i2);
        for( ; lh; ++lh) {
            if (i1->ori==HORIZONTAL && lh->ori==HORIZONTAL)
                continue;   //  ignore H-H
            IndexRange r2 = model->map((Pointer)lh);
            Q_ASSERT(lh->ori==r2.ori);
            if (lh->type != NON_ACCESSIBLE)
                add_range(r1,r2);
        }
    }
}


Graph::ptr Graph::merge2(const Graph* B, MatrixPool* pool) const
{
    Graph::ptr C = newGraph(model, model->mergedHMask(this->hmask(),B->hmask()));
    //  new Graph without data
    C->merge2(this,B,pool);
    return C;
}

void Graph::merge2(const Graph* A, const Graph* B, MatrixPool* pool)
{
    Q_ASSERT(this->mtx[VERTICAL][VERTICAL]==nullptr);
    Q_ASSERT(is_adjacent_to(*A, *B));

    mzd_t* A_VV = A->mtx[VERTICAL][VERTICAL];
    mzd_t* B_VV = B->mtx[VERTICAL][VERTICAL];
    if (! A_VV || ! B_VV)
        return;

    mzd_t* C_VV = new_mzd(A_VV->nrows,B_VV->ncols,pool);

    //  both matrices are upper triangular
    int k, blocksize;
    int threads = ConcurrencyContext::countThreads();
    mzd_optimal_parameters(A_VV, B_VV,
                           ConcurrencyContext::cacheSize(2), // right ?
                           ConcurrencyContext::cacheSize(3), // right ?
                           &k, &blocksize);

    mzd_bool_mul_uptri(C_VV, A_VV, B_VV, k,blocksize,threads);

    if (! mzd_is_zero(C_VV)) {
        this->mtx[VERTICAL][VERTICAL] = C_VV;
        Q_ASSERT(C_VV->nrows == this->mask[VERTICAL].len());
        Q_ASSERT(C_VV->ncols == this->mask[VERTICAL].len());
    }
    else {
        reclaim(C_VV,pool);
    }
}


Graph::ptr Graph::merge3(const Graph* B, MatrixPool* pool) const
{
    Graph::ptr C = newGraph(model,hmask());
    C->merge3(this,B,pool);
    return C;
}

void Graph::merge3(const Graph* A, const Graph* B, MatrixPool* pool)
{
    Q_ASSERT(this->mtx[HORIZONTAL][HORIZONTAL]==nullptr);
    Q_ASSERT(is_adjacent_to(*A, *B));
    Q_ASSERT(is_adjacent_to(*B, *A));

    mzd_t* A_HV = A->mtx[HORIZONTAL][VERTICAL];
    mzd_t* B_VV = B->mtx[VERTICAL][VERTICAL];
    mzd_t* A_VH = A->mtx[VERTICAL][HORIZONTAL];

    if (! A_HV || ! B_VV || ! A_VH)
        return;   //  shortcut for empty

    int k, blocksize;
    int threads = ConcurrencyContext::countThreads();
    mzd_optimal_parameters(A_HV, B_VV,
                           ConcurrencyContext::cacheSize(2), // right ?
                           ConcurrencyContext::cacheSize(3), // right ?
                           &k, &blocksize);

    mzd_t* temp = new_mzd(A_HV->nrows,B_VV->ncols,nullptr); // don't pool odd-sized matrices
    mzd_bool_mul_m4rm(temp, A_HV, B_VV, k,blocksize,threads);

    mzd_optimal_parameters(temp, A_VH,
                           ConcurrencyContext::cacheSize(2), // right ?
                           ConcurrencyContext::cacheSize(3), // right ?
                           &k, &blocksize);

    mzd_t* C_HH = new_mzd(A_HV->nrows,A_VH->ncols,pool);

    Q_ASSERT(C_HH->nrows==A_HV->nrows);
    Q_ASSERT(C_HH->ncols==A_VH->ncols);
    mzd_bool_mul_m4rm(C_HH, temp, A_VH, k,blocksize,threads);
    reclaim(temp,nullptr);

    if (! mzd_is_zero(C_HH)) {
        this->mtx[HORIZONTAL][HORIZONTAL] = C_HH;
        Q_ASSERT(C_HH->nrows == mask[HORIZONTAL].len());
        Q_ASSERT(C_HH->ncols == mask[HORIZONTAL].len());
    }
    else {
        reclaim(C_HH,pool);
    }
}


void Graph::combine (const Graph* B)
{
    mzd_t* A_VV = this->mtx[VERTICAL][VERTICAL];
    mzd_t* B_VV = B->mtx[VERTICAL][VERTICAL];
    if (!A_VV) return;
    if (!B_VV) {
        release(VERTICAL,VERTICAL);
        return;
    }

    //  VV are upper triangular matrices
    int threads = ConcurrencyContext::countThreads();
    mzd_bool_and_uptri(A_VV, A_VV, B_VV, threads);

    if (mzd_is_zero(A_VV))
        release(VERTICAL,VERTICAL);
}

void Graph::queryDiagonalElement() const
{
    diagonalElement = -1;
    mzd_t* M = mtx[HORIZONTAL][HORIZONTAL];
    if (!M) return;

    for (int k = 0; k < std::min(M->ncols,M->nrows); ++k)
        if (mzd_read_bit(M, k, k)) {
            diagonalElement = k + mask[HORIZONTAL].lower;
            return;
        }
}

int Graph::foundDiagonalElement() const {
    return diagonalElement;
}

void Graph::add_edge(Orientation ori1, int i1,
                                     Orientation ori2, int i2)
{
    IndexRange r1 (ori1,i1,i1+1);
    IndexRange r2 (ori2,i2,i2+1);
    add_range(r1,r2);
}

bool Graph::contains_edge(Orientation ori1, int i1, Orientation ori2, int i2) const
{
    mzd_t* M = mtx[ori1][ori2];
    if (!M) return false;

    if (!mask[ori1].contains(i1) || !mask[ori2].contains(i2))
        return false;

    int row = i1-mask[ori1].lower;
    int col = i2-mask[ori2].lower;
    return mzd_read_bit(M,row,col);
}


bool Graph::contains_edge(IndexRange qfrom, IndexRange qto) const
{   
    mzd_t* M = mtx[qfrom.ori][qto.ori];
    if (!M) return false;

    qfrom &= mask[qfrom.ori];
    qto &= mask[qto.ori];

    qfrom -= mask[qfrom.ori].lower;
    qto -= mask[qto.ori].lower;

    Q_ASSERT(qfrom.lower >= 0 && qfrom.upper <= mask[qfrom.ori].len());
    Q_ASSERT(qto.lower >= 0 && qto.upper <= mask[qto.ori].len());

    for(int row=qfrom.lower; row < qfrom.upper; ++row)
        if (has_bits(M,row,qto.lower,qto.upper))
            return true;
    return false;
}

bool Graph::empty(Orientation o1, Orientation o2) const
{
    return ! mtx[o1][o2];
}

bool Graph::zero(Orientation o1, Orientation o2) const
{
    return mtx[o1][o2] && mzd_is_zero(mtx[o1][o2]);
}

void Graph::add_range(IndexRange ri, IndexRange rj)
{
    if (ri.empty() || rj.empty())
        return;

    //  Clip to horizontal upper bound.
    //  Reason: free-space boundaries are represented by a LOWER_UPPER interval.
    //  but we do not want RGs to overlap. That's why this (and only this) interval is clipped!

    if (ri.ori==HORIZONTAL && ri.upper > mask[HORIZONTAL].upper) {
        Q_ASSERT(ri.upper == mask[HORIZONTAL].upper+1);
        // assert bias(ri) == LOWER_UPPER
        ri.upper = mask[HORIZONTAL].upper;
    }
    if (rj.ori==HORIZONTAL && rj.upper > mask[HORIZONTAL].upper) {
        Q_ASSERT(rj.upper == mask[HORIZONTAL].upper+1);
        // assert bias(ri) == LOWER_UPPER
        rj.upper = mask[HORIZONTAL].upper;
    }

    Q_ASSERT(mask[ri.ori].contains(ri));
    Q_ASSERT(mask[rj.ori].contains(rj));

    //  map to mask
    ri -= mask[ri.ori].lower;
    rj -= mask[rj.ori].lower;

    mzd_t* M = mtx[ri.ori][rj.ori];
    Q_ASSERT(M);
    if (ri.ori==rj.ori)
    {
        //  HH is ignored
        Q_ASSERT(ri.ori==VERTICAL && rj.ori==VERTICAL);
        //  VV is an upper triangular matrix !!
        for(int i=ri.lower; i < ri.upper; ++i)
        {
            int j = std::max(i,rj.lower);
            if (j < rj.upper)
                set_bits(M, i, j, rj.upper);
        }
    }
    else
    {
        //  HV and VH are regular matrices
        for(int i=ri.lower; i < ri.upper; ++i)
            set_bits(M, i,rj.lower, rj.upper);
    }

}

void Graph::set_bits(mzd_t* M, int r, int i, int j)
{
    Q_ASSERT(i >= 0 && i <= M->ncols);
    Q_ASSERT(j >= i && j <= M->ncols);
    Q_ASSERT(r >= 0 && r < M->nrows);

    wi_t w1 = i / m4ri_radix;
    wi_t w2 = j / m4ri_radix;

    word mask1 = ((word)-1) << (i % m4ri_radix);
    word mask2 = (((word)1) << (j % m4ri_radix)) - 1;
    word* row = M->rows[r];
    Q_ASSERT(mzd_row(M,r)==row);
    Q_ASSERT(w1>=0 && w1<M->width);
    Q_ASSERT(w2>=0 && w2<M->width);

    if (w1==w2) {
        row[w1] |= mask1 & mask2;
    }
    else {
        row[w1] |= mask1;
        row[w2] |= mask2;

        //memset(row+w1+1, (word)-1, w2-w1-1);
        for(wi_t w=w1+1; w < w2; ++w)
            row[w] = ((word)-1);
    }
}

bool Graph::has_bits(mzd_t* M, int r, int i, int j)
{
    Q_ASSERT(i >= 0 && i <= M->ncols);
    Q_ASSERT(j >= i && j <= M->ncols);
    Q_ASSERT(r >= 0 && r < M->nrows);

    wi_t w1 = i / m4ri_radix;
    wi_t w2 = j / m4ri_radix;

    word mask1 = ((word)-1) << (i % m4ri_radix);
    word mask2 = (((word)1) << (j % m4ri_radix)) - 1;
    word* row = M->rows[r];
    Q_ASSERT(mzd_row(M,r)==row);

    if (w1==w2)
        return (row[w1] & mask1 & mask2);

    if (row[w1] & mask1) return true;
    if (row[w2] & mask2) return true;

    for(wi_t w=w1+1; w < w2; ++w)
        if (row[w]) return true;

    return false;
}


double Graph::memory(Orientation ori1, Orientation ori2) const
{
    mzd_t* m = mtx[ori1][ori2];
    Q_ASSERT(m);
    return m->nrows * m->rowstride * sizeof(word);
}

double Graph::density(Orientation ori1, Orientation ori2) const
{
    mzd_t* m = mtx[ori1][ori2];
    Q_ASSERT(m);
    return mzd_density(m,0);
    //  If res = 0 then 100 samples per row are made, if res > 0 the function takes res sized steps within each row (res = 1 uses every word).
}

double Graph::memory() const
{
    double result=0.0;
    for(Orientation o1=HORIZONTAL; o1 <= VERTICAL; ++o1)
        for(Orientation o2=HORIZONTAL; o2 <= VERTICAL; ++o2)
            result += memory(o1,o2);
    return result;
}

double Graph::density() const
{
    double result=0.0;
    for(Orientation o1=HORIZONTAL; o1 <= VERTICAL; ++o1)
        for(Orientation o2=HORIZONTAL; o2 <= VERTICAL; ++o2)
        {
            mzd_t* m = mtx[o1][o2];
            result += density(o1,o2) * m->nrows*m->ncols;
        }
    int nrows = model->count2(HORIZONTAL) + model->count2(VERTICAL);
    return result / (nrows*nrows);
}

void Graph::clear()
{
    for(Orientation o1=HORIZONTAL; o1 <= VERTICAL; ++o1)
        for(Orientation o2=HORIZONTAL; o2 <= VERTICAL; ++o2)
            clear(o1,o2);
}

void Graph::clear(Orientation o1, Orientation o2)
{
    if (mtx[o1][o2])
        mzd_set_ui(mtx[o1][o2],0);
}

bool Graph::is_upper_triangular(Orientation ori1, Orientation ori2) const
{
    return is_upper_triangular(mtx[ori1][ori2]);
}

bool Graph::is_upper_triangular(const mzd_t* m)
{
    if (m->nrows != m->ncols)
        return false;

    for(int i=1; i < m->nrows; ++i)
        for (int j=0; j < i; ++j)
            if (mzd_read_bit(m,i,j))
                return false;
    return true;
}

int Graph::find_next1(Orientation o1, int i, Orientation o2, int j) const
{
    mzd_t* M = mtx[o1][o2];
    if (!M) return INT_MAX;

    if ((i >= mask[o1].upper) || (j >= mask[o2].upper))
        return INT_MAX;

    i = std::max(i, mask[o1].lower);
    j = std::max(j, mask[o2].lower);

    int p = find_next1(M, i-mask[o1].lower, j-mask[o2].lower);
    if (p!=INT_MAX) p += mask[o2].lower;
    return p;
}

int Graph::find_next0(Orientation o1, int i, Orientation o2, int j) const
{
    mzd_t* M = mtx[o1][o2];
    if (!M) return j;

    if (! mask[o1].contains(i) || ! mask[o2].contains(j))
        return j;

    int p = find_next0(M, i-mask[o1].lower, j-mask[o2].lower);
    if (p==INT_MAX)
        return mask[o2].upper;
    else
        return p + mask[o2].lower;
}

IndexRange Graph::find_next_range(Orientation o1, int i, Orientation o2, int j) const
{
    int j1 = find_next1(o1,i,o2,j);
    if (j1==INT_MAX)
        return IndexRange(o2,-1,-1);
    int j2 = find_next0(o1,i,o2,j1+1);
    Q_ASSERT(j2 <= mask[o2].upper);
    return IndexRange(o2,j1,j2);
}


#define INVERT(w) ((w) ^ ((uint64_t)-1))

inline int lsb_pos(uint64_t x)
{
    Q_ASSERT(x != 0);
    return boost::multiprecision::lsb(x);
}

int Graph::find_next1(mzd_t* M, int r, int i)
{
    Q_ASSERT(r>=0 && r<M->nrows);
    Q_ASSERT(i>=0);

    word* row = M->rows[r];
    wi_t w = i / m4ri_radix;
    if (w >= M->width)
        return INT_MAX;

    word v = row[w] >> (i % m4ri_radix);
    if (v != 0)
        return i + lsb_pos(v);

    while(++w < M->width)
        if (row[w]!=0)
            return w*m4ri_radix + lsb_pos(row[w]);

   return INT_MAX;
}

int Graph::find_next0(mzd_t* M, int r, int i)
{
    Q_ASSERT(r>=0 && r<M->nrows);
    Q_ASSERT(i>=0);

    word* row = M->rows[r];
    wi_t w = i / m4ri_radix;
    if (w >= M->width)
        return INT_MAX;

    word v = (INVERT(row[w]) >> (i % m4ri_radix));
    if (v != 0)
        return i + lsb_pos(v) ;

    while(++w < M->width)
        if (row[w]!=(uint64_t)-1)
            return w*m4ri_radix + lsb_pos(INVERT(row[w]));

    return INT_MAX;
}

bool Graph::is_adjacent_to(const Graph& A, const Graph& B)
{
    return (A.hmask().upper==B.hmask().lower)
        || (A.hmask().upper==B.graphModel()->lowerShiftedHorizontally(B.hmask().lower));
}

void Graph::print(std::ostream& out, Orientation o1, Orientation o2) const {
    const mzd_t* M = mtx[o1][o2];
    print(out,M);
}

void Graph::print(std::ostream& out, const mzd_t* M)
{
    for(int row=0; row < M->nrows; ++row) {
        for (int col = 0; col < M->ncols; col += 64) {
            word w = mzd_read_bits(M, row, col, 64);
            out << std::hex << std::setw(16) << w << " ";
        }
        out << std::endl;
    }
    out << std::endl;
}

void Graph::release(Orientation o1, Orientation o2, MatrixPool* pool) {
    if (mtx[o1][o2]) {
        reclaim(mtx[o1][o2],pool);
        mtx[o1][o2] = nullptr;
    }
}