text: // Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
  // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
  // LICENSE and NOTICE for details. LLNL-CODE-806117.
  //
  // This file is part of the MFEM library. For more information and source code
  // availability visit https://mfem.org.
  //
  // MFEM is free software; you can redistribute it and/or modify it under the
  // terms of the BSD-3 license. We welcome feedback and contributions, see file
  // CONTRIBUTING.md for details.
  
  #include "gecko.hpp"
  
  // This file collects the sources of the Gecko library as a single module.
  // The original library can be found at https://github.com/LLNL/gecko
  // Used here with permission.
  
  // ------------------------------------------------------------------------------
  
  // BSD 3-Clause License
  //
  // Copyright (c) 2019-2021, Lawrence Livermore National Security, LLC
  // All rights reserved.
  //
  // Redistribution and use in source and binary forms, with or without
  // modification, are permitted provided that the following conditions are met:
  //
  // * Redistributions of source code must retain the above copyright notice, this
  //   list of conditions and the following disclaimer.
  //
  // * Redistributions in binary form must reproduce the above copyright notice,
  //   this list of conditions and the following disclaimer in the documentation
  //   and/or other materials provided with the distribution.
  //
  // * Neither the name of the copyright holder nor the names of its
  //   contributors may be used to endorse or promote products derived from
  //   this software without specific prior written permission.
  //
  // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
  // AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
  // DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
  // FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
  // DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
  // SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
  // CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
  // OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
  // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  
  // ------------------------------------------------------------------------------
  
  // Copyright (c) 2019-2021, Lawrence Livermore National Security, LLC and other
  // gecko project contributors. See the above license for details.
  // SPDX-License-Identifier: BSD-3-Clause
  // LLNL-CODE-800597
  
  #include <algorithm>
  #include <functional>
  #include <cstdlib>
  #include <cstddef>
  #include <iomanip>
  #include <iostream>
  #include <sstream>
  #include <stdexcept>
  #include <map>
  
  // ----- options.h --------------------------------------------------------------
  
  // ratio of max to min weight for aggregation
  #ifndef GECKO_PART_FRAC
  #define GECKO_PART_FRAC 4
  #endif
  
  // number of compatible relaxation sweeps
  #ifndef GECKO_CR_SWEEPS
  #define GECKO_CR_SWEEPS 1
  #endif
  
  // number of Gauss-Seidel relaxation sweeps
  #ifndef GECKO_GS_SWEEPS
  #define GECKO_GS_SWEEPS 1
  #endif
  
  // max number of nodes in subgraph
  #ifndef GECKO_WINDOW_MAX
  #define GECKO_WINDOW_MAX 16
  #endif
  
  // use adjacency list (1) or adjacency matrix (0)
  #ifndef GECKO_WITH_ADJLIST
  #define GECKO_WITH_ADJLIST 0
  #endif
  
  // use nonrecursive permutation algorithm
  #ifndef GECKO_WITH_NONRECURSIVE
  #define GECKO_WITH_NONRECURSIVE 0
  #endif
  
  // use double-precision computations
  #ifndef GECKO_WITH_DOUBLE_PRECISION
  #define GECKO_WITH_DOUBLE_PRECISION 0
  #endif
  
  // ----- heap.h -----------------------------------------------------------------
  
  template <
     typename T,                             // data type
     typename P,                             // priority type
     class    C = std::less<P>,              // comparator for priorities
     class    M = std::map<T, unsigned int>  // maps type T to unsigned integer
     >
  class DynamicHeap
  {
  public:
     DynamicHeap(size_t count = 0);
     ~DynamicHeap() {}
     void insert(T data, P priority);
     void update(T data, P priority);
     bool top(T& data);
     bool top(T& data, P& priority);
     bool pop();
     bool extract(T& data);
     bool extract(T& data, P& priority);
     bool erase(T data);
     bool find(T data) const;
     bool find(T data, P& priority) const;
     bool empty() const { return heap.empty(); }
     size_t size() const { return heap.size(); }
  private:
     struct HeapEntry
     {
        HeapEntry(P p, T d) : priority(p), data(d) {}
        P priority;
        T data;
     };
     std::vector<HeapEntry> heap;
     M index;
     C lower;
     void ascend(unsigned int i);
     void descend(unsigned int i);
     void swap(unsigned int i, unsigned int j);
     bool ordered(unsigned int i, unsigned int j) const
     {
        return !lower(heap[i].priority, heap[j].priority);
     }
     unsigned int parent(unsigned int i) const { return (i - 1) / 2; }
     unsigned int left(unsigned int i) const { return 2 * i + 1; }
     unsigned int right(unsigned int i) const { return 2 * i + 2; }
  };
  
  template < typename T, typename P, class C, class M >
  DynamicHeap<T, P, C, M>::DynamicHeap(size_t count)
  {
     heap.reserve(count);
  }
  
  template < typename T, typename P, class C, class M >
  void
  DynamicHeap<T, P, C, M>::insert(T data, P priority)
  {
     if (index.find(data) != index.end())
     {
        update(data, priority);
     }
     else
     {
        unsigned int i = (unsigned int)heap.size();
        heap.push_back(HeapEntry(priority, data));
        ascend(i);
     }
  }
  
  template < typename T, typename P, class C, class M >
  void
  DynamicHeap<T, P, C, M>::update(T data, P priority)
  {
     unsigned int i = index[data];
     heap[i].priority = priority;
     ascend(i);
     descend(i);
  }
  
  template < typename T, typename P, class C, class M >
  bool
  DynamicHeap<T, P, C, M>::top(T& data)
  {
     if (!heap.empty())
     {
        data = heap[0].data;
        return true;
     }
     else
     {
        return false;
     }
  }
  
  template < typename T, typename P, class C, class M >
  bool
  DynamicHeap<T, P, C, M>::top(T& data, P& priority)
  {
     if (!heap.empty())
     {
        data = heap[0].data;
        priority = heap[0].priority;
        return true;
     }
     else
     {
        return false;
     }
  }
  
  template < typename T, typename P, class C, class M >
  bool
  DynamicHeap<T, P, C, M>::pop()
  {
     if (!heap.empty())
     {
        T data = heap[0].data;
        swap(0, (unsigned int)heap.size() - 1);
        index.erase(data);
        heap.pop_back();
        if (!heap.empty())
        {
           descend(0);
        }
        return true;
     }
     else
     {
        return false;
     }
  }
  
  template < typename T, typename P, class C, class M >
  bool
  DynamicHeap<T, P, C, M>::extract(T& data)
  {
     if (!heap.empty())
     {
        data = heap[0].data;
        return pop();
     }
     else
     {
        return false;
     }
  }
  
  template < typename T, typename P, class C, class M >
  bool
  DynamicHeap<T, P, C, M>::extract(T& data, P& priority)
  {
     if (!heap.empty())
     {
        data = heap[0].data;
        priority = heap[0].priority;
        return pop();
     }
     else
     {
        return false;
     }
  }
  
  template < typename T, typename P, class C, class M >
  bool
  DynamicHeap<T, P, C, M>::erase(T data)
  {
     if (index.find(data) == index.end())
     {
        return false;
     }
     unsigned int i = index[data];
     swap(i, heap.size() - 1);
     index.erase(data);
     heap.pop_back();
     if (i < heap.size())
     {
        ascend(i);
        descend(i);
     }
     return true;
  }
  
  template < typename T, typename P, class C, class M >
  bool
  DynamicHeap<T, P, C, M>::find(T data) const
  {
     return index.find(data) != index.end();
  }
  
  template < typename T, typename P, class C, class M >
  bool
  DynamicHeap<T, P, C, M>::find(T data, P& priority) const
  {
     typename M::const_iterator p;
     if ((p = index.find(data)) == index.end())
     {
        return false;
     }
     unsigned int i = p->second;
     priority = heap[i].priority;
     return true;
  }
  
  template < typename T, typename P, class C, class M >
  void
  DynamicHeap<T, P, C, M>::ascend(unsigned int i)
  {
     for (unsigned int j; i && !ordered(j = parent(i), i); i = j)
     {
        swap(i, j);
     }
     index[heap[i].data] = i;
  }
  
  template < typename T, typename P, class C, class M >
  void
  DynamicHeap<T, P, C, M>::descend(unsigned int i)
  {
     for (unsigned int j, k;
          (j = ((k =  left(i)) < heap.size() && !ordered(i, k) ? k : i),
           j = ((k = right(i)) < heap.size() && !ordered(j, k) ? k : j)) != i;
          i = j)
     {
        swap(i, j);
     }
     index[heap[i].data] = i;
  }
  
  template < typename T, typename P, class C, class M >
  void
  DynamicHeap<T, P, C, M>::swap(unsigned int i, unsigned int j)
  {
     std::swap(heap[i], heap[j]);
     index[heap[i].data] = i;
  }
  
  // ----- subgraph.h -------------------------------------------------------------
  
  namespace Gecko
  {
  
  // Node in a subgraph.
  class Subnode
  {
  public:
     typedef unsigned char Index;
     Float pos;        // node position
     WeightedSum cost; // external cost at this position
  };
  
  class Subgraph
  {
  public:
     Subgraph(Graph* g, uint n);
     ~Subgraph() { delete[] cache; }
     void optimize(uint k);
  
  private:
     Graph* const g;                        // full graph
     const uint n;                          // number of subgraph nodes
     Functional* const f;                   // ordering functional
     WeightedSum min;                       // minimum cost so far
     Subnode::Index best[GECKO_WINDOW_MAX]; // best permutation so far
     Subnode::Index perm[GECKO_WINDOW_MAX]; // current permutation
     const Subnode* node[GECKO_WINDOW_MAX]; // pointers to precomputed nodes
     Subnode* cache;                        // precomputed node positions and costs
  #if GECKO_WITH_ADJLIST
     Subnode::Index
     adj[GECKO_WINDOW_MAX][GECKO_WINDOW_MAX]; // internal adjacency list
  #else
     uint adj[GECKO_WINDOW_MAX];            // internal adjacency matrix
  #endif
     Float weight[GECKO_WINDOW_MAX][GECKO_WINDOW_MAX]; // internal arc weights
     WeightedSum cost(uint k) const;
     void swap(uint k);
     void swap(uint k, uint l);
     void optimize(WeightedSum c, uint i);
  };
  
  }
  
  // ----- subgraph.cpp -----------------------------------------------------------
  
  using namespace Gecko;
  
  // Constructor.
  Subgraph::Subgraph(Graph* g, uint n) : g(g), n(n), f(g->functional)
  {
     if (n > GECKO_WINDOW_MAX)
     {
        throw std::out_of_range("optimization window too large");
     }
     cache = new Subnode[n << n];
  }
  
  // Cost of k'th node's edges to external nodes and nodes at {k+1, ..., n-1}.
  WeightedSum
  Subgraph::cost(uint k) const
  {
     Subnode::Index i = perm[k];
     WeightedSum c = node[i]->cost;
     Float p = node[i]->pos;
  #if GECKO_WITH_ADJLIST
     for (k = 0; adj[i][k] != i; k++)
     {
        Subnode::Index j = adj[i][k];
        Float l = node[j]->pos - p;
        if (l > 0)
        {
           Float w = weight[i][k];
           f->accumulate(c, WeightedValue(l, w));
        }
     }
  #else
     uint m = adj[i];
     while (++k < n)
     {
        Subnode::Index j = perm[k];
        if (m & (1u << j))
        {
           Float l = node[j]->pos - p;
           Float w = weight[i][j];
           f->accumulate(c, WeightedValue(l, w));
        }
     }
  #endif
     return c;
  }
  
  // Swap the two nodes in positions k and k + 1.
  void
  Subgraph::swap(uint k)
  {
     uint l = k + 1;
     Subnode::Index i = perm[k];
     Subnode::Index j = perm[l];
     perm[k] = j;
     perm[l] = i;
     node[i] -= ptrdiff_t(1) << j;
     node[j] += ptrdiff_t(1) << i;
  }
  
  // Swap the two nodes in positions k and l, k <= l.
  void
  Subgraph::swap(uint k, uint l)
  {
     Subnode::Index i = perm[k];
     Subnode::Index j = perm[l];
     perm[k] = j;
     perm[l] = i;
     // Update node positions.
     uint m = 0;
     while (++k < l)
     {
        Subnode::Index h = perm[k];
        node[h] += ptrdiff_t(1) << i;
        node[h] -= ptrdiff_t(1) << j;
        m += 1u << h;
     }
     node[i] -= (1u << j) + m;
     node[j] += (1u << i) + m;
  }
  
  #if GECKO_WITH_NONRECURSIVE
  // Evaluate all permutations generated by Heap's nonrecursive algorithm.
  void
  Subgraph::optimize(WeightedSum, uint)
  {
     WeightedSum c[GECKO_WINDOW_MAX + 1];
     uint j[GECKO_WINDOW_MAX + 1];
     j[n] = 1;
     c[n] = 0;
     uint i = n;
     do
     {
        i--;
        j[i] = i;
     loop:
        c[i] = f->sum(c[i + 1], cost(i));
     }
     while (i);
     if (f->less(c[0], min))
     {
        min = c[0];
        for (uint k = 0; k < n; k++)
        {
           best[k] = perm[k];
        }
     }
     do
     {
        if (++i == n)
        {
           return;
        }
        swap(i & 1 ? i - j[i] : 0, i);
     }
     while (!j[i]--);
     goto loop;
  }
  #else
  // Apply branch-and-bound to permutations generated by Heap's algorithm.
  void
  Subgraph::optimize(WeightedSum c, uint i)
  {
     i--;
     if (f->less(c, min))
     {
        if (i)
        {
           uint j = i;
           do
           {
              optimize(f->sum(c, cost(i)), i);
              swap(i & 1 ? i - j : 0, i);
           }
           while (j--);
        }
        else
        {
           f->accumulate(c, cost(0));
           if (f->less(c, min))
           {
              min = c;
              for (uint j = 0; j < n; j++)
              {
                 best[j] = perm[j];
              }
           }
        }
     }
     else if (i & 1)
        do { swap(--i); }
        while (i);
  }
  #endif
  
  // Optimize layout of nodes {p, ..., p + n - 1}.
  void
  Subgraph::optimize(uint p)
  {
     // Initialize subgraph.
     const Float q = g->node[g->perm[p]].pos - g->node[g->perm[p]].hlen;
     min = WeightedSum(GECKO_FLOAT_MAX, 1);
     for (Subnode::Index k = 0; k < n; k++)
     {
        best[k] = perm[k] = k;
        Node::Index i = g->perm[p + k];
        // Copy i's outgoing arcs.  We distinguish between internal
        // and external arcs to nodes within and outside the subgraph,
        // respectively.
  #if GECKO_WITH_ADJLIST
        uint m = 0;
  #else
        adj[k] = 0;
  #endif
        std::vector<Arc::Index> external;
        for (Arc::Index a = g->node_begin(i); a < g->node_end(i); a++)
        {
           Node::Index j = g->adj[a];
           Subnode::Index l;
           for (l = 0; l < n && g->perm[p + l] != j; l++);
           if (l == n)
           {
              external.push_back(a);
           }
           else
           {
              // Copy internal arc to subgraph.
  #if GECKO_WITH_ADJLIST
              adj[k][m] = l;
              weight[k][m] = g->weight[a];
              m++;
  #else
              adj[k] += 1u << l;
              weight[k][l] = g->weight[a];
  #endif
           }
        }
  #if GECKO_WITH_ADJLIST
        adj[k][m] = k;
  #endif
        // Precompute external costs associated with all possible positions
        // of this node.  Since node lengths can be arbitrary, there are as
        // many as 2^(n-1) possible positions, each corresponding to an
        // (n-1)-bit string that specifies whether the remaining n-1 nodes
        // succeed this node or not.  Caching the
        //                n
        //   2^(n-1) n = sum k C(n, k) = A001787
        //               k=1
        // external costs is exponentially cheaper than recomputing the
        //      n-1         n
        //   n! sum 1/k! = sum k! C(n, k) = A007526
        //      k=0        k=1
        // costs associated with all permutations.
        node[k] = cache + (k << n);
        for (uint m = 0; m < (1u << n); m++)
           if (!(m & (1u << k)))
           {
              Subnode* s = cache + (k << n) + m;
              s->pos = q + g->node[i].hlen;
              for (Subnode::Index l = 0; l < n; l++)
                 if (l != k && !(m & (1u << l)))
                 {
                    s->pos += 2 * g->node[g->perm[p + l]].hlen;
                 }
              s->cost = g->cost(external, s->pos);
           }
           else
           {
              m += (1u << k) - 1;
           }
        node[k] += (1u << n) - (2u << k);
     }
  
     // Find optimal permutation of the n nodes.
     optimize(0, n);
  
     // Apply permutation to original graph.
     for (uint i = 0; i < n; i++)
     {
        g->swap(p + i, p + best[i]);
        for (uint j = i + 1; j < n; j++)
           if (best[j] == i)
           {
              best[j] = best[i];
           }
     }
  }
  
  // ----- graph.cpp --------------------------------------------------------------
  
  using namespace std;
  using namespace Gecko;
  
  // Constructor.
  void
  Graph::init(uint nodes)
  {
     node.push_back(Node(-1, 0, 1, Node::null));
     adj.push_back(Node::null);
     weight.push_back(0);
     bond.push_back(0);
     while (nodes--)
     {
        insert_node();
     }
  }
  
  // Insert node.
  Node::Index
  Graph::insert_node(Float length)
  {
     Node::Index p = Node::Index(node.size());
     perm.push_back(p);
     node.push_back(Node(-1, length));
     return p;
  }
  
  // Return nodes adjacent to i.
  std::vector<Node::Index>
  Graph::node_neighbors(Node::Index i) const
  {
     std::vector<Node::Index> neighbor;
     for (Arc::Index a = node_begin(i); a < node_end(i); a++)
     {
        neighbor.push_back(adj[a]);
     }
     return neighbor;
  }
  
  // Insert directed edge (i, j).
  Arc::Index
  Graph::insert_arc(Node::Index i, Node::Index j, Float w, Float b)
  {
     if (!i || !j || i == j || !(last_node <= i && i <= nodes()))
     {
        return Arc::null;
     }
     last_node = i;
     for (Node::Index k = i - 1; node[k].arc == Arc::null; k--)
     {
        node[k].arc = Arc::Index(adj.size());
     }
     adj.push_back(j);
     weight.push_back(w);
     bond.push_back(b);
     node[i].arc = Arc::Index(adj.size());
     return Arc::Index(adj.size() - 1);
  }
  
  // Remove arc a.
  bool
  Graph::remove_arc(Arc::Index a)
  {
     if (a == Arc::null)
     {
        return false;
     }
     Node::Index i = arc_source(a);
     adj.erase(adj.begin() + a);
     weight.erase(weight.begin() + a);
     bond.erase(bond.begin() + a);
     for (Node::Index k = i; k < node.size(); k++)
     {
        node[k].arc--;
     }
     return true;
  }
  
  // Remove directed edge (i, j).
  bool
  Graph::remove_arc(Node::Index i, Node::Index j)
  {
     return remove_arc(arc_index(i, j));
  }
  
  // Remove edge {i, j}.
  bool
  Graph::remove_edge(Node::Index i, Node::Index j)
  {
     bool success = remove_arc(i, j);
     if (success)
     {
        success = remove_arc(j, i);
     }
     return success;
  }
  
  // Index of arc (i, j) or null if not a valid arc.
  Arc::Index
  Graph::arc_index(Node::Index i, Node::Index j) const
  {
     for (Arc::Index a = node_begin(i); a < node_end(i); a++)
        if (adj[a] == j)
        {
           return a;
        }
     return Arc::null;
  }
  
  // Return source node i in arc a = (i, j).
  Node::Index
  Graph::arc_source(Arc::Index a) const
  {
     Node::Index j = adj[a];
     for (Arc::Index b = node_begin(j); b < node_end(j); b++)
     {
        Node::Index i = adj[b];
        if (node_begin(i) <= a && a < node_end(i))
        {
           return i;
        }
     }
     // should never get here
     throw std::runtime_error("internal data structure corrupted");
  }
  
  // Return reverse arc (j, i) of arc a = (i, j).
  Arc::Index
  Graph::reverse_arc(Arc::Index a) const
  {
     Node::Index j = adj[a];
     for (Arc::Index b = node_begin(j); b < node_end(j); b++)
     {
        Node::Index i = adj[b];
        if (node_begin(i) <= a && a < node_end(i))
        {
           return b;
        }
     }
     return Arc::null;
  }
  
  // Return first directed arc if one exists or null otherwise.
  Arc::Index
  Graph::directed() const
  {
     for (Node::Index i = 1; i < node.size(); i++)
        for (Arc::Index a = node_begin(i); a < node_end(i); a++)
        {
           Node::Index j = adj[a];
           if (!arc_index(j, i))
           {
              return a;
           }
        }
     return Arc::null;
  }
  
  // Add contribution of fine arc to coarse graph.
  void
  Graph::update(Node::Index i, Node::Index j, Float w, Float b)
  {
     Arc::Index a = arc_index(i, j);
     if (a == Arc::null)
     {
        insert_arc(i, j, w, b);
     }
     else
     {
        weight[a] += w;
        bond[a] += b;
     }
  }
  
  // Transfer contribution of fine arc a to coarse node p.
  void
  Graph::transfer(Graph* g, const vector<Float>& part, Node::Index p,
                  Arc::Index a, Float f) const
  {
     Float w = f * weight[a];
     Float m = f * bond[a];
     Node::Index j = arc_target(a);
     Node::Index q = node[j].parent;
     if (q == Node::null)
     {
        for (Arc::Index b = node_begin(j); b < node_end(j); b++)
           if (part[b] > 0)
           {
              q = node[adj[b]].parent;
              if (q != p)
              {
                 g->update(p, q, w * part[b], m * part[b]);
              }
           }
     }
     else
     {
        g->update(p, q, w, m);
     }
  }
  
  // Compute cost of a subset of arcs incident on node placed at pos.
  WeightedSum
  Graph::cost(const vector<Arc::Index>& subset, Float pos) const
  {
     WeightedSum c;
     for (Arc::ConstPtr ap = subset.begin(); ap != subset.end(); ap++)
     {
        Arc::Index a = *ap;
        Node::Index j = arc_target(a);
        Float l = fabs(node[j].pos - pos);
        Float w = weight[a];
        functional->accumulate(c, WeightedValue(l, w));
     }
     return c;
  }
  
  // Compute cost of graph layout.
  Float
  Graph::cost() const
  {
     if (edges())
     {
        WeightedSum c;
        Node::Index i = 1;
        for (Arc::Index a = 1; a < adj.size(); a++)
        {
           while (node_end(i) <= a)
           {
              i++;
           }
           Node::Index j = arc_target(a);
           Float l = length(i, j);
           Float w = weight[a];
           functional->accumulate(c, WeightedValue(l, w));
        }
        return functional->mean(c);
     }
     else
     {
        return Float(0);
     }
  }
  
  // Swap the two nodes in positions k and l, k <= l.
  void
  Graph::swap(uint k, uint l)
  {
     Node::Index i = perm[k];
     perm[k] = perm[l];
     perm[l] = i;
     Float p = node[i].pos - node[i].hlen;
     do
     {
        i = perm[k];
        p += node[i].hlen;
        node[i].pos = p;
        p += node[i].hlen;
     }
     while (k++ != l);
  }
  
  // Optimize continuous position of a single node.
  Float
  Graph::optimal(Node::Index i) const
  {
     vector<WeightedValue> v;
     for (Arc::Index a = node_begin(i); a < node_end(i); a++)
     {
        Node::Index j = adj[a];
        if (placed(j))
        {
           v.push_back(WeightedValue(node[j].pos, weight[a]));
        }
     }
     return v.empty() ? -1 : functional->optimum(v);
  }
  
  // Compute coarse graph with roughly half the number of nodes.
  Graph*
  Graph::coarsen()
  {
     progress->beginphase(this, string("coarse"));
     Graph* g = new Graph(0, level - 1);
     g->functional = functional;
     g->progress = progress;
  
     // Compute importance of nodes in fine graph.
     DynamicHeap<Node::Index, Float> heap;
     for (Node::Index i = 1; i < node.size(); i++)
     {
        node[i].parent = Node::null;
        Float w = 0;
        for (Arc::Index a = node_begin(i); a < node_end(i); a++)
        {
           w += bond[a];
        }
        heap.insert(i, w);
     }
  
     // Select set of important nodes from fine graph that will remain in
     // coarse graph.
     vector<Node::Index> child(1, Node::null);
     while (!heap.empty())
     {
        Node::Index i;
        Float w = 0;
        heap.extract(i, w);
        if (w < 0)
        {
           break;
        }
        child.push_back(i);
        node[i].parent = g->insert_node(2 * node[i].hlen);
  
        // Reduce importance of neighbors.
        for (Arc::Index a = node_begin(i); a < node_end(i); a++)
        {
           Node::Index j = adj[a];
           if (heap.find(j, w))
           {
              heap.update(j, w - 2 * bond[a]);
           }
        }
     }
  
     // Assign parts of remaining nodes to aggregates.
     vector<Float> part = bond;
     for (Node::Index i = 1; i < node.size(); i++)
        if (!persistent(i))
        {
           // Find all connections to coarse nodes.
           Float w = 0;
           Float max = 0;
           for (Arc::Index a = node_begin(i); a < node_end(i); a++)
           {
              Node::Index j = adj[a];
              if (persistent(j))
              {
                 w += part[a];
                 if (max < part[a])
                 {
                    max = part[a];
                 }
              }
              else
              {
                 part[a] = -1;
              }
           }
           max /= GECKO_PART_FRAC;
  
           // Weed out insignificant connections.
           for (Arc::Index a = node_begin(i); a < node_end(i); a++)
              if (0 < part[a] && part[a] < max)
              {
                 w -= part[a];
                 part[a] = -1;
              }
  
           // Compute node fractions (interpolation matrix) and assign
           // partial nodes to aggregates.
           for (Arc::Index a = node_begin(i); a < node_end(i); a++)
              if (part[a] > 0)
              {
                 part[a] /= w;
                 Node::Index p = node[adj[a]].parent;
                 g->node[p].hlen += part[a] * node[i].hlen;
              }
        }
  
     // Transfer arcs to coarse graph.
     for (Node::Index p = 1; p < g->node.size(); p++)
     {
        Node::Index i = child[p];
        for (Arc::Index a = node_begin(i); a < node_end(i); a++)
        {
           transfer(g, part, p, a);
           Node::Index j = adj[a];
           if (!persistent(j))
           {
              Arc::Index b = arc_index(j, i);
              if (part[b] > 0)
                 for (Arc::Index c = node_begin(j); c < node_end(j); c++)
                 {
                    Node::Index k = adj[c];
                    if (k != i)
                    {
                       transfer(g, part, p, c, part[b]);
                    }
                 }
           }
        }
     }
  
  #if DEBUG
     if (g->directed())
     {
        throw runtime_error("directed edge found");
     }
  #endif
  
     // Free memory.
     vector<Float> t = bond;
     bond.swap(t);
  
     progress->endphase(this, false);
  
     return g;
  }
  
  // Order nodes according to coarsened graph layout.
  void
  Graph::refine(const Graph* graph)
  {
     progress->beginphase(this, string("refine"));
  
     // Place persistent nodes.
     DynamicHeap<Node::Index, Float> heap;
     for (Node::Index i = 1; i < node.size(); i++)
        if (persistent(i))
        {
           Node::Index p = node[i].parent;
           node[i].pos = graph->node[p].pos;
        }
        else
        {
           node[i].pos = -1;
           Float w = 0;
           for (Arc::Index a = node_begin(i); a < node_end(i); a++)
           {
              Node::Index j = adj[a];
              if (persistent(j))
              {
                 w += weight[a];
              }
           }
           heap.insert(i, w);
        }
  
     // Place remaining nodes in order of decreasing connectivity with
     // already placed nodes.
     while (!heap.empty())
     {
        Node::Index i = 0;
        heap.extract(i);
        node[i].pos = optimal(i);
        for (Arc::Index a = node_begin(i); a < node_end(i); a++)
        {
           Node::Index j = adj[a];
           Float w;
           if (heap.find(j, w))
           {
              heap.update(j, w + weight[a]);
           }
        }
     }
  
     place(true);
     progress->endphase(this, true);
  }
  
  // Perform m sweeps of compatible or Gauss-Seidel relaxation.
  void
  Graph::relax(bool compatible, uint m)
  {
     progress->beginphase(this, compatible ? string("crelax") : string("frelax"));
     while (m--)
        for (uint k = 0; k < perm.size() && !progress->quit(); k++)
        {
           Node::Index i = perm[k];
           if (!compatible || !persistent(i))
           {
              node[i].pos = optimal(i);
           }
        }
     place(true);
     progress->endphase(this, true);
  }
  
  // Optimize successive n-node subgraphs.
  void
  Graph::optimize(uint n)
  {
     if (n > perm.size())
     {
        n = uint(perm.size());
     }
     ostringstream count;
     count << setw(2) << n;
     progress->beginphase(this, string("perm") + count.str());
     Subgraph* subgraph = new Subgraph(this, n);
     for (uint k = 0; k <= perm.size() - n && !progress->quit(); k++)
     {
        subgraph->optimize(k);
     }
     delete subgraph;
     progress->endphase(this, true);
  }
  
  // Place all nodes according to their positions.
  void
  Graph::place(bool sort)
  {
     place(sort, 0, uint(perm.size()));
  }
  
  // Place nodes {k, ..., k + n - 1} according to their positions.
  void
  Graph::place(bool sort, uint k, uint n)
  {
     // Place nodes.
     if (sort)
     {
        stable_sort(perm.begin() + k, perm.begin() + k + n,
                    Node::Comparator(node.begin()));
     }
  
     // Assign node positions according to permutation.
     for (Float p = k ? node[perm[k - 1]].pos + node[perm[k - 1]].hlen : 0; n--;
          k++)
     {
        Node::Index i = perm[k];
        p += node[i].hlen;
        node[i].pos = p;
        p += node[i].hlen;
     }
  }
  
  // Perform one V-cycle.
  void
  Graph::vcycle(uint n, uint work)
  {
     if (n < nodes() && nodes() < edges() && level && !progress->quit())
     {
        Graph* graph = coarsen();
        graph->vcycle(n, work + edges());
        refine(graph);
        delete graph;
     }
     else
     {
        place();
     }
     if (edges())
     {
        relax(true, GECKO_CR_SWEEPS);
        relax(false, GECKO_GS_SWEEPS);
        for (uint w = edges(); w * (n + 1) < work; w *= ++n);
        n = std::min(n, uint(GECKO_WINDOW_MAX));
        if (n)
        {
           optimize(n);
        }
     }
  }
  
  // Custom random-number generator for reproducibility.
  // LCG from doi:10.1090/S0025-5718-99-00996-5.
  uint
  Graph::random(uint seed)
  {
     static uint state = 1;
     state = (seed ? seed : 0x1ed0675 * state + 0xa14f);
     return state;
  }
  
  // Generate a random permutation of the nodes.
  void
  Graph::shuffle(uint seed)
  {
     random(seed);
     for (uint k = 0; k < perm.size(); k++)
     {
        uint r = random() >> 8;
        uint l = k + r % (uint(perm.size()) - k);
        std::swap(perm[k], perm[l]);
     }
     place();
  }
  
  // Recompute bonds for k'th V-cycle.
  void
  Graph::reweight(uint k)
  {
     bond.resize(weight.size());
     for (Arc::Index a = 1; a < adj.size(); a++)
     {
        bond[a] = functional->bond(weight[a], length(a), k);
     }
  }
  
  // Linearly order graph.
  void
  Graph::order(Functional* functional, uint iterations, uint window, uint period,
               uint seed, Progress* progress)
  {
     // Initialize graph.
     this->functional = functional;
     progress = this->progress = progress ? progress : new Progress;
     for (level = 0; (1u << level) < nodes(); level++);
     place();
     Float mincost = cost();
     vector<Node::Index> minperm = perm;
     if (seed)
     {
        shuffle(seed);
     }
  
     progress->beginorder(this, mincost);
     if (edges())
     {
        // Perform specified number of V-cycles.
        for (uint k = 1; k <= iterations && !progress->quit(); k++)
        {
           progress->beginiter(this, k, iterations, window);
           reweight(k);
           vcycle(window);
           Float c = cost();
           if (c < mincost)
           {
              mincost = c;
              minperm = perm;
           }
           progress->enditer(this, mincost, c);
           if (period && !(k % period))
           {
              window++;
           }
        }
        perm = minperm;
        place();
     }
     progress->endorder(this, mincost);
  
     if (!progress)
     {
        delete this->progress;
        this->progress = 0;
     }
  }
  