Finished topological sort

sbg/sbg_algorithms.cpp

@@ -849,71 +849,93 @@ std::ostream &operator<<(std::ostream &out, const VertexOrder<Set> &vo)
 
 template<typename Set>
 SBGTopSort<Set>::SBGTopSort()
-  : dsbg_(), mapB_(), mapD_(), disordered_(), deletedE_(), debug_(false) {}
+  : dsbg_(), mapB_(), mapD_(), disordered_(), rec_map_(), debug_(false) {}
 template<typename Set>
 SBGTopSort<Set>::SBGTopSort(DSBGraph<Set> dsbg, bool debug)
   : dsbg_(dsbg), mapB_(dsbg.mapB()), mapD_(dsbg.mapD()), disordered_(dsbg.V())
-    , deletedE_(), debug_(debug) {}
+    , rec_map_(), debug_(debug) {}
 
 member_imp_temp(template<typename Set>, SBGTopSort<Set>, DSBGraph<Set>, dsbg);
 member_imp_temp(template<typename Set>, SBGTopSort<Set>, PWMap<Set>, mapB);
 member_imp_temp(template<typename Set>, SBGTopSort<Set>, PWMap<Set>, mapD);
 member_imp_temp(template<typename Set>, SBGTopSort<Set>, Set, disordered);
-member_imp_temp(template<typename Set>, SBGTopSort<Set>, Set, deletedE);
+member_imp_temp(template<typename Set>, SBGTopSort<Set>, PWMap<Set>, rec_map);
 member_imp_temp(template<typename Set>, SBGTopSort<Set>, bool, debug);
 
 template<typename Set>
