Fréchet View  1.6.0
A Tool for Exploring Fréchet Distance Algorithms
graph_m4ri.cpp
Go to the documentation of this file.
1 
3 #include <concurrency.h>
4 // m4ri extension to Boolean Algebra
5 #if defined(__cplusplus) && ! defined(_MSC_VER)
6 extern "C" {
7 #endif
8 # include <m4ri/m4ri.h>
9 # include <m4ri/mzd_bool.h>
10 # include <m4ri/russian_bool.h>
11 #if defined(__cplusplus) && !defined(_MSC_VER)
12 }
13 #endif
14 
15 #include <iomanip>
16 #include <boost/cstdint.hpp>
17 #include <boost/multiprecision/integer.hpp>
18 
19 using namespace frechet;
20 using namespace reach;
21 using namespace app;
22 
24  : model(amodel),
25  origin{Origin::RG,-1,0,-1,-1,nullptr,nullptr,nullptr},
26  diagonalElement(-1)
27 {
28  for(Orientation o1=HORIZONTAL; o1 <= VERTICAL; ++o1)
29  mask[o1] = IndexRange(o1, 0, model->count2(o1));
30  for(Orientation o1=HORIZONTAL; o1 <= VERTICAL; ++o1)
31  for(Orientation o2=HORIZONTAL; o2 <= VERTICAL; ++o2)
32  mtx[o1][o2] = nullptr;
33 }
34 
36  IndexRange range)
37  : Graph(amodel)
38 {
39  mask[HORIZONTAL] = range;
40 
41  //Q_ASSERT(mask[HORIZONTAL].lower < model->count1(HORIZONTAL));
42  Q_ASSERT(mask[HORIZONTAL].len() <= model->count2(HORIZONTAL));
43 }
44 
46  Structure& str)
47  : Graph(amodel)
48 {
50  mask[HORIZONTAL].lower = model->map_lower(HORIZONTAL,str.bottom().first()->lower());
51  mask[HORIZONTAL].upper = model->map_lower(HORIZONTAL,str.bottom().last()->upper());
52  // Note: upper bound is exclusive.
53  // Don't overlap with next RG.
54  //Q_ASSERT(mask[HORIZONTAL].lower < model->count1(HORIZONTAL));
55  Q_ASSERT(mask[HORIZONTAL].len() <= model->count2(HORIZONTAL));
56 
60 
61  // insert edges into the graph (EXCEPT H-H egdes!)
62  copy(str,VERTICAL);
65 }
66 
67 Graph::Graph(const Graph& that)
68  : Graph(that.model)
69 {
70  copy(that);
71 }
72 
74  : Graph(that.model)
75 {
76  swap(that);
77 }
78 
80  for(Orientation o1=HORIZONTAL; o1 <= VERTICAL; ++o1)
81  for(Orientation o2=HORIZONTAL; o2 <= VERTICAL; ++o2)
82  allocate(o1,o2);
83 }
84 
86  Q_ASSERT(mtx[o1][o2]==nullptr);
87  mtx[o1][o2] = mzd_init(mask[o1].len(),mask[o2].len());
88 
89  Q_ASSERT(mtx[o1][o2]->nrows==mask[o1].len());
90  Q_ASSERT(mtx[o1][o2]->ncols==mask[o2].len());
91  return mtx[o1][o2];
92 }
93 
95  return Rect(mask[o1].lower,mask[o2].lower,
96  mask[o1].upper,mask[o2].upper);
97 }
98 
100 {
101  release();
102 }
103 
104 void Graph::setOriginPlacement(int di, int dj)
105 {
107  origin.blevel = -1;
108  origin.succesors = 0;
109  origin.i = di;
110  origin.j = dj;
111  origin.A.reset();
112  origin.B.reset();
113  origin.P.reset();
114 }
115 
117 {
119  origin.blevel = 0;
120  origin.succesors = 0;
121  origin.i = i;
122  origin.j = i+1;
123  origin.A.reset();
124  origin.B.reset();
125  origin.P.reset();
126 }
127 
129 {
131  origin.blevel = std::max(A->origin.blevel, B->origin.blevel) + 1;
132  origin.succesors = 0;
133  origin.i = A->origin.i;
134  origin.j = A->origin.j + (B->origin.j - B->origin.i);
135  origin.A = A;
136  origin.B = B;
137  origin.P.reset();
138  A->origin.succesors++; // count dependencies on A_VV
139  B->origin.succesors++; // count dependencies on B_VV
140 }
141 
143 {
145  Q_ASSERT(A->origin.blevel == 0);
146  origin.blevel = B->origin.blevel + 1;
147  origin.succesors = 0;
148  origin.i = A->origin.i;
149  origin.j = A->origin.i;
150  origin.A = A;
151  origin.B = B;
152  origin.P.reset();
153  //A->origin.succesors++; don't count dependecies on A_HV and A_VH
154  B->origin.succesors++; // count dependencies on B_VV
155 }
156 
158 {
159  origin.P = P;
160 }
161 
162 void Graph::printOrigin(std::ostream& out) const
163 {
164  if (origin.P)
165  out << "COMBINE(";
166  switch(origin.operation)
167  {
168  case Origin::RG:
169  out << "RG("<<origin.i<<","<<origin.j<<")";
170  break;
171  case Origin::MERGE2:
172  Q_ASSERT(origin.A.get() != this);
173  Q_ASSERT(origin.B.get() != this);
174  // otherwise we have infinit recursion right here
175  out << "MERGE(";
176  origin.A->printOrigin(out);
177  out << ",";
178  origin.B->printOrigin(out);
179  out << ")";
180  break;
181  case Origin::MERGE3:
182  Q_ASSERT(origin.A.get() != this);
183  Q_ASSERT(origin.B.get() != this);
184  // otherwise we have infinit recursion right here
185  out << "MERGE(";
186  origin.A->printOrigin(out);
187  out << ",";
188  origin.B->printOrigin(out);
189  out << ",";
190  origin.A->printOrigin(out);
191  out << ")";
192  break; default:
193  Q_ASSERT(false);
194  }
195  if (origin.P)
196  out << ")"; // COMBINE
197 }
198 
200  releaseIfZero();
201 }
202 
204 {
205  for(Orientation o1=HORIZONTAL; o1 <= VERTICAL; ++o1)
206  for(Orientation o2=HORIZONTAL; o2 <= VERTICAL; ++o2)
207  release(o1,o2);
208 }
209 
211 {
212  for(Orientation o1=HORIZONTAL; o1 <= VERTICAL; ++o1)
213  for(Orientation o2=HORIZONTAL; o2 <= VERTICAL; ++o2)
214  if (zero(o1,o2))
215  release(o1,o2);
216 }
217 
219 {
220  if (mtx[o1][o2]) {
221  mzd_free(mtx[o1][o2]);
222  mtx[o1][o2] = nullptr;
223  }
224 }
225 
226 
228 {
229  Q_ASSERT(model == that.model);
230  copy(that);
231  return *this;
232 }
233 
235 {
236  swap(that);
237  return *this;
238 }
239 
240 bool Graph::operator== (const Graph& that) const
241 {
242  for(Orientation o1=HORIZONTAL; o1 <= VERTICAL; ++o1)
243  {
244  if (mask[o1] != that.mask[o1])
245  return false;
246 
247  for(Orientation o2=HORIZONTAL; o2 <= VERTICAL; ++o2)
248  {
249  mzd_t* m1 = mtx[o1][o2];
250  mzd_t* m2 = that.mtx[o1][o2];
251  if (!mzd_equal(m1,m2))
252  return false;
253  }
254  }
255  return true;
256 }
257 
258 void Graph::copy(const Graph& that)
259 {
260  for(Orientation o1=HORIZONTAL; o1 <= VERTICAL; ++o1)
261  {
262  this->mask[o1] = that.mask[o1];
263  for (Orientation o2 = HORIZONTAL; o2 <= VERTICAL; ++o2)
264  {
265  mzd_t* m2 = that.mtx[o1][o2];
266  mzd_t*& m1 = mtx[o1][o2];
267 
268  if (m1 && m2
269  && (m1->nrows == m2->nrows)
270  && (m1->ncols == m2->ncols))
271  {
272  mzd_copy(m1, m2);
273  }
274  else
275  {
276  if (m1)
277  mzd_free(m1);
278  if (m2)
279  m1 = mzd_copy(nullptr,m2);
280  else
281  m1 = nullptr;
282  }
283  }
284  }
285 
286  this->origin = that.origin;
287  this->diagonalElement = that.diagonalElement;
288 }
289 
290 void Graph::swap(Graph& that)
291 {
292  Q_ASSERT(&model == &that.model);
293  std::swap(this->mtx, that.mtx);
294  std::swap(this->origin, that.origin);
295  std::swap(this->mask, that.mask);
297 }
298 
300 {
301  /* Copy from Reachability Structure
302  */
303  Pointer i1 = str.first(ori).first();
304  Pointer i2 = str.second(ori).first();
305  for( ; i1; i1=i1->next(), i2=i2->next())
306  {
307  if (i1->type==NON_ACCESSIBLE) continue;
308  // traverse all reachable intervals
309  IndexRange r1 = model->map(i1);
310 
311  // map [l..h] to node indexes
312  StructureIterator lh(str,*i1,i2);
313  for( ; lh; ++lh) {
314  if (i1->ori==HORIZONTAL && lh->ori==HORIZONTAL)
315  continue; // ignore H-H
316  IndexRange r2 = model->map((Pointer)lh);
317  Q_ASSERT(lh->ori==r2.ori);
318  if (lh->type != NON_ACCESSIBLE)
319  add_range(r1,r2);
320  }
321  }
322 }
323 
324 
325 Graph::ptr Graph::merge2(const Graph* B, MatrixPool* pool) const
326 {
327  Graph::ptr C = newGraph(model, model->mergedHMask(this->hmask(),B->hmask()));
328  // new Graph without data
329  C->merge2(this,B,pool);
330  return C;
331 }
332 
333 void Graph::merge2(const Graph* A, const Graph* B, MatrixPool* pool)
334 {
335  Q_ASSERT(this->mtx[VERTICAL][VERTICAL]==nullptr);
336  Q_ASSERT(is_adjacent_to(*A, *B));
337 
338  mzd_t* A_VV = A->mtx[VERTICAL][VERTICAL];
339  mzd_t* B_VV = B->mtx[VERTICAL][VERTICAL];
340  if (! A_VV || ! B_VV)
341  return;
342 
343  mzd_t* C_VV = new_mzd(A_VV->nrows,B_VV->ncols,pool);
344 
345  // both matrices are upper triangular
346  int k, blocksize;
347  int threads = ConcurrencyContext::countThreads();
348  mzd_optimal_parameters(A_VV, B_VV,
349  ConcurrencyContext::cacheSize(2), // right ?
350  ConcurrencyContext::cacheSize(3), // right ?
351  &k, &blocksize);
352 
353  mzd_bool_mul_uptri(C_VV, A_VV, B_VV, k,blocksize,threads);
354 
355  if (! mzd_is_zero(C_VV)) {
356  this->mtx[VERTICAL][VERTICAL] = C_VV;
357  Q_ASSERT(C_VV->nrows == this->mask[VERTICAL].len());
358  Q_ASSERT(C_VV->ncols == this->mask[VERTICAL].len());
359  }
360  else {
361  reclaim(C_VV,pool);
362  }
363 }
364 
365 
366 Graph::ptr Graph::merge3(const Graph* B, MatrixPool* pool) const
367 {
369  C->merge3(this,B,pool);
370  return C;
371 }
372 
373 void Graph::merge3(const Graph* A, const Graph* B, MatrixPool* pool)
374 {
375  Q_ASSERT(this->mtx[HORIZONTAL][HORIZONTAL]==nullptr);
376  Q_ASSERT(is_adjacent_to(*A, *B));
377  Q_ASSERT(is_adjacent_to(*B, *A));
378 
379  mzd_t* A_HV = A->mtx[HORIZONTAL][VERTICAL];
380  mzd_t* B_VV = B->mtx[VERTICAL][VERTICAL];
381  mzd_t* A_VH = A->mtx[VERTICAL][HORIZONTAL];
382 
383  if (! A_HV || ! B_VV || ! A_VH)
384  return; // shortcut for empty
385 
386  int k, blocksize;
387  int threads = ConcurrencyContext::countThreads();
388  mzd_optimal_parameters(A_HV, B_VV,
389  ConcurrencyContext::cacheSize(2), // right ?
390  ConcurrencyContext::cacheSize(3), // right ?
391  &k, &blocksize);
392 
393  mzd_t* temp = new_mzd(A_HV->nrows,B_VV->ncols,nullptr); // don't pool odd-sized matrices
394  mzd_bool_mul_m4rm(temp, A_HV, B_VV, k,blocksize,threads);
395 
396  mzd_optimal_parameters(temp, A_VH,
397  ConcurrencyContext::cacheSize(2), // right ?
398  ConcurrencyContext::cacheSize(3), // right ?
399  &k, &blocksize);
400 
401  mzd_t* C_HH = new_mzd(A_HV->nrows,A_VH->ncols,pool);
402 
403  Q_ASSERT(C_HH->nrows==A_HV->nrows);
404  Q_ASSERT(C_HH->ncols==A_VH->ncols);
405  mzd_bool_mul_m4rm(C_HH, temp, A_VH, k,blocksize,threads);
406  reclaim(temp,nullptr);
407 
408  if (! mzd_is_zero(C_HH)) {
409  this->mtx[HORIZONTAL][HORIZONTAL] = C_HH;
410  Q_ASSERT(C_HH->nrows == mask[HORIZONTAL].len());
411  Q_ASSERT(C_HH->ncols == mask[HORIZONTAL].len());
412  }
413  else {
414  reclaim(C_HH,pool);
415  }
416 }
417 
418 
419 void Graph::combine (const Graph* B)
420 {
421  mzd_t* A_VV = this->mtx[VERTICAL][VERTICAL];
422  mzd_t* B_VV = B->mtx[VERTICAL][VERTICAL];
423  if (!A_VV) return;
424  if (!B_VV) {
426  return;
427  }
428 
429  // VV are upper triangular matrices
430  int threads = ConcurrencyContext::countThreads();
431  mzd_bool_and_uptri(A_VV, A_VV, B_VV, threads);
432 
433  if (mzd_is_zero(A_VV))
435 }
436 
438 {
439  diagonalElement = -1;
440  mzd_t* M = mtx[HORIZONTAL][HORIZONTAL];
441  if (!M) return;
442 
443  for (int k = 0; k < std::min(M->ncols,M->nrows); ++k)
444  if (mzd_read_bit(M, k, k)) {
446  return;
447  }
448 }
449 
451  return diagonalElement;
452 }
453 
454 void Graph::add_edge(Orientation ori1, int i1,
455  Orientation ori2, int i2)
456 {
457  IndexRange r1 (ori1,i1,i1+1);
458  IndexRange r2 (ori2,i2,i2+1);
459  add_range(r1,r2);
460 }
461 
462 bool Graph::contains_edge(Orientation ori1, int i1, Orientation ori2, int i2) const
463 {
464  mzd_t* M = mtx[ori1][ori2];
465  if (!M) return false;
466 
467  if (!mask[ori1].contains(i1) || !mask[ori2].contains(i2))
468  return false;
469 
470  int row = i1-mask[ori1].lower;
471  int col = i2-mask[ori2].lower;
472  return mzd_read_bit(M,row,col);
473 }
474 
475 
476 bool Graph::contains_edge(IndexRange qfrom, IndexRange qto) const
477 {
478  mzd_t* M = mtx[qfrom.ori][qto.ori];
479  if (!M) return false;
480 
481  qfrom &= mask[qfrom.ori];
482  qto &= mask[qto.ori];
483 
484  qfrom -= mask[qfrom.ori].lower;
485  qto -= mask[qto.ori].lower;
486 
487  Q_ASSERT(qfrom.lower >= 0 && qfrom.upper <= mask[qfrom.ori].len());
488  Q_ASSERT(qto.lower >= 0 && qto.upper <= mask[qto.ori].len());
489 
490  for(int row=qfrom.lower; row < qfrom.upper; ++row)
491  if (has_bits(M,row,qto.lower,qto.upper))
492  return true;
493  return false;
494 }
495 
497 {
498  return ! mtx[o1][o2];
499 }
500 
502 {
503  return mtx[o1][o2] && mzd_is_zero(mtx[o1][o2]);
504 }
505 
507 {
508  if (ri.empty() || rj.empty())
509  return;
510 
511  // Clip to horizontal upper bound.
512  // Reason: free-space boundaries are represented by a LOWER_UPPER interval.
513  // but we do not want RGs to overlap. That's why this (and only this) interval is clipped!
514 
515  if (ri.ori==HORIZONTAL && ri.upper > mask[HORIZONTAL].upper) {
516  Q_ASSERT(ri.upper == mask[HORIZONTAL].upper+1);
517  // assert bias(ri) == LOWER_UPPER
518  ri.upper = mask[HORIZONTAL].upper;
519  }
520  if (rj.ori==HORIZONTAL && rj.upper > mask[HORIZONTAL].upper) {
521  Q_ASSERT(rj.upper == mask[HORIZONTAL].upper+1);
522  // assert bias(ri) == LOWER_UPPER
523  rj.upper = mask[HORIZONTAL].upper;
524  }
525 
526  Q_ASSERT(mask[ri.ori].contains(ri));
527  Q_ASSERT(mask[rj.ori].contains(rj));
528 
529  // map to mask
530  ri -= mask[ri.ori].lower;
531  rj -= mask[rj.ori].lower;
532 
533  mzd_t* M = mtx[ri.ori][rj.ori];
534  Q_ASSERT(M);
535  if (ri.ori==rj.ori)
536  {
537  // HH is ignored
538  Q_ASSERT(ri.ori==VERTICAL && rj.ori==VERTICAL);
539  // VV is an upper triangular matrix !!
540  for(int i=ri.lower; i < ri.upper; ++i)
541  {
542  int j = std::max(i,rj.lower);
543  if (j < rj.upper)
544  set_bits(M, i, j, rj.upper);
545  }
546  }
547  else
548  {
549  // HV and VH are regular matrices
550  for(int i=ri.lower; i < ri.upper; ++i)
551  set_bits(M, i,rj.lower, rj.upper);
552  }
553 
554 }
555 
556 void Graph::set_bits(mzd_t* M, int r, int i, int j)
557 {
558  Q_ASSERT(i >= 0 && i <= M->ncols);
559  Q_ASSERT(j >= i && j <= M->ncols);
560  Q_ASSERT(r >= 0 && r < M->nrows);
561 
562  wi_t w1 = i / m4ri_radix;
563  wi_t w2 = j / m4ri_radix;
564 
565  word mask1 = ((word)-1) << (i % m4ri_radix);
566  word mask2 = (((word)1) << (j % m4ri_radix)) - 1;
567  word* row = M->rows[r];
568  Q_ASSERT(mzd_row(M,r)==row);
569  Q_ASSERT(w1>=0 && w1<M->width);
570  Q_ASSERT(w2>=0 && w2<M->width);
571 
572  if (w1==w2) {
573  row[w1] |= mask1 & mask2;
574  }
575  else {
576  row[w1] |= mask1;
577  row[w2] |= mask2;
578 
579  //memset(row+w1+1, (word)-1, w2-w1-1);
580  for(wi_t w=w1+1; w < w2; ++w)
581  row[w] = ((word)-1);
582  }
583 }
584 
585 bool Graph::has_bits(mzd_t* M, int r, int i, int j)
586 {
587  Q_ASSERT(i >= 0 && i <= M->ncols);
588  Q_ASSERT(j >= i && j <= M->ncols);
589  Q_ASSERT(r >= 0 && r < M->nrows);
590 
591  wi_t w1 = i / m4ri_radix;
592  wi_t w2 = j / m4ri_radix;
593 
594  word mask1 = ((word)-1) << (i % m4ri_radix);
595  word mask2 = (((word)1) << (j % m4ri_radix)) - 1;
596  word* row = M->rows[r];
597  Q_ASSERT(mzd_row(M,r)==row);
598 
599  if (w1==w2)
600  return (row[w1] & mask1 & mask2);
601 
602  if (row[w1] & mask1) return true;
603  if (row[w2] & mask2) return true;
604 
605  for(wi_t w=w1+1; w < w2; ++w)
606  if (row[w]) return true;
607 
608  return false;
609 }
610 
611 
612 double Graph::memory(Orientation ori1, Orientation ori2) const
613 {
614  mzd_t* m = mtx[ori1][ori2];
615  Q_ASSERT(m);
616  return m->nrows * m->rowstride * sizeof(word);
617 }
618 
619 double Graph::density(Orientation ori1, Orientation ori2) const
620 {
621  mzd_t* m = mtx[ori1][ori2];
622  Q_ASSERT(m);
623  return mzd_density(m,0);
624  // 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).
625 }
626 
627 double Graph::memory() const
628 {
629  double result=0.0;
630  for(Orientation o1=HORIZONTAL; o1 <= VERTICAL; ++o1)
631  for(Orientation o2=HORIZONTAL; o2 <= VERTICAL; ++o2)
632  result += memory(o1,o2);
633  return result;
634 }
635 
636 double Graph::density() const
637 {
638  double result=0.0;
639  for(Orientation o1=HORIZONTAL; o1 <= VERTICAL; ++o1)
640  for(Orientation o2=HORIZONTAL; o2 <= VERTICAL; ++o2)
641  {
642  mzd_t* m = mtx[o1][o2];
643  result += density(o1,o2) * m->nrows*m->ncols;
644  }
645  int nrows = model->count2(HORIZONTAL) + model->count2(VERTICAL);
646  return result / (nrows*nrows);
647 }
648 
650 {
651  for(Orientation o1=HORIZONTAL; o1 <= VERTICAL; ++o1)
652  for(Orientation o2=HORIZONTAL; o2 <= VERTICAL; ++o2)
653  clear(o1,o2);
654 }
655 
657 {
658  if (mtx[o1][o2])
659  mzd_set_ui(mtx[o1][o2],0);
660 }
661 
663 {
664  return is_upper_triangular(mtx[ori1][ori2]);
665 }
666 
667 bool Graph::is_upper_triangular(const mzd_t* m)
668 {
669  if (m->nrows != m->ncols)
670  return false;
671 
672  for(int i=1; i < m->nrows; ++i)
673  for (int j=0; j < i; ++j)
674  if (mzd_read_bit(m,i,j))
675  return false;
676  return true;
677 }
678 
679 int Graph::find_next1(Orientation o1, int i, Orientation o2, int j) const
680 {
681  mzd_t* M = mtx[o1][o2];
682  if (!M) return INT_MAX;
683 
684  if ((i >= mask[o1].upper) || (j >= mask[o2].upper))
685  return INT_MAX;
686 
687  i = std::max(i, mask[o1].lower);
688  j = std::max(j, mask[o2].lower);
689 
690  int p = find_next1(M, i-mask[o1].lower, j-mask[o2].lower);
691  if (p!=INT_MAX) p += mask[o2].lower;
692  return p;
693 }
694 
695 int Graph::find_next0(Orientation o1, int i, Orientation o2, int j) const
696 {
697  mzd_t* M = mtx[o1][o2];
698  if (!M) return j;
699 
700  if (! mask[o1].contains(i) || ! mask[o2].contains(j))
701  return j;
702 
703  int p = find_next0(M, i-mask[o1].lower, j-mask[o2].lower);
704  if (p==INT_MAX)
705  return mask[o2].upper;
706  else
707  return p + mask[o2].lower;
708 }
709 
711 {
712  int j1 = find_next1(o1,i,o2,j);
713  if (j1==INT_MAX)
714  return IndexRange(o2,-1,-1);
715  int j2 = find_next0(o1,i,o2,j1+1);
716  Q_ASSERT(j2 <= mask[o2].upper);
717  return IndexRange(o2,j1,j2);
718 }
719 
720 
721 #define INVERT(w) ((w) ^ ((uint64_t)-1))
722 
727 inline int lsb_pos(uint64_t x)
728 {
729  Q_ASSERT(x != 0);
730  return boost::multiprecision::lsb(x);
731 }
732 
733 int Graph::find_next1(mzd_t* M, int r, int i)
734 {
735  Q_ASSERT(r>=0 && r<M->nrows);
736  Q_ASSERT(i>=0);
737 
738  word* row = M->rows[r];
739  wi_t w = i / m4ri_radix;
740  if (w >= M->width)
741  return INT_MAX;
742 
743  word v = row[w] >> (i % m4ri_radix);
744  if (v != 0)
745  return i + lsb_pos(v);
746 
747  while(++w < M->width)
748  if (row[w]!=0)
749  return w*m4ri_radix + lsb_pos(row[w]);
750 
751  return INT_MAX;
752 }
753 
754 int Graph::find_next0(mzd_t* M, int r, int i)
755 {
756  Q_ASSERT(r>=0 && r<M->nrows);
757  Q_ASSERT(i>=0);
758 
759  word* row = M->rows[r];
760  wi_t w = i / m4ri_radix;
761  if (w >= M->width)
762  return INT_MAX;
763 
764  word v = (INVERT(row[w]) >> (i % m4ri_radix));
765  if (v != 0)
766  return i + lsb_pos(v) ;
767 
768  while(++w < M->width)
769  if (row[w]!=(uint64_t)-1)
770  return w*m4ri_radix + lsb_pos(INVERT(row[w]));
771 
772  return INT_MAX;
773 }
774 
775 bool Graph::is_adjacent_to(const Graph& A, const Graph& B)
776 {
777  return (A.hmask().upper==B.hmask().lower)
778  || (A.hmask().upper==B.graphModel()->lowerShiftedHorizontally(B.hmask().lower));
779 }
780 
781 void Graph::print(std::ostream& out, Orientation o1, Orientation o2) const {
782  const mzd_t* M = mtx[o1][o2];
783  print(out,M);
784 }
785 
786 void Graph::print(std::ostream& out, const mzd_t* M)
787 {
788  for(int row=0; row < M->nrows; ++row) {
789  for (int col = 0; col < M->ncols; col += 64) {
790  word w = mzd_read_bits(M, row, col, 64);
791  out << std::hex << std::setw(16) << w << " ";
792  }
793  out << std::endl;
794  }
795  out << std::endl;
796 }
797 
799  if (mtx[o1][o2]) {
800  reclaim(mtx[o1][o2],pool);
801  mtx[o1][o2] = nullptr;
802  }
803 }
int find_next0(Orientation o1, int i, Orientation o2, int j) const
find the next missing edge in a sub-graph
Definition: graph_m4ri.cpp:695
bool contains(int i) const
containment test
Definition: graph_model.cpp:26
void swap(gpuword **A, gpuword **B)
int count2(Orientation ori) const
void print(std::ostream &out, Orientation o1, Orientation o2) const
print a sub-adjacancy matrix
Definition: graph_m4ri.cpp:781
Graph::ptr P
valid placement applied to the graph
Definition: graph_m4ri.h:110
mzd_t * allocate(Orientation o1, Orientation o2)
allocate an adjacancy sub-matrix (one of four parts)
Definition: graph_m4ri.cpp:85
int i
if this is a reachability graph: column range in free-space. this = RG(i,j), or CRG(i,...
Definition: graph_m4ri.h:102
void copy(const Graph &)
virtual void combine(const Graph *P)
apply the COMBINE operation, filtering edges with valid placements. Effectively performs a Boolean AN...
Definition: graph_m4ri.cpp:419
global definitions for all algorithms.
mzd_t * mtx[2][2]
adjacency matrix (M4RI structure) split into four parts to allow for memory savings.
Definition: graph_m4ri.h:45
const Orientation ori
horizontal or vertical
Definition: boundary.h:318
non-accesible. no interval is reachable from this one
Definition: boundary.h:19
int map_lower(Pointer p) const
map the lower bound of a reachability structure interval
boost::shared_ptr< GraphModel > ptr
smart pointer to a GraphModel object
Definition: graph_model.h:307
int diagonalElement
result of call to searchDiagonalElement
Definition: graph_m4ri.h:78
The StructureIterator class.
Definition: structure.h:488
void releaseIfZero()
release memory that is not neeeded (empty sub-graphs)
Definition: graph_m4ri.cpp:210
bool zero(Orientation o1, Orientation o2) const
empty test for one of the four sub-graphs
Definition: graph_m4ri.cpp:501
void setOriginMerge3(Graph::ptr A, Graph::ptr B)
this graph is the transitive closure of two other graphs. It is the result of a final MERGE operation...
Definition: graph_m4ri.cpp:142
boundary interval in the reachability structure. Represents an interval on the boundary of the FreeSp...
Definition: boundary.h:285
mzd_t * new_mzd(int rows, int cols, MatrixPool *pool)
allocate a new mzd_t structure (a matrix for the M4RI algorithms)
bool empty(Orientation o1, Orientation o2) const
empty test for one of the four sub-graphs
Definition: graph_m4ri.cpp:496
bool contains_edge(Orientation ori_from, int i_from, Orientation ori_to, int i_to) const
check if an edge is present
bool is_upper_triangular(Orientation ori1, Orientation ori2) const
Definition: graph_m4ri.cpp:662
int lsb_pos(uint64_t x)
Definition: graph_m4ri.cpp:727
int upper
upper index (exclusive)
Definition: graph_model.h:23
int lower
lower index (inclusive)
Definition: graph_model.h:21
#define INVERT(w)
Definition: graph_m4ri.cpp:721
GraphPtr newGraph(const GraphModel::ptr model)
Definition: graph_cl.cpp:12
Type type
reachability label
Definition: boundary.h:317
virtual void release()
release memory for all parts of the adjacancy matrix
Definition: graph_m4ri.cpp:203
a very simple Rectangle structure, with integer boundaries.
Definition: types.h:35
Graph::ptr merge2(const Graph *B, MatrixPool *pool=nullptr) const
apply the MERGE operation, computing the transitive closure of two graphs. Effectively performs a mat...
double memory() const
Definition: graph_m4ri.cpp:627
Graph::ptr merge3(const Graph *B, MatrixPool *pool=nullptr) const
apply the final MERGE operation, computing the transitive closure of two graphs. Effectively performs...
void allocateAll()
allocate data for adjacenty matrixes (all four parts)
Definition: graph_m4ri.cpp:79
void setOriginRG(int i)
this graph represents a reachability graph, it was created from a reachability structure
Definition: graph_m4ri.cpp:116
void setOriginPlacement(int di, int dj)
this graph represents a set of valid placements
Definition: graph_m4ri.cpp:104
a range of node indices in a Reachability Graph
Definition: graph_model.h:17
void setOriginCombine(Graph::ptr P)
a sub-set of valid placements is applied to the graph
Definition: graph_m4ri.cpp:157
Orientation
Segment Orientation.
Definition: boundary.h:31
IndexRange map(Orientation ori, const Interval &, bool upperBoundInclusive) const
map a free-space interval to a range of graph nodes
boost::shared_ptr< Graph > ptr
Definition: graph_boost.h:59
this is a rechability graph, constructed from a reachability structure
Definition: graph_m4ri.h:92
memory pool for matrix objects (M4RI matrices mzd_t* and OpenCL matrices clm4rm_t*)
Definition: matrix_pool.h:26
double min(double a, double b)
minimum function with checks for NAN
Definition: numeric.h:222
void reclaim(mzd_t *m, MatrixPool *pool)
reclaim an object (i.e. put it into the recycling list)
Graph(const GraphModel &model, graph_t *ag)
Definition: graph_boost.cpp:10
Represents a Reachability Graph. Vertices correspond to intervals in the reachability structure,...
Definition: graph_boost.h:39
void add_range(IndexRange r_from, IndexRange r_to)
add a set of edges to the graph, covering a range of nodes
Definition: graph_m4ri.cpp:506
this graph is the result of a MERGE operation (transitive closure) on two graphs
Definition: graph_m4ri.h:93
IndexRange mergedHMask(const IndexRange &m1, const IndexRange &m2) const
merged two node ranges representing the horizontal masks of reachability graph (see frechet::reach::G...
#define str(S)
double density() const
Definition: graph_m4ri.cpp:636
Orientation ori
horizontal or vertical part of the reachability structure?
Definition: graph_model.h:19
IndexRange find_next_range(Orientation o1, int i, Orientation o2, int j) const
find a range of consecutive edges in a sub-graph
Definition: graph_m4ri.cpp:710
void setOriginMerge2(Graph::ptr A, Graph::ptr B)
this graph is the transitive closure of two other graphs. It is the result of a MERGE operation.
Definition: graph_m4ri.cpp:128
void add_edge(int i, int j)
this graph represents a set of valid placements
Definition: graph_m4ri.h:91
struct frechet::reach::Graph::Origin origin
int find_next1(Orientation o1, int i, Orientation o2, int j) const
find the next edge in a sub-graph
Definition: graph_m4ri.cpp:679
static bool has_bits(mzd_t *M, int row, int i, int j)
query bits in a Boolean matrix
Definition: graph_m4ri.cpp:585
Graph::ptr A
if this graph was constructed as transitive closure: original graphs this = MERGE(A,...
Definition: graph_m4ri.h:106
int blevel
bottom-level (for topologial sorting; leaves = blevel 0
Definition: graph_m4ri.h:97
virtual void queryDiagonalElement() const
find an edge on the diagonal of the adjacancy matrix. Does not return a result. To query the result o...
Definition: graph_m4ri.cpp:437
const GraphModel & model
Definition: graph_boost.h:55
static void set_bits(mzd_t *M, int row, int i, int j)
set some bits in a Boolean matrix
Definition: graph_m4ri.cpp:556
const GraphModel::ptr graphModel() const
Definition: graph_m4ri.h:198
The Reachability Structure; maintains a list of intervals on the border of Free Space,...
Definition: structure.h:32
const IndexRange & hmask() const
horizontal range covered by this graph
Definition: graph_m4ri.h:187
int succesors
number of successors
Definition: graph_m4ri.h:99
void printOrigin(std::ostream &out) const
debug output
Definition: graph_m4ri.cpp:162
virtual void finalize()
release memory that is not neeeded (empty sub-graphs)
Definition: graph_m4ri.cpp:199
Rect rect(Orientation o1, Orientation o2) const
range covered by this graph
Definition: graph_m4ri.cpp:94
static bool is_adjacent_to(const Graph &A, const Graph &B)
Definition: graph_m4ri.cpp:775
Graph & operator=(const Graph &)
Definition: graph_boost.cpp:64
this graph is the result of a MERGE operation (transitive closure) on three graphs
Definition: graph_m4ri.h:94
double max(double a, double b)
maximum function with checks for NAN
Definition: numeric.h:233
enum frechet::reach::Graph::Origin::@7 operation
IndexRange mask[2]
Definition: graph_m4ri.h:72
virtual int foundDiagonalElement() const
Definition: graph_m4ri.cpp:450
bool operator==(const Graph &that) const
comparison operator
Definition: graph_m4ri.cpp:240