Branch data     Line data    Source code

       1                 :            : #include <mutable/util/AdjacencyMatrix.hpp>
       2                 :            : 
       3                 :            : #include <queue>
       4                 :            : 
       5                 :            : 
       6                 :            : using namespace m;
       7                 :            : 
       8                 :            : 
       9                 :            : /*======================================================================================================================
      10                 :            :  * AdjacencyMatrix
      11                 :            :  *====================================================================================================================*/
      12                 :            : 
      13                 :          3 : AdjacencyMatrix AdjacencyMatrix::transitive_closure_directed() const
      14                 :            : {
      15                 :          3 :     AdjacencyMatrix closure(*this); // copy
      16                 :            : 
      17                 :            :     bool changed;
      18                 :          3 :     do {
      19                 :          4 :         changed = false;
      20         [ +  + ]:         22 :         for (std::size_t i = 0; i != num_vertices_; ++i) {
      21         [ +  + ]:        102 :             for (std::size_t j = 0; j != num_vertices_; ++j) {
      22         [ +  + ]:         84 :                 if (not closure(i, j)) {
      23                 :         67 :                     const SmallBitset row = closure.m_[i];
      24                 :         67 :                     const SmallBitset col_mask = SmallBitset::Singleton(j);
      25         [ +  + ]:        117 :                     for (std::size_t k : row) {
      26                 :         54 :                         M_insist(row[k]);
      27         [ +  + ]:         54 :                         if (closure.m_[k] & col_mask) {
      28                 :          4 :                             closure(i, j) = true;
      29                 :          4 :                             changed = true;
      30                 :          4 :                             break;
      31                 :            :                         }
      32                 :            :                     }
      33                 :         67 :                 }
      34                 :         84 :             }
      35                 :         18 :         }
      36         [ +  + ]:          4 :     } while (changed);
      37                 :            : 
      38                 :          3 :     return closure;
      39                 :            : }
      40                 :            : 
      41                 :          4 : AdjacencyMatrix AdjacencyMatrix::transitive_closure_undirected() const
      42                 :            : {
      43                 :          4 :     AdjacencyMatrix closure(*this); // copy
      44                 :            : 
      45                 :            :     bool changed;
      46                 :          4 :     do {
      47                 :          6 :         changed = false;
      48         [ +  + ]:         30 :         for (std::size_t i = 0; i != num_vertices_; ++i) {
      49         [ +  + ]:         87 :             for (std::size_t j = i; j != num_vertices_; ++j) { // exploit symmetry
      50                 :         63 :                 M_insist(closure(i, j) == closure(j, i), "not symmetric");
      51                 :         63 :                 const bool before = closure(i, j);
      52                 :         63 :                 const SmallBitset row = closure.m_[i];
      53                 :         63 :                 const SmallBitset col = closure.m_[j]; // exploit symmetry
      54                 :         63 :                 const bool dot = not (row & col).empty(); // compute connected-ness as "dot product"
      55                 :         63 :                 closure(i, j) = closure(j, i) = dot;
      56         [ +  + ]:         63 :                 changed = changed or before != dot;
      57                 :         63 :             }
      58                 :         24 :         }
      59         [ +  + ]:          6 :     } while (changed);
      60                 :            : 
      61                 :          4 :     return closure;
      62                 :            : }
      63                 :            : 
      64                 :          8 : AdjacencyMatrix AdjacencyMatrix::minimum_spanning_forest(std::function<double(std::size_t, std::size_t)> weight) const
      65                 :            : {
      66                 :          8 :     AdjacencyMatrix MSF(num_vertices_);
      67                 :            : 
      68         [ +  + ]:          8 :     if (num_vertices_ == 0)
      69                 :          1 :         return MSF;
      70                 :            : 
      71                 :            :     struct weighted_edge
      72                 :            :     {
      73                 :            :         std::size_t source, sink;
      74                 :          1 :         double weight;
      75                 :            : 
      76                 :            :         weighted_edge() = default;
      77                 :         22 :         weighted_edge(std::size_t source, std::size_t sink, double weight)
      78                 :         22 :             : source(source), sink(sink), weight(weight)
      79                 :         22 :         { }
      80                 :            : 
      81                 :            :         bool operator<(const weighted_edge &other) const { return this->weight < other.weight; }
      82                 :         23 :         bool operator>(const weighted_edge &other) const { return this->weight > other.weight; }
      83                 :            :     };
      84                 :          7 :     std::priority_queue<weighted_edge, std::vector<weighted_edge>, std::greater<weighted_edge>> Q;
      85                 :            : 
      86                 :            :     /* Run Prim's algorithm for each remaining vertex to compute a MSF. */
      87         [ +  - ]:          7 :     SmallBitset vertices_remaining = SmallBitset::All(num_vertices_);
      88                 :            : 
      89   [ +  -  +  + ]:         21 :     while (vertices_remaining) {
      90   [ +  -  +  - ]:         14 :         SmallBitset next_vertex = vertices_remaining.begin().as_set();
      91         [ +  - ]:         14 :         vertices_remaining = vertices_remaining - next_vertex;
      92                 :            : 
      93                 :            :         /* Prim's algorithm for finding a MST. */
      94   [ +  -  +  - ]:         29 :         while (next_vertex) {
      95                 :            :             /* Explore edges of `next_vertex`. */
      96   [ +  -  +  - ]:         29 :             M_insist(next_vertex.is_singleton());
      97   [ +  -  +  - ]:         29 :             const SmallBitset N = neighbors(next_vertex) & vertices_remaining;
      98   [ +  -  +  - ]:         29 :             const std::size_t u = *next_vertex.begin();
      99   [ +  -  +  -  :         51 :             for (std::size_t v : N)
          +  -  +  +  +  
                -  +  - ]
     100   [ +  -  +  - ]:         22 :                 Q.emplace(u, v, weight(u, v));
     101                 :            : 
     102                 :            :             /* Search for the cheapest edge not within the MSF. */
     103                 :            :             weighted_edge E;
     104                 :         29 :             do {
     105   [ +  -  +  + ]:         36 :                 if (Q.empty())
     106                 :         14 :                     goto MST_done;
     107         [ +  - ]:         22 :                 E = Q.top();
     108         [ +  - ]:         22 :                 Q.pop();
     109   [ +  -  +  -  :         22 :             } while (not vertices_remaining[E.sink]);
                   +  + ]
     110                 :            : 
     111                 :            :             /* Add edge to MSF. */
     112   [ +  -  +  -  :         15 :             M_insist(vertices_remaining[E.sink], "sink must not yet be in the MSF");
                   +  - ]
     113   [ +  -  +  - ]:         15 :             vertices_remaining[E.sink] = false;
     114   [ +  -  +  -  :         15 :             MSF(E.source, E.sink) = MSF(E.sink, E.source) = true; // add undirected edge to MSF
             +  -  +  - ]
     115         [ +  - ]:         15 :             next_vertex = SmallBitset::Singleton(E.sink);
     116                 :            :         }
     117                 :            : MST_done: /* MST is complete */;
     118                 :            :     }
     119                 :            :     return MSF;
     120                 :          8 : }
     121                 :            : 
     122                 :         11 : AdjacencyMatrix AdjacencyMatrix::tree_directed_away_from(SmallBitset root)
     123                 :            : {
     124                 :         11 :     M_insist(root.is_singleton());
     125                 :         11 :     AdjacencyMatrix directed_tree(num_vertices_);
     126                 :         11 :     SmallBitset V = root;
     127                 :         11 :     SmallBitset X;
     128         [ +  + ]:         41 :     while (not V.empty()) {
     129                 :         30 :         X |= V;
     130         [ +  + ]:         67 :         for (auto v : V)
     131                 :         37 :             directed_tree.m_[v] = this->m_[v] - X;
     132                 :         30 :         V = directed_tree.neighbors(V);
     133                 :            :     }
     134                 :            : 
     135                 :         11 :     return directed_tree;
     136                 :            : }
     137                 :            : 
     138                 :            : M_LCOV_EXCL_START
     139                 :            : void AdjacencyMatrix::dump(std::ostream &out) const { out << *this << std::endl; }
     140                 :            : void AdjacencyMatrix::dump() const { dump(std::cerr); }
     141                 :            : M_LCOV_EXCL_STOP