-Set SBGTopSort<Set>::topSortStep()
+void SBGTopSort<Set>::findRecursions()
 {
+  PW res;
+
   DSBGraph<Set> aux_dsbg = dsbg();
   PW subE_map = aux_dsbg.subE_map();
 
-  // Vertices without outgoing edges
-  Set nth = disordered().difference(mapB().image());
+  Set unvisited = subE_map.dom(); 
+  Util::NAT dims = subE_map.nmbrDims();
+  unsigned int k = 1;
+  for (const SBGMap<Set> &map : subE_map) {
+    Set first = map.dom().intersection(unvisited), last = first;
 
-  // Ingoing edges to nth
-  Set ingoing = mapD().preImage(nth), ingoing_plus = ingoing;
-  // Subset edges id's with an ingoing edge to nth
-  Set ith_ids = subE_map.image(ingoing);
-  Set repeated = ith_ids.intersection(subE_map.image(deletedE()));
-  // Already visited subset edge
-  if (!repeated.isEmpty()) {
-    Set start_rec = ingoing.intersection(subE_map.preImage(repeated));
-    Set end_rec = deletedE().intersection(subE_map.preImage(repeated));
-
-    for (const SBGMap<Set> &map : subE_map) {
-      Set sub_first = map.dom().intersection(start_rec);
-      Set sub_last = map.dom().intersection(end_rec);
-
-      Set v_start = aux_dsbg.mapD().image(sub_first); // Ending vertex of start
+    if (!first.isEmpty()) {
+      Set v_start = aux_dsbg.mapD().image(first); // Ending vertex of start
       Set start_ith = aux_dsbg.mapB().preImage(v_start); // Following edge from the start
-      Set v_end = aux_dsbg.mapB().image(sub_last); // Starting vertex of end
+      Set v_end = aux_dsbg.mapB().image(last); // Starting vertex of end
       Set end_ith = aux_dsbg.mapD().preImage(v_end); // Following edge from the end
-      Set start_path = sub_first, end_path = sub_last; 
-      for (unsigned int j = 0; j < 3; ++j) {
-        start_path = start_path.cup(start_ith);
-        end_path = end_path.cup(end_ith);
-        if (!start_ith.intersection(sub_last).isEmpty())
-          break;
-
-        v_start = aux_dsbg.mapD().image(start_ith); // Ending vertex of start_ith
-        start_ith = aux_dsbg.mapB().preImage(v_start); // Following edge from the start_ith
-        v_end = aux_dsbg.mapB().image(end_ith); // Starting vertex of end_ith
-        end_ith = aux_dsbg.mapD().preImage(v_end); // Following edge from the end_ith
+      Set start_path = first, end_path = last, rec_path; 
+
+      if (!start_ith.intersection(first).isEmpty()) {
+        rec_path = start_path.intersection(end_path);
+        rec_path = subE_map.preImage(subE_map.image(rec_path));
+
+        unvisited = unvisited.difference(rec_path);
+        Util::MD_NAT id(dims, k);
+        rec_map_ref().emplaceBack(SBGMap<Set>(rec_path, Exp(id)));
+        ++k;
       }
 
-      Set rec_path = start_path.intersection(end_path);
-      ingoing_plus = ingoing_plus.cup(subE_map.preImage(subE_map.image(rec_path)));
+      else {
+        for (unsigned int j = 0; j < aux_dsbg.Vmap().size(); ++j) {
+          start_path = start_path.cup(start_ith);
+          end_path = end_path.cup(end_ith);
+          if (!start_ith.intersection(first).isEmpty())
+            break;
+
+          v_start = aux_dsbg.mapD().image(start_ith); // Ending vertex of start_ith
+          start_ith = aux_dsbg.mapB().preImage(v_start); // Following edge from the start_ith
+          v_end = aux_dsbg.mapB().image(end_ith); // Starting vertex of end_ith
+          end_ith = aux_dsbg.mapD().preImage(v_end); // Following edge from end_ith
+        }
+
+        rec_path = start_path.intersection(end_path);
+        rec_path = subE_map.preImage(subE_map.image(rec_path));
+
+        if (first != rec_path) {
+          unvisited = unvisited.difference(rec_path);
+          Util::MD_NAT id(dims, k);
+          rec_map_ref().emplaceBack(SBGMap<Set>(rec_path, Exp(id)));
+          ++k;
+        }
+      }
     }
   }
+}
+
+template<typename Set>
+Set SBGTopSort<Set>::topSortStep()
+{
+  // Vertices without outgoing edges
+  Set nth = disordered().difference(mapB().image());
+
+  // Ingoing edges to nth
+  Set ingoing = mapD().preImage(nth);
+  Set ingoing_plus = ingoing.cup(rec_map().preImage(rec_map().image(ingoing)));
 
   Set domB_minus = mapB().dom().difference(ingoing_plus);
   Set domD_minus = mapD().dom().difference(ingoing_plus);
   set_mapB(mapB().restrict(domB_minus));
   set_mapD(mapD().restrict(domD_minus));
 
   set_disordered(disordered().difference(nth));
-  set_deletedE(deletedE().cup(ingoing_plus));
 
   return nth;
 }
@@ -926,6 +948,9 @@ VertexOrder<Set> SBGTopSort<Set>::calculate()
 
   VertexOrder<Set> res;
 
+  findRecursions();
+  if (debug())
+    Util::SBG_LOG << "rec_map: " << rec_map() << "\n\n";
   do {
     Set nth = topSortStep();
     if (!nth.isEmpty())

sbg/sbg_algorithms.hpp

@@ -255,7 +255,10 @@ struct SBGTopSort {
   member_class(PW, mapD);
 
   member_class(Set, disordered);
-  member_class(Set, deletedE); // Deleted edges
+
+  // Auxiliary map where recursive paths are grouped as pertaining to the same
+  // subset-edge
+  member_class(PW, rec_map);
 
   member_class(bool, debug);
 
@@ -265,6 +268,7 @@ struct SBGTopSort {
   VertexOrder<Set> calculate(); 
 
   private:
+  void findRecursions(); // Find all possible recursions
   Set topSortStep();
 };
 

test/test1_ts.test

@@ -11,6 +11,20 @@ equation
 end Test1;
 */
 
+/* Ordered model:
+
+model Test1
+  constant Integer N = 100;
+  Real a[N], b[N], x[N];
+equation
+  for i in 1:N loop
+    a[i] = 2*x[i] - b[i]; // (1)
+    b[i] = 0.5*(a[i] + x[i]); // (1)
+  end for;
+  der(x[i]) = 1 - a[i];
+end Test1;
+*/
+
 N = 100000
 
 F1 = N

test/test2_ts.test

@@ -37,7 +37,7 @@ equation
 end Test2;
 */
 
-N = 10
+N = 100000
 
 F1 = N-2
 F2 = N-2+F1

test/test3_scc.test

@@ -17,7 +17,7 @@ equation
 end Test3;
 */
 
-N = 1000
+N = 10
 
 F1 = N-2
 F2 = N-2+F1

test/test3_ts.test

@@ -29,15 +29,18 @@ equation
   b[1] = 0;
   b[N] = x[N-1];
   a[1] = 1;
-  der(x[1]) = a[1] - x[1];
   a[2] = a[1] + b[1];
-  der(x[2]) = a[2] - x[2]; 
+  der(x[1]) = a[1] - x[1];
+  der(x[2]) = a[2] - x[2];
   a[3] = a[2] + b[2];
-  for i in 3:N-2 loop
-    der(x[i]) = a[i] - x[i];
+  der(x[3]) = a[3] - x[3];
+  for i in 3:N-1 loop
     a[i+1] = a[i] + b[i];
   end for;
-  der(x[N]) = a[N] - x[N];
+  for i in loop
+    der(x[i]) = a[i] - x[i];
+  end for;
+  der(x[N]) = a[N] - x[N]; 
 end Test3;
 */
 

test/topsort1.test

@@ -1,6 +1,6 @@
 // Topological sort test, with a simple recursion
 
-N = 10
+N = 10000
 
 V1 = 1
 V2 = N-1+V1

test/topsort2.test

@@ -0,0 +1,38 @@
+// Topological sort test, with a simple recursion, and a self recursion on the right
+
+N = 100000
+
+V1 = 1
+V2 = N-1+V1
+V3 = N+V2
+
+E1 = 1
+E2 = N-1+E1
+E3 = N-1+E2
+E4 = N-1+E3
+
+off3b = N-1
+off4b = r(V3, 1)-E4-1
+
+off1d = N
+off2d = N-1
+off3d = 1
+off4d = r(V3, 1)-E4
+
+V %= {[1:1:V1], [V1+1:1:V2], [V2+1:1:V3]};
+Vmap %= <<{[1:1:V1]} -> 0*x+1, {[V1+1:1:V2]} -> 0*x+2, {[V2+1:1:V3]} -> 0*x+3>>; 
+mapB %= <<{[1:1:E1]} -> 1*x+0, {[E1+1:1:E2]} -> 1*x+off2d, {[E2+1:1:E3]} -> 1*x-off3b 
+         , {[E3+1:1:E4]} -> 1*x+off4b>>;
+mapD %= <<{[1:1:E1]} -> 1*x+off1d, {[E1+1:1:E2]} -> 1*x+0, {[E2+1:1:E3]} -> 1*x+off3d
+          , {[E3+1:1:E4]} -> 1*x+off4d>>;
+Emap %= <<{[1:1:E1]} -> 0*x+1, {[E1+1:1:E2], [E2+1:1:E3]} -> 0*x+2, {[E3+1:1:E4]} -> 0*x+3>>;
+
+sort(
+V %= {[1:1:V1], [V1+1:1:V2], [V2+1:1:V3]};
+Vmap %= <<{[1:1:V1]} -> 0*x+1, {[V1+1:1:V2]} -> 0*x+2, {[V2+1:1:V3]} -> 0*x+3>>; 
+mapB %= <<{[1:1:E1]} -> 1*x+0, {[E1+1:1:E2]} -> 1*x+off2d, {[E2+1:1:E3]} -> 1*x-off3b 
+         , {[E3+1:1:E4]} -> 1*x+off4b>>;
+mapD %= <<{[1:1:E1]} -> 1*x+off1d, {[E1+1:1:E2]} -> 1*x+0, {[E2+1:1:E3]} -> 1*x+off3d
+          , {[E3+1:1:E4]} -> 1*x+off4d>>;
+Emap %= <<{[1:1:E1]} -> 0*x+1, {[E1+1:1:E2], [E2+1:1:E3]} -> 0*x+2, {[E3+1:1:E4]} -> 0*x+3>>;
+)