MFEM  v4.2.0
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
gecko.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2020, Lawrence Livermore National Security, LLC. Produced
2 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3 // LICENSE and NOTICE for details. LLNL-CODE-806117.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability visit https://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the BSD-3 license. We welcome feedback and contributions, see file
10 // CONTRIBUTING.md for details.
11 
12 #include "gecko.hpp"
13 
14 // This file collects the sources of the Gecko library as a single module.
15 // The original library can be found at https://github.com/LLNL/gecko
16 // Used here with permission.
17 
18 // ------------------------------------------------------------------------------
19 
20 // BSD 3-Clause License
21 //
22 // Copyright (c) 2019-2020, Lawrence Livermore National Security, LLC
23 // All rights reserved.
24 //
25 // Redistribution and use in source and binary forms, with or without
26 // modification, are permitted provided that the following conditions are met:
27 //
28 // * Redistributions of source code must retain the above copyright notice, this
29 // list of conditions and the following disclaimer.
30 //
31 // * Redistributions in binary form must reproduce the above copyright notice,
32 // this list of conditions and the following disclaimer in the documentation
33 // and/or other materials provided with the distribution.
34 //
35 // * Neither the name of the copyright holder nor the names of its
36 // contributors may be used to endorse or promote products derived from
37 // this software without specific prior written permission.
38 //
39 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
40 // AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
41 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
42 // DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
43 // FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
44 // DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
45 // SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
46 // CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
47 // OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
48 // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
49 
50 // ------------------------------------------------------------------------------
51 
52 // Copyright (c) 2019-2020, Lawrence Livermore National Security, LLC and other
53 // gecko project contributors. See the above license for details.
54 // SPDX-License-Identifier: BSD-3-Clause
55 // LLNL-CODE-800597
56 
57 #include <algorithm>
58 #include <functional>
59 #include <cstdlib>
60 #include <cstddef>
61 #include <iomanip>
62 #include <iostream>
63 #include <sstream>
64 #include <stdexcept>
65 #include <map>
66 
67 // ----- options.h --------------------------------------------------------------
68 
69 // ratio of max to min weight for aggregation
70 #ifndef GECKO_PART_FRAC
71 #define GECKO_PART_FRAC 4
72 #endif
73 
74 // number of compatible relaxation sweeps
75 #ifndef GECKO_CR_SWEEPS
76 #define GECKO_CR_SWEEPS 1
77 #endif
78 
79 // number of Gauss-Seidel relaxation sweeps
80 #ifndef GECKO_GS_SWEEPS
81 #define GECKO_GS_SWEEPS 1
82 #endif
83 
84 // max number of nodes in subgraph
85 #ifndef GECKO_WINDOW_MAX
86 #define GECKO_WINDOW_MAX 16
87 #endif
88 
89 // use adjacency list (1) or adjacency matrix (0)
90 #ifndef GECKO_WITH_ADJLIST
91 #define GECKO_WITH_ADJLIST 0
92 #endif
93 
94 // use nonrecursive permutation algorithm
95 #ifndef GECKO_WITH_NONRECURSIVE
96 #define GECKO_WITH_NONRECURSIVE 0
97 #endif
98 
99 // use double-precision computations
100 #ifndef GECKO_WITH_DOUBLE_PRECISION
101 #define GECKO_WITH_DOUBLE_PRECISION 0
102 #endif
103 
104 // ----- heap.h -----------------------------------------------------------------
105 
106 template <
107  typename T, // data type
108  typename P, // priority type
109  class C = std::less<P>, // comparator for priorities
110  class M = std::map<T, unsigned int> // maps type T to unsigned integer
111  >
112 class DynamicHeap
113 {
114 public:
115  DynamicHeap(size_t count = 0);
116  ~DynamicHeap() {}
117  void insert(T data, P priority);
118  void update(T data, P priority);
119  bool top(T& data);
120  bool top(T& data, P& priority);
121  bool pop();
122  bool extract(T& data);
123  bool extract(T& data, P& priority);
124  bool erase(T data);
125  bool find(T data) const;
126  bool find(T data, P& priority) const;
127  bool empty() const { return heap.empty(); }
128  size_t size() const { return heap.size(); }
129 private:
130  struct HeapEntry
131  {
132  HeapEntry(P p, T d) : priority(p), data(d) {}
133  P priority;
134  T data;
135  };
136  std::vector<HeapEntry> heap;
137  M index;
138  C lower;
139  void ascend(unsigned int i);
140  void descend(unsigned int i);
141  void swap(unsigned int i, unsigned int j);
142  bool ordered(unsigned int i, unsigned int j) const
143  {
144  return !lower(heap[i].priority, heap[j].priority);
145  }
146  unsigned int parent(unsigned int i) const { return (i - 1) / 2; }
147  unsigned int left(unsigned int i) const { return 2 * i + 1; }
148  unsigned int right(unsigned int i) const { return 2 * i + 2; }
149 };
150 
151 template < typename T, typename P, class C, class M >
152 DynamicHeap<T, P, C, M>::DynamicHeap(size_t count)
153 {
154  heap.reserve(count);
155 }
156 
157 template < typename T, typename P, class C, class M >
158 void
159 DynamicHeap<T, P, C, M>::insert(T data, P priority)
160 {
161  if (index.find(data) != index.end())
162  {
163  update(data, priority);
164  }
165  else
166  {
167  unsigned int i = (unsigned int)heap.size();
168  heap.push_back(HeapEntry(priority, data));
169  ascend(i);
170  }
171 }
172 
173 template < typename T, typename P, class C, class M >
174 void
175 DynamicHeap<T, P, C, M>::update(T data, P priority)
176 {
177  unsigned int i = index[data];
178  heap[i].priority = priority;
179  ascend(i);
180  descend(i);
181 }
182 
183 template < typename T, typename P, class C, class M >
184 bool
185 DynamicHeap<T, P, C, M>::top(T& data)
186 {
187  if (!heap.empty())
188  {
189  data = heap[0].data;
190  return true;
191  }
192  else
193  {
194  return false;
195  }
196 }
197 
198 template < typename T, typename P, class C, class M >
199 bool
200 DynamicHeap<T, P, C, M>::top(T& data, P& priority)
201 {
202  if (!heap.empty())
203  {
204  data = heap[0].data;
205  priority = heap[0].priority;
206  return true;
207  }
208  else
209  {
210  return false;
211  }
212 }
213 
214 template < typename T, typename P, class C, class M >
215 bool
216 DynamicHeap<T, P, C, M>::pop()
217 {
218  if (!heap.empty())
219  {
220  T data = heap[0].data;
221  swap(0, (unsigned int)heap.size() - 1);
222  index.erase(data);
223  heap.pop_back();
224  if (!heap.empty())
225  {
226  descend(0);
227  }
228  return true;
229  }
230  else
231  {
232  return false;
233  }
234 }
235 
236 template < typename T, typename P, class C, class M >
237 bool
238 DynamicHeap<T, P, C, M>::extract(T& data)
239 {
240  if (!heap.empty())
241  {
242  data = heap[0].data;
243  return pop();
244  }
245  else
246  {
247  return false;
248  }
249 }
250 
251 template < typename T, typename P, class C, class M >
252 bool
253 DynamicHeap<T, P, C, M>::extract(T& data, P& priority)
254 {
255  if (!heap.empty())
256  {
257  data = heap[0].data;
258  priority = heap[0].priority;
259  return pop();
260  }
261  else
262  {
263  return false;
264  }
265 }
266 
267 template < typename T, typename P, class C, class M >
268 bool
269 DynamicHeap<T, P, C, M>::erase(T data)
270 {
271  if (index.find(data) == index.end())
272  {
273  return false;
274  }
275  unsigned int i = index[data];
276  swap(i, heap.size() - 1);
277  index.erase(data);
278  heap.pop_back();
279  if (i < heap.size())
280  {
281  ascend(i);
282  descend(i);
283  }
284  return true;
285 }
286 
287 template < typename T, typename P, class C, class M >
288 bool
289 DynamicHeap<T, P, C, M>::find(T data) const
290 {
291  return index.find(data) != index.end();
292 }
293 
294 template < typename T, typename P, class C, class M >
295 bool
296 DynamicHeap<T, P, C, M>::find(T data, P& priority) const
297 {
298  typename M::const_iterator p;
299  if ((p = index.find(data)) == index.end())
300  {
301  return false;
302  }
303  unsigned int i = p->second;
304  priority = heap[i].priority;
305  return true;
306 }
307 
308 template < typename T, typename P, class C, class M >
309 void
310 DynamicHeap<T, P, C, M>::ascend(unsigned int i)
311 {
312  for (unsigned int j; i && !ordered(j = parent(i), i); i = j)
313  {
314  swap(i, j);
315  }
316  index[heap[i].data] = i;
317 }
318 
319 template < typename T, typename P, class C, class M >
320 void
321 DynamicHeap<T, P, C, M>::descend(unsigned int i)
322 {
323  for (unsigned int j, k;
324  (j = ((k = left(i)) < heap.size() && !ordered(i, k) ? k : i),
325  j = ((k = right(i)) < heap.size() && !ordered(j, k) ? k : j)) != i;
326  i = j)
327  {
328  swap(i, j);
329  }
330  index[heap[i].data] = i;
331 }
332 
333 template < typename T, typename P, class C, class M >
334 void
335 DynamicHeap<T, P, C, M>::swap(unsigned int i, unsigned int j)
336 {
337  std::swap(heap[i], heap[j]);
338  index[heap[i].data] = i;
339 }
340 
341 // ----- subgraph.h -------------------------------------------------------------
342 
343 namespace Gecko
344 {
345 
346 // Node in a subgraph.
347 class Subnode
348 {
349 public:
350  typedef unsigned char Index;
351  Float pos; // node position
352  WeightedSum cost; // external cost at this position
353 };
354 
355 class Subgraph
356 {
357 public:
358  Subgraph(Graph* g, uint n);
359  ~Subgraph() { delete[] cache; }
360  void optimize(uint k);
361 
362 private:
363  Graph* const g; // full graph
364  const uint n; // number of subgraph nodes
365  Functional* const f; // ordering functional
366  WeightedSum min; // minimum cost so far
367  Subnode::Index best[GECKO_WINDOW_MAX]; // best permutation so far
368  Subnode::Index perm[GECKO_WINDOW_MAX]; // current permutation
369  const Subnode* node[GECKO_WINDOW_MAX]; // pointers to precomputed nodes
370  Subnode* cache; // precomputed node positions and costs
371 #if GECKO_WITH_ADJLIST
372  Subnode::Index
373  adj[GECKO_WINDOW_MAX][GECKO_WINDOW_MAX]; // internal adjacency list
374 #else
375  uint adj[GECKO_WINDOW_MAX]; // internal adjacency matrix
376 #endif
377  Float weight[GECKO_WINDOW_MAX][GECKO_WINDOW_MAX]; // internal arc weights
378  WeightedSum cost(uint k) const;
379  void swap(uint k);
380  void swap(uint k, uint l);
381  void optimize(WeightedSum c, uint i);
382 };
383 
384 }
385 
386 // ----- subgraph.cpp -----------------------------------------------------------
387 
388 using namespace Gecko;
389 
390 // Constructor.
391 Subgraph::Subgraph(Graph* g, uint n) : g(g), n(n), f(g->functional)
392 {
393  if (n > GECKO_WINDOW_MAX)
394  {
395  throw std::out_of_range("optimization window too large");
396  }
397  cache = new Subnode[n << n];
398 }
399 
400 // Cost of k'th node's edges to external nodes and nodes at {k+1, ..., n-1}.
402 Subgraph::cost(uint k) const
403 {
404  Subnode::Index i = perm[k];
405  WeightedSum c = node[i]->cost;
406  Float p = node[i]->pos;
407 #if GECKO_WITH_ADJLIST
408  for (k = 0; adj[i][k] != i; k++)
409  {
410  Subnode::Index j = adj[i][k];
411  Float l = node[j]->pos - p;
412  if (l > 0)
413  {
414  Float w = weight[i][k];
415  f->accumulate(c, WeightedValue(l, w));
416  }
417  }
418 #else
419  uint m = adj[i];
420  while (++k < n)
421  {
422  Subnode::Index j = perm[k];
423  if (m & (1u << j))
424  {
425  Float l = node[j]->pos - p;
426  Float w = weight[i][j];
427  f->accumulate(c, WeightedValue(l, w));
428  }
429  }
430 #endif
431  return c;
432 }
433 
434 // Swap the two nodes in positions k and k + 1.
435 void
436 Subgraph::swap(uint k)
437 {
438  uint l = k + 1;
439  Subnode::Index i = perm[k];
440  Subnode::Index j = perm[l];
441  perm[k] = j;
442  perm[l] = i;
443  node[i] -= ptrdiff_t(1) << j;
444  node[j] += ptrdiff_t(1) << i;
445 }
446 
447 // Swap the two nodes in positions k and l, k <= l.
448 void
449 Subgraph::swap(uint k, uint l)
450 {
451  Subnode::Index i = perm[k];
452  Subnode::Index j = perm[l];
453  perm[k] = j;
454  perm[l] = i;
455  // Update node positions.
456  uint m = 0;
457  while (++k < l)
458  {
459  Subnode::Index h = perm[k];
460  node[h] += ptrdiff_t(1) << i;
461  node[h] -= ptrdiff_t(1) << j;
462  m += 1u << h;
463  }
464  node[i] -= (1u << j) + m;
465  node[j] += (1u << i) + m;
466 }
467 
468 #if GECKO_WITH_NONRECURSIVE
469 // Evaluate all permutations generated by Heap's nonrecursive algorithm.
470 void
471 Subgraph::optimize(WeightedSum, uint)
472 {
473  WeightedSum c[GECKO_WINDOW_MAX + 1];
474  uint j[GECKO_WINDOW_MAX + 1];
475  j[n] = 1;
476  c[n] = 0;
477  uint i = n;
478  do
479  {
480  i--;
481  j[i] = i;
482  loop:
483  c[i] = f->sum(c[i + 1], cost(i));
484  }
485  while (i);
486  if (f->less(c[0], min))
487  {
488  min = c[0];
489  for (uint k = 0; k < n; k++)
490  {
491  best[k] = perm[k];
492  }
493  }
494  do
495  {
496  if (++i == n)
497  {
498  return;
499  }
500  swap(i & 1 ? i - j[i] : 0, i);
501  }
502  while (!j[i]--);
503  goto loop;
504 }
505 #else
506 // Apply branch-and-bound to permutations generated by Heap's algorithm.
507 void
508 Subgraph::optimize(WeightedSum c, uint i)
509 {
510  i--;
511  if (f->less(c, min))
512  {
513  if (i)
514  {
515  uint j = i;
516  do
517  {
518  optimize(f->sum(c, cost(i)), i);
519  swap(i & 1 ? i - j : 0, i);
520  }
521  while (j--);
522  }
523  else
524  {
525  f->accumulate(c, cost(0));
526  if (f->less(c, min))
527  {
528  min = c;
529  for (uint j = 0; j < n; j++)
530  {
531  best[j] = perm[j];
532  }
533  }
534  }
535  }
536  else if (i & 1)
537  do { swap(--i); }
538  while (i);
539 }
540 #endif
541 
542 // Optimize layout of nodes {p, ..., p + n - 1}.
543 void
544 Subgraph::optimize(uint p)
545 {
546  // Initialize subgraph.
547  const Float q = g->node[g->perm[p]].pos - g->node[g->perm[p]].hlen;
548  min = WeightedSum(GECKO_FLOAT_MAX, 1);
549  for (Subnode::Index k = 0; k < n; k++)
550  {
551  best[k] = perm[k] = k;
552  Node::Index i = g->perm[p + k];
553  // Copy i's outgoing arcs. We distinguish between internal
554  // and external arcs to nodes within and outside the subgraph,
555  // respectively.
556 #if GECKO_WITH_ADJLIST
557  uint m = 0;
558 #else
559  adj[k] = 0;
560 #endif
561  std::vector<Arc::Index> external;
562  for (Arc::Index a = g->node_begin(i); a < g->node_end(i); a++)
563  {
564  Node::Index j = g->adj[a];
565  Subnode::Index l;
566  for (l = 0; l < n && g->perm[p + l] != j; l++);
567  if (l == n)
568  {
569  external.push_back(a);
570  }
571  else
572  {
573  // Copy internal arc to subgraph.
574 #if GECKO_WITH_ADJLIST
575  adj[k][m] = l;
576  weight[k][m] = g->weight[a];
577  m++;
578 #else
579  adj[k] += 1u << l;
580  weight[k][l] = g->weight[a];
581 #endif
582  }
583  }
584 #if GECKO_WITH_ADJLIST
585  adj[k][m] = k;
586 #endif
587  // Precompute external costs associated with all possible positions
588  // of this node. Since node lengths can be arbitrary, there are as
589  // many as 2^(n-1) possible positions, each corresponding to an
590  // (n-1)-bit string that specifies whether the remaining n-1 nodes
591  // succeed this node or not. Caching the
592  // n
593  // 2^(n-1) n = sum k C(n, k) = A001787
594  // k=1
595  // external costs is exponentially cheaper than recomputing the
596  // n-1 n
597  // n! sum 1/k! = sum k! C(n, k) = A007526
598  // k=0 k=1
599  // costs associated with all permutations.
600  node[k] = cache + (k << n);
601  for (uint m = 0; m < (1u << n); m++)
602  if (!(m & (1u << k)))
603  {
604  Subnode* s = cache + (k << n) + m;
605  s->pos = q + g->node[i].hlen;
606  for (Subnode::Index l = 0; l < n; l++)
607  if (l != k && !(m & (1u << l)))
608  {
609  s->pos += 2 * g->node[g->perm[p + l]].hlen;
610  }
611  s->cost = g->cost(external, s->pos);
612  }
613  else
614  {
615  m += (1u << k) - 1;
616  }
617  node[k] += (1u << n) - (2u << k);
618  }
619 
620  // Find optimal permutation of the n nodes.
621  optimize(0, n);
622 
623  // Apply permutation to original graph.
624  for (uint i = 0; i < n; i++)
625  {
626  g->swap(p + i, p + best[i]);
627  for (uint j = i + 1; j < n; j++)
628  if (best[j] == i)
629  {
630  best[j] = best[i];
631  }
632  }
633 }
634 
635 // ----- graph.cpp --------------------------------------------------------------
636 
637 using namespace std;
638 using namespace Gecko;
639 
640 // Constructor.
641 void
642 Graph::init(uint nodes)
643 {
644  node.push_back(Node(-1, 0, 1, Node::null));
645  adj.push_back(Node::null);
646  weight.push_back(0);
647  bond.push_back(0);
648  while (nodes--)
649  {
650  insert_node();
651  }
652 }
653 
654 // Insert node.
656 Graph::insert_node(Float length)
657 {
658  Node::Index p = Node::Index(node.size());
659  perm.push_back(p);
660  node.push_back(Node(-1, length));
661  return p;
662 }
663 
664 // Return nodes adjacent to i.
665 std::vector<Node::Index>
666 Graph::node_neighbors(Node::Index i) const
667 {
668  std::vector<Node::Index> neighbor;
669  for (Arc::Index a = node_begin(i); a < node_end(i); a++)
670  {
671  neighbor.push_back(adj[a]);
672  }
673  return neighbor;
674 }
675 
676 // Insert directed edge (i, j).
678 Graph::insert_arc(Node::Index i, Node::Index j, Float w, Float b)
679 {
680  if (!i || !j || i == j || !(last_node <= i && i <= nodes()))
681  {
682  return Arc::null;
683  }
684  last_node = i;
685  for (Node::Index k = i - 1; node[k].arc == Arc::null; k--)
686  {
687  node[k].arc = Arc::Index(adj.size());
688  }
689  adj.push_back(j);
690  weight.push_back(w);
691  bond.push_back(b);
692  node[i].arc = Arc::Index(adj.size());
693  return Arc::Index(adj.size() - 1);
694 }
695 
696 // Remove arc a.
697 bool
698 Graph::remove_arc(Arc::Index a)
699 {
700  if (a == Arc::null)
701  {
702  return false;
703  }
704  Node::Index i = arc_source(a);
705  adj.erase(adj.begin() + a);
706  weight.erase(weight.begin() + a);
707  bond.erase(bond.begin() + a);
708  for (Node::Index k = i; k < node.size(); k++)
709  {
710  node[k].arc--;
711  }
712  return true;
713 }
714 
715 // Remove directed edge (i, j).
716 bool
717 Graph::remove_arc(Node::Index i, Node::Index j)
718 {
719  return remove_arc(arc_index(i, j));
720 }
721 
722 // Remove edge {i, j}.
723 bool
724 Graph::remove_edge(Node::Index i, Node::Index j)
725 {
726  bool success = remove_arc(i, j);
727  if (success)
728  {
729  success = remove_arc(j, i);
730  }
731  return success;
732 }
733 
734 // Index of arc (i, j) or null if not a valid arc.
736 Graph::arc_index(Node::Index i, Node::Index j) const
737 {
738  for (Arc::Index a = node_begin(i); a < node_end(i); a++)
739  if (adj[a] == j)
740  {
741  return a;
742  }
743  return Arc::null;
744 }
745 
746 // Return source node i in arc a = (i, j).
748 Graph::arc_source(Arc::Index a) const
749 {
750  Node::Index j = adj[a];
751  for (Arc::Index b = node_begin(j); b < node_end(j); b++)
752  {
753  Node::Index i = adj[b];
754  if (node_begin(i) <= a && a < node_end(i))
755  {
756  return i;
757  }
758  }
759  // should never get here
760  throw std::runtime_error("internal data structure corrupted");
761 }
762 
763 // Return reverse arc (j, i) of arc a = (i, j).
765 Graph::reverse_arc(Arc::Index a) const
766 {
767  Node::Index j = adj[a];
768  for (Arc::Index b = node_begin(j); b < node_end(j); b++)
769  {
770  Node::Index i = adj[b];
771  if (node_begin(i) <= a && a < node_end(i))
772  {
773  return b;
774  }
775  }
776  return Arc::null;
777 }
778 
779 // Return first directed arc if one exists or null otherwise.
781 Graph::directed() const
782 {
783  for (Node::Index i = 1; i < node.size(); i++)
784  for (Arc::Index a = node_begin(i); a < node_end(i); a++)
785  {
786  Node::Index j = adj[a];
787  if (!arc_index(j, i))
788  {
789  return a;
790  }
791  }
792  return Arc::null;
793 }
794 
795 // Add contribution of fine arc to coarse graph.
796 void
797 Graph::update(Node::Index i, Node::Index j, Float w, Float b)
798 {
799  Arc::Index a = arc_index(i, j);
800  if (a == Arc::null)
801  {
802  insert_arc(i, j, w, b);
803  }
804  else
805  {
806  weight[a] += w;
807  bond[a] += b;
808  }
809 }
810 
811 // Transfer contribution of fine arc a to coarse node p.
812 void
813 Graph::transfer(Graph* g, const vector<Float>& part, Node::Index p,
814  Arc::Index a, Float f) const
815 {
816  Float w = f * weight[a];
817  Float m = f * bond[a];
818  Node::Index j = arc_target(a);
819  Node::Index q = node[j].parent;
820  if (q == Node::null)
821  {
822  for (Arc::Index b = node_begin(j); b < node_end(j); b++)
823  if (part[b] > 0)
824  {
825  q = node[adj[b]].parent;
826  if (q != p)
827  {
828  g->update(p, q, w * part[b], m * part[b]);
829  }
830  }
831  }
832  else
833  {
834  g->update(p, q, w, m);
835  }
836 }
837 
838 // Compute cost of a subset of arcs incident on node placed at pos.
840 Graph::cost(const vector<Arc::Index>& subset, Float pos) const
841 {
842  WeightedSum c;
843  for (Arc::ConstPtr ap = subset.begin(); ap != subset.end(); ap++)
844  {
845  Arc::Index a = *ap;
846  Node::Index j = arc_target(a);
847  Float l = fabs(node[j].pos - pos);
848  Float w = weight[a];
849  functional->accumulate(c, WeightedValue(l, w));
850  }
851  return c;
852 }
853 
854 // Compute cost of graph layout.
855 Float
856 Graph::cost() const
857 {
858  if (edges())
859  {
860  WeightedSum c;
861  Node::Index i = 1;
862  for (Arc::Index a = 1; a < adj.size(); a++)
863  {
864  while (node_end(i) <= a)
865  {
866  i++;
867  }
868  Node::Index j = arc_target(a);
869  Float l = length(i, j);
870  Float w = weight[a];
871  functional->accumulate(c, WeightedValue(l, w));
872  }
873  return functional->mean(c);
874  }
875  else
876  {
877  return Float(0);
878  }
879 }
880 
881 // Swap the two nodes in positions k and l, k <= l.
882 void
883 Graph::swap(uint k, uint l)
884 {
885  Node::Index i = perm[k];
886  perm[k] = perm[l];
887  perm[l] = i;
888  Float p = node[i].pos - node[i].hlen;
889  do
890  {
891  i = perm[k];
892  p += node[i].hlen;
893  node[i].pos = p;
894  p += node[i].hlen;
895  }
896  while (k++ != l);
897 }
898 
899 // Optimize continuous position of a single node.
900 Float
901 Graph::optimal(Node::Index i) const
902 {
903  vector<WeightedValue> v;
904  for (Arc::Index a = node_begin(i); a < node_end(i); a++)
905  {
906  Node::Index j = adj[a];
907  if (placed(j))
908  {
909  v.push_back(WeightedValue(node[j].pos, weight[a]));
910  }
911  }
912  return v.empty() ? -1 : functional->optimum(v);
913 }
914 
915 // Compute coarse graph with roughly half the number of nodes.
916 Graph*
917 Graph::coarsen()
918 {
919  progress->beginphase(this, string("coarse"));
920  Graph* g = new Graph(0, level - 1);
921  g->functional = functional;
922  g->progress = progress;
923 
924  // Compute importance of nodes in fine graph.
925  DynamicHeap<Node::Index, Float> heap;
926  for (Node::Index i = 1; i < node.size(); i++)
927  {
928  node[i].parent = Node::null;
929  Float w = 0;
930  for (Arc::Index a = node_begin(i); a < node_end(i); a++)
931  {
932  w += bond[a];
933  }
934  heap.insert(i, w);
935  }
936 
937  // Select set of important nodes from fine graph that will remain in
938  // coarse graph.
939  vector<Node::Index> child(1, Node::null);
940  while (!heap.empty())
941  {
942  Node::Index i;
943  Float w = 0;
944  heap.extract(i, w);
945  if (w < 0)
946  {
947  break;
948  }
949  child.push_back(i);
950  node[i].parent = g->insert_node(2 * node[i].hlen);
951 
952  // Reduce importance of neighbors.
953  for (Arc::Index a = node_begin(i); a < node_end(i); a++)
954  {
955  Node::Index j = adj[a];
956  if (heap.find(j, w))
957  {
958  heap.update(j, w - 2 * bond[a]);
959  }
960  }
961  }
962 
963  // Assign parts of remaining nodes to aggregates.
964  vector<Float> part = bond;
965  for (Node::Index i = 1; i < node.size(); i++)
966  if (!persistent(i))
967  {
968  // Find all connections to coarse nodes.
969  Float w = 0;
970  Float max = 0;
971  for (Arc::Index a = node_begin(i); a < node_end(i); a++)
972  {
973  Node::Index j = adj[a];
974  if (persistent(j))
975  {
976  w += part[a];
977  if (max < part[a])
978  {
979  max = part[a];
980  }
981  }
982  else
983  {
984  part[a] = -1;
985  }
986  }
987  max /= GECKO_PART_FRAC;
988 
989  // Weed out insignificant connections.
990  for (Arc::Index a = node_begin(i); a < node_end(i); a++)
991  if (0 < part[a] && part[a] < max)
992  {
993  w -= part[a];
994  part[a] = -1;
995  }
996 
997  // Compute node fractions (interpolation matrix) and assign
998  // partial nodes to aggregates.
999  for (Arc::Index a = node_begin(i); a < node_end(i); a++)
1000  if (part[a] > 0)
1001  {
1002  part[a] /= w;
1003  Node::Index p = node[adj[a]].parent;
1004  g->node[p].hlen += part[a] * node[i].hlen;
1005  }
1006  }
1007 
1008  // Transfer arcs to coarse graph.
1009  for (Node::Index p = 1; p < g->node.size(); p++)
1010  {
1011  Node::Index i = child[p];
1012  for (Arc::Index a = node_begin(i); a < node_end(i); a++)
1013  {
1014  transfer(g, part, p, a);
1015  Node::Index j = adj[a];
1016  if (!persistent(j))
1017  {
1018  Arc::Index b = arc_index(j, i);
1019  if (part[b] > 0)
1020  for (Arc::Index c = node_begin(j); c < node_end(j); c++)
1021  {
1022  Node::Index k = adj[c];
1023  if (k != i)
1024  {
1025  transfer(g, part, p, c, part[b]);
1026  }
1027  }
1028  }
1029  }
1030  }
1031 
1032 #if DEBUG
1033  if (g->directed())
1034  {
1035  throw runtime_error("directed edge found");
1036  }
1037 #endif
1038 
1039  // Free memory.
1040  vector<Float> t = bond;
1041  bond.swap(t);
1042 
1043  progress->endphase(this, false);
1044 
1045  return g;
1046 }
1047 
1048 // Order nodes according to coarsened graph layout.
1049 void
1050 Graph::refine(const Graph* graph)
1051 {
1052  progress->beginphase(this, string("refine"));
1053 
1054  // Place persistent nodes.
1055  DynamicHeap<Node::Index, Float> heap;
1056  for (Node::Index i = 1; i < node.size(); i++)
1057  if (persistent(i))
1058  {
1059  Node::Index p = node[i].parent;
1060  node[i].pos = graph->node[p].pos;
1061  }
1062  else
1063  {
1064  node[i].pos = -1;
1065  Float w = 0;
1066  for (Arc::Index a = node_begin(i); a < node_end(i); a++)
1067  {
1068  Node::Index j = adj[a];
1069  if (persistent(j))
1070  {
1071  w += weight[a];
1072  }
1073  }
1074  heap.insert(i, w);
1075  }
1076 
1077  // Place remaining nodes in order of decreasing connectivity with
1078  // already placed nodes.
1079  while (!heap.empty())
1080  {
1081  Node::Index i = 0;
1082  heap.extract(i);
1083  node[i].pos = optimal(i);
1084  for (Arc::Index a = node_begin(i); a < node_end(i); a++)
1085  {
1086  Node::Index j = adj[a];
1087  Float w;
1088  if (heap.find(j, w))
1089  {
1090  heap.update(j, w + weight[a]);
1091  }
1092  }
1093  }
1094 
1095  place(true);
1096  progress->endphase(this, true);
1097 }
1098 
1099 // Perform m sweeps of compatible or Gauss-Seidel relaxation.
1100 void
1101 Graph::relax(bool compatible, uint m)
1102 {
1103  progress->beginphase(this, compatible ? string("crelax") : string("frelax"));
1104  while (m--)
1105  for (uint k = 0; k < perm.size() && !progress->quit(); k++)
1106  {
1107  Node::Index i = perm[k];
1108  if (!compatible || !persistent(i))
1109  {
1110  node[i].pos = optimal(i);
1111  }
1112  }
1113  place(true);
1114  progress->endphase(this, true);
1115 }
1116 
1117 // Optimize successive n-node subgraphs.
1118 void
1119 Graph::optimize(uint n)
1120 {
1121  if (n > perm.size())
1122  {
1123  n = uint(perm.size());
1124  }
1125  ostringstream count;
1126  count << setw(2) << n;
1127  progress->beginphase(this, string("perm") + count.str());
1128  Subgraph* subgraph = new Subgraph(this, n);
1129  for (uint k = 0; k <= perm.size() - n && !progress->quit(); k++)
1130  {
1131  subgraph->optimize(k);
1132  }
1133  delete subgraph;
1134  progress->endphase(this, true);
1135 }
1136 
1137 // Place all nodes according to their positions.
1138 void
1139 Graph::place(bool sort)
1140 {
1141  place(sort, 0, uint(perm.size()));
1142 }
1143 
1144 // Place nodes {k, ..., k + n - 1} according to their positions.
1145 void
1146 Graph::place(bool sort, uint k, uint n)
1147 {
1148  // Place nodes.
1149  if (sort)
1150  {
1151  stable_sort(perm.begin() + k, perm.begin() + k + n,
1152  Node::Comparator(node.begin()));
1153  }
1154 
1155  // Assign node positions according to permutation.
1156  for (Float p = k ? node[perm[k - 1]].pos + node[perm[k - 1]].hlen : 0; n--;
1157  k++)
1158  {
1159  Node::Index i = perm[k];
1160  p += node[i].hlen;
1161  node[i].pos = p;
1162  p += node[i].hlen;
1163  }
1164 }
1165 
1166 // Perform one V-cycle.
1167 void
1168 Graph::vcycle(uint n, uint work)
1169 {
1170  if (n < nodes() && nodes() < edges() && level && !progress->quit())
1171  {
1172  Graph* graph = coarsen();
1173  graph->vcycle(n, work + edges());
1174  refine(graph);
1175  delete graph;
1176  }
1177  else
1178  {
1179  place();
1180  }
1181  if (edges())
1182  {
1183  relax(true, GECKO_CR_SWEEPS);
1184  relax(false, GECKO_GS_SWEEPS);
1185  for (uint w = edges(); w * (n + 1) < work; w *= ++n);
1186  n = std::min(n, uint(GECKO_WINDOW_MAX));
1187  if (n)
1188  {
1189  optimize(n);
1190  }
1191  }
1192 }
1193 
1194 // Custom random-number generator for reproducibility.
1195 // LCG from doi:10.1090/S0025-5718-99-00996-5.
1196 uint
1197 Graph::random(uint seed)
1198 {
1199  static uint state = 1;
1200  state = (seed ? seed : 0x1ed0675 * state + 0xa14f);
1201  return state;
1202 }
1203 
1204 // Generate a random permutation of the nodes.
1205 void
1206 Graph::shuffle(uint seed)
1207 {
1208  random(seed);
1209  for (uint k = 0; k < perm.size(); k++)
1210  {
1211  uint r = random() >> 8;
1212  uint l = k + r % (uint(perm.size()) - k);
1213  std::swap(perm[k], perm[l]);
1214  }
1215  place();
1216 }
1217 
1218 // Recompute bonds for k'th V-cycle.
1219 void
1220 Graph::reweight(uint k)
1221 {
1222  bond.resize(weight.size());
1223  for (Arc::Index a = 1; a < adj.size(); a++)
1224  {
1225  bond[a] = functional->bond(weight[a], length(a), k);
1226  }
1227 }
1228 
1229 // Linearly order graph.
1230 void
1231 Graph::order(Functional* functional, uint iterations, uint window, uint period,
1232  uint seed, Progress* progress)
1233 {
1234  // Initialize graph.
1235  this->functional = functional;
1236  progress = this->progress = progress ? progress : new Progress;
1237  for (level = 0; (1u << level) < nodes(); level++);
1238  place();
1239  Float mincost = cost();
1240  vector<Node::Index> minperm = perm;
1241  if (seed)
1242  {
1243  shuffle(seed);
1244  }
1245 
1246  progress->beginorder(this, mincost);
1247  if (edges())
1248  {
1249  // Perform specified number of V-cycles.
1250  for (uint k = 1; k <= iterations && !progress->quit(); k++)
1251  {
1252  progress->beginiter(this, k, iterations, window);
1253  reweight(k);
1254  vcycle(window);
1255  Float c = cost();
1256  if (c < mincost)
1257  {
1258  mincost = c;
1259  minperm = perm;
1260  }
1261  progress->enditer(this, mincost, c);
1262  if (period && !(k % period))
1263  {
1264  window++;
1265  }
1266  }
1267  perm = minperm;
1268  place();
1269  }
1270  progress->endorder(this, mincost);
1271 
1272  if (!progress)
1273  {
1274  delete this->progress;
1275  this->progress = 0;
1276  }
1277 }
Arc::Index node_begin(Node::Index i) const
Definition: gecko.hpp:635
std::vector< Index >::const_iterator ConstPtr
Definition: gecko.hpp:587
std::vector< Float > weight
Definition: gecko.hpp:735
unsigned int uint
Definition: gecko.hpp:204
virtual void beginiter(const Graph *, uint, uint, uint) const
Definition: gecko.hpp:566
std::vector< Node::Index > adj
Definition: gecko.hpp:734
Arc::Index node_end(Node::Index i) const
Definition: gecko.hpp:636
virtual void endorder(const Graph *, Float) const
Definition: gecko.hpp:565
std::vector< Node > node
Definition: gecko.hpp:733
virtual void enditer(const Graph *, Float, Float) const
Definition: gecko.hpp:568
uint Index
Definition: gecko.hpp:595
std::vector< Node::Index > perm
Definition: gecko.hpp:732
double b
Definition: lissajous.cpp:42
Float cost() const
Definition: gecko.cpp:856
double Float
Definition: gecko.hpp:208
Arc::Index directed() const
Definition: gecko.cpp:781
Functional * functional
Definition: gecko.hpp:730
uint Index
Definition: gecko.hpp:586
void vcycle(uint n, uint work=0)
Definition: gecko.cpp:1168
double a
Definition: lissajous.cpp:41
virtual void beginorder(const Graph *, Float) const
Definition: gecko.hpp:564
int index(int i, int j, int nx, int ny)
Definition: life.cpp:241
Node::Index insert_node(Float length=1)
Definition: gecko.cpp:656
virtual bool quit() const
Definition: gecko.hpp:572
Progress * progress
Definition: gecko.hpp:731
double f(const Vector &p)