|
| 1 | +/******************************************************************************* |
| 2 | +
|
| 3 | + This file is part of Set--Based Graph Library. |
| 4 | +
|
| 5 | + SBG Library is free software: you can redistribute it and/or modify |
| 6 | + it under the terms of the GNU General Public License as published by |
| 7 | + the Free Software Foundation, either version 3 of the License, or |
| 8 | + (at your option) any later version. |
| 9 | +
|
| 10 | + SBG Library is distributed in the hope that it will be useful, |
| 11 | + but WITHOUT ANY WARRANTY; without even the implied warranty of |
| 12 | + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
| 13 | + GNU General Public License for more details. |
| 14 | +
|
| 15 | + You should have received a copy of the GNU General Public License |
| 16 | + along with SBG Library. If not, see <http://www.gnu.org/licenses/>. |
| 17 | +
|
| 18 | + ******************************************************************************/ |
| 19 | + |
| 20 | +#include <chrono> |
| 21 | + |
| 22 | +#include "algorithms/tearing/tearing.hpp" |
| 23 | +#include "util/logger.hpp" |
| 24 | + |
| 25 | +namespace SBG { |
| 26 | + |
| 27 | +namespace LIB { |
| 28 | + |
| 29 | +//////////////////////////////////////////////////////////////////////////////// |
| 30 | +// Tearing --------------------------------------------------------------------- |
| 31 | +//////////////////////////////////////////////////////////////////////////////// |
| 32 | + |
| 33 | +bool sccNotEqId(const Map &sbgmap) { return !(sbgmap.isId()); } |
| 34 | + |
| 35 | +Tearing::Tearing(const DSBG &dsbg, bool debug) |
| 36 | + : fact_(dsbg.fact()), |
| 37 | + dsbg_(dsbg), |
| 38 | + V_(fact_.createSet()), |
| 39 | + Vmap_(fact_.createPWMap()), |
| 40 | + Emap_(fact_.createPWMap()), |
| 41 | + subEmap_(fact_.createPWMap()), |
| 42 | + E_(fact_.createSet()), |
| 43 | + Ediff_(fact_.createSet()), |
| 44 | + mapB_(fact_.createPWMap()), |
| 45 | + mapD_(fact_.createPWMap()), |
| 46 | + rmap_(fact_.createPWMap()), |
| 47 | + debug_(debug) |
| 48 | +{ |
| 49 | + DSBG dg = dsbg; |
| 50 | + dsbg_ = dg; |
| 51 | + |
| 52 | + V_ = dg.V(); |
| 53 | + Vmap_ = dg.Vmap(); |
| 54 | + |
| 55 | + E_ = dg.E(); |
| 56 | + Emap_ = dg.Emap(); |
| 57 | + subEmap_ = dg.subEmap(); |
| 58 | + |
| 59 | + mapB_ = dg.mapB(); |
| 60 | + mapD_ = dg.mapD(); |
| 61 | + |
| 62 | + rmap_ = fact_.createPWMap(V_); |
| 63 | +} |
| 64 | + |
| 65 | +member_imp(Tearing, DSBG, dsbg); |
| 66 | +member_imp(Tearing, Set, V); |
| 67 | +member_imp(Tearing, PWMap, Vmap); |
| 68 | +member_imp(Tearing, Set, E); |
| 69 | +member_imp(Tearing, PWMap, Emap); |
| 70 | +member_imp(Tearing, PWMap, subEmap); |
| 71 | + |
| 72 | +member_imp(Tearing, PWMap, mapB); |
| 73 | +member_imp(Tearing, PWMap, mapD); |
| 74 | + |
| 75 | +member_imp(Tearing, Set, Ediff); |
| 76 | + |
| 77 | +member_imp(Tearing, PWMap, rmap); |
| 78 | + |
| 79 | +member_imp(Tearing, bool, debug); |
| 80 | + |
| 81 | +PWMap Tearing::sccMinReach(const DSBG &dg) const |
| 82 | +{ |
| 83 | + if (debug()) |
| 84 | + Util::SBG_LOG << "Min reach graph:\n" << dg << "\n\n"; |
| 85 | + |
| 86 | + Set V = dg.V(), E = dg.E(); |
| 87 | + PWMap mapB = dg.mapB(), mapD = dg.mapD(), subEmap = dg.subEmap(); |
| 88 | + if (!V.isEmpty()) { |
| 89 | + unsigned int copies = V.arity(); |
| 90 | + PWMap rmap = fact_.createPWMap(V), old_rmap = fact_.createPWMap(); |
| 91 | + |
| 92 | + if (E.isEmpty()) |
| 93 | + return rmap; |
| 94 | + |
| 95 | + Set Vc = fact_.createSet(); |
| 96 | + do { |
| 97 | + old_rmap = rmap; |
| 98 | + |
| 99 | + PWMap ermapD = rmap.composition(mapD); |
| 100 | + |
| 101 | + PWMap new_rmap = mapB.minAdjMap(ermapD); |
| 102 | + rmap = rmap.minMap(new_rmap).combine(rmap); |
| 103 | + if (debug()) |
| 104 | + Util::SBG_LOG << "rmap before rec: " << rmap << "\n\n"; |
| 105 | + |
| 106 | + Set positive = fact_.createSet(SetPiece(copies, Interval(1, 1, Inf))); |
| 107 | + PWMap rec_rmap = fact_.createPWMap(); |
| 108 | + Vc = V.difference(old_rmap.equalImage(rmap)); |
| 109 | + if (!Vc.isEmpty()) { |
| 110 | + for (const Map &subv : dg.Vmap()) { |
| 111 | + Set vs = subv.dom(); |
| 112 | + if (!vs.intersection(Vc).isEmpty()) { |
| 113 | + // If the mrv is in the same SV, the algorithm would detect a false |
| 114 | + // recursion, i.e. if we have a cycle 1 -> 2 -> ... -> 10 -> 1, |
| 115 | + // where SV = [1:10], then it detects a recursion when mrv(10) |
| 116 | + // becomes "1" (false recursion). So we take self mrvs out. |
| 117 | + auto other_rep = rmap.filterMap([](const Map &sbgmap) { |
| 118 | + return sccNotEqId(sbgmap); |
| 119 | + }).dom(); |
| 120 | + // Vertices in the set-vertex that share its rep with other vertex |
| 121 | + // in the set-vertex |
| 122 | + Set VR = rmap.restrict(vs.intersection(other_rep)).sharedImage(); |
| 123 | + // There is a recursive vertex that changed its rep in the last step |
| 124 | + // (to avoid computing again an already found recursion) |
| 125 | + if (!VR.intersection(Vc).isEmpty()) { |
| 126 | + // Vertices that reach the shared representative |
| 127 | + Set repV = rmap.preImage(rmap.image(VR)); |
| 128 | + Set end = VR.difference(Vc); |
| 129 | + |
| 130 | + // Edges with both endings in VR (path to a minimum rep) |
| 131 | + Set ERB = mapB.preImage(repV), ERD = mapD.preImage(repV); |
| 132 | + Set ER = ERB.intersection(ERD); |
| 133 | + if (!end.isEmpty() && !ER.isEmpty()) { |
| 134 | + // Distance map |
| 135 | + PWMap dmap = fact_.createPWMap(); |
| 136 | + Set ith = end; |
| 137 | + NAT dist = 0; |
| 138 | + // Calculate distance for vertices in same_rep that reach reps |
| 139 | + for (; dmap.dom().intersection(Vc.intersection(VR)).isEmpty();) { |
| 140 | + Set dom = ith.difference(dmap.dom()); |
| 141 | + Exp exp(MD_NAT(copies, dist)); |
| 142 | + dmap.emplaceBack(fact_.createMap(dom, exp)); |
| 143 | + // Update ith to vertices that have outgoing edges entering ith |
| 144 | + ith = mapB.image(mapD.preImage(ith)); |
| 145 | + ++dist; |
| 146 | + } |
| 147 | + PWMap dmapB = dmap.composition(mapB), dmapD = dmap.composition(mapD); |
| 148 | + // Get edges where the end is closer to the rep than the beginning |
| 149 | + Set not_cycle_edges = (dmapB - dmapD).preImage(positive); |
| 150 | + ER = ER.intersection(not_cycle_edges); |
| 151 | + |
| 152 | + // Extend to subset-edge |
| 153 | + Set ER_plus = subEmap.preImage(subEmap.image(ER)); |
| 154 | + // Calculate a succesor |
| 155 | + PWMap auxB = mapB.restrict(ER_plus), auxD = mapD.restrict(ER_plus); |
| 156 | + PWMap smap_plus = rmap.restrict(VR); |
| 157 | + smap_plus = smap_plus.combine(auxB.minAdjMap(auxD)); |
| 158 | + |
| 159 | + // Update rmap for recursion, and leave the rest unchanged |
| 160 | + rec_rmap = smap_plus.combine(rec_rmap).compact(); |
| 161 | + |
| 162 | + if (debug()) { |
| 163 | + Util::SBG_LOG << "VR: " << VR << "\n"; |
| 164 | + Util::SBG_LOG << "repV: " << repV << "\n"; |
| 165 | + Util::SBG_LOG << "ER: " << ER << "\n"; |
| 166 | + Util::SBG_LOG << "dmap: " << dmap << "\n"; |
| 167 | + Util::SBG_LOG << "not_cycle_edges: " << not_cycle_edges << "\n"; |
| 168 | + Util::SBG_LOG << "ER_plus: " << ER_plus << "\n"; |
| 169 | + Util::SBG_LOG << "smap_plus: " << smap_plus << "\n"; |
| 170 | + Util::SBG_LOG << "rec_rmap: " << rec_rmap << "\n\n"; |
| 171 | + } |
| 172 | + } |
| 173 | + } |
| 174 | + } |
| 175 | + } |
| 176 | + PWMap rmap_plus = rec_rmap.combine(rmap); |
| 177 | + rmap_plus = rmap_plus.mapInf(); |
| 178 | + rmap = rmap.minMap(rmap_plus).compact(); |
| 179 | + |
| 180 | + if (debug()) |
| 181 | + Util::SBG_LOG << "rmap after rec: " << rmap << "\n\n"; |
| 182 | + } |
| 183 | + } while (!Vc.isEmpty()); |
| 184 | + |
| 185 | + return rmap; |
| 186 | + } |
| 187 | + |
| 188 | + return fact_.createPWMap(); |
| 189 | +} |
| 190 | + |
| 191 | +PWMap Tearing::sccStep() |
| 192 | +{ |
| 193 | + PWMap id_V = fact_.createPWMap(V()); |
| 194 | + |
| 195 | + DSBG aux_dsbg( |
| 196 | + fact_, V(), Vmap() |
| 197 | + , mapB().restrict(E()), mapD().restrict(E()), Emap().restrict(E()).compact() |
| 198 | + , subEmap().restrict(E()) |
| 199 | + ); |
| 200 | + PWMap new_rmap = sccMinReach(aux_dsbg); |
| 201 | + if (debug()) |
| 202 | + Util::SBG_LOG << "SCC new_rmap: " << new_rmap << "\n"; |
| 203 | + |
| 204 | + PWMap rmap_B = new_rmap.composition(mapB()); |
| 205 | + PWMap rmap_D = new_rmap.composition(mapD()); |
| 206 | + Set Esame = rmap_B.equalImage(rmap_D); // Edges in the same SCC |
| 207 | + |
| 208 | + // Leave edges in the same SCC |
| 209 | + Ediff_ = E().difference(Esame); |
| 210 | + E_ = Esame; |
| 211 | + if (debug()) |
| 212 | + Util::SBG_LOG << "SCC erased edges: " << Ediff() << "\n\n"; |
| 213 | + |
| 214 | + // Swap directions |
| 215 | + PWMap aux_B = mapB(); |
| 216 | + mapB_ = mapD(); |
| 217 | + mapD_ = aux_B; |
| 218 | + |
| 219 | + return new_rmap; |
| 220 | +} |
| 221 | + |
| 222 | +void Tearing::restoreSBG() |
| 223 | +{ |
| 224 | + V_ = dsbg_.V(); |
| 225 | + Vmap_ = dsbg_.Vmap(); |
| 226 | + |
| 227 | + E_ = dsbg_.E(); |
| 228 | + Emap_ = dsbg_.Emap(); |
| 229 | + subEmap_ = dsbg_.subEmap(); |
| 230 | + |
| 231 | + mapB_ = dsbg_.mapB(); |
| 232 | + mapD_ = dsbg_.mapD(); |
| 233 | +} |
| 234 | + |
| 235 | +PWMap Tearing::calculate() |
| 236 | +{ |
| 237 | + if (debug()) Util::SBG_LOG << "Tearing dsbg: \n" << dsbg() << "\n\n"; |
| 238 | + |
| 239 | + auto begin = std::chrono::high_resolution_clock::now(); |
| 240 | + PWMap rmap = sccStep(); |
| 241 | + Set vrem = fact_.createSet(); |
| 242 | + Set e_scc = fact_.createSet(); |
| 243 | + do { |
| 244 | + // SCC |
| 245 | + do { |
| 246 | + rmap = sccStep(); |
| 247 | + } while (Ediff() != fact_.createSet()); |
| 248 | + // |
| 249 | + e_scc = E(); |
| 250 | + if (E() != fact_.createSet()) { |
| 251 | + Set erem_b = mapB().preImage(rmap.image()); |
| 252 | + Set erem_d = mapD().preImage(rmap.image()); |
| 253 | + vrem = vrem.cup(mapB().image(erem_b)).cup(mapD().image(erem_d)); |
| 254 | + E_ = E().difference(erem_b.cup(erem_d)); |
| 255 | + } |
| 256 | + rmap_ = rmap.compact(); |
| 257 | + } while (E() != fact_.createSet()); |
| 258 | + restoreSBG(); |
| 259 | + Set e_notscc = E().difference(e_scc); |
| 260 | + Set e_to_R = mapD().preImage(rmap.image()); |
| 261 | + PWMap mapD_notscc = mapD().restrict(e_to_R); |
| 262 | + |
| 263 | + Set erem = e_scc.intersection(mapB().preImage(vrem)); |
| 264 | + Set mapD_aux = mapD().image(erem); |
| 265 | + |
| 266 | + auto end = std::chrono::high_resolution_clock::now(); |
| 267 | + |
| 268 | + auto total = std::chrono::duration_cast<std::chrono::microseconds>(end - begin); |
| 269 | + Util::SBG_LOG << "Total Tearing exec time: " << total.count() << " [μs]\n\n"; |
| 270 | + std::cout << vrem << std::endl; |
| 271 | + std::cout << dsbg().E() << std::endl; |
| 272 | + |
| 273 | + if (debug()) { |
| 274 | + Util::SBG_LOG << "Tearing result: " << rmap.compact() << "\n\n"; |
| 275 | + E().print(std::cout); |
| 276 | + } |
| 277 | + |
| 278 | + return rmap.compact(); |
| 279 | +} |
| 280 | + |
| 281 | +const PWMapAF &Tearing::fact() const { return fact_; } |
| 282 | + |
| 283 | +} // namespace LIB |
| 284 | + |
| 285 | +} // namespace SBG |
0 commit comments