From 4c8bad0ccd93f0f16df36f421566063228721167 Mon Sep 17 00:00:00 2001
From: Gard Spreemann <gspreemann@gmail.com>
Date: Wed, 7 Nov 2018 17:17:34 +0100
Subject: Massively speed up the MM version (reasonable sparse MM).

---
 include/gudhi_laplacian/simplex_tree.hpp | 46 ++++++++++++--------------------
 1 file changed, 17 insertions(+), 29 deletions(-)

diff --git a/include/gudhi_laplacian/simplex_tree.hpp b/include/gudhi_laplacian/simplex_tree.hpp
index 6d5aff5..179362f 100644
--- a/include/gudhi_laplacian/simplex_tree.hpp
+++ b/include/gudhi_laplacian/simplex_tree.hpp
@@ -1,5 +1,10 @@
 #pragma once
 
+
+/*
+WARNING: DO NOT USE FOR WEIGHTED STUFF YET! ALL WEIGHTS MUST BE 1!
+*/
+
 #include <vector>
 #include <algorithm>
 #include <cassert>
@@ -89,6 +94,7 @@ namespace Gudhi_laplacian
     std::vector<Sparse_vector<int> > cbd_same;
     std::vector<Sparse_vector<int> > cbd_same_t;
     std::vector<Sparse_vector<int> > cbd_below;
+    std::vector<Sparse_vector<int> > cbd_below_t;
     
     for (int d = from_degree; d < to_degree; ++d)
     {
@@ -99,12 +105,15 @@ namespace Gudhi_laplacian
       if (d == 0)
       {
         cbd_below = std::vector<Sparse_vector<int> >(n); // An nx0 matrix to avoid special-casing further down.
+        cbd_below_t = Sparse_vector<int>::transpose(cbd_below, 0);
       }
       else
       {
         cbd_below = cbd_same;
+        cbd_below_t = cbd_same_t;
       }
       assert(cbd_below.size() == n);
+      assert(cbd_below_t.size() == o);
 
       cbd_same.clear();
       for (unsigned int i = 0; i < m; ++i)
@@ -112,6 +121,8 @@ namespace Gudhi_laplacian
         cbd_same.push_back(assemble_cbd_row(st, stratification[d+1][i]));
       }
       cbd_same_t = Sparse_vector<int>::transpose(cbd_same, n);
+      assert(cbd_same.size() == m);
+      assert(cbd_same_t.size() == n);
 
       std::vector<Sparse_vector<double> > lap;
       unsigned int nnz = 0;
@@ -121,43 +132,20 @@ namespace Gudhi_laplacian
         lap.push_back(Sparse_vector<double>());
 
         // Up part.
-        for (unsigned int j = 0; j < n; ++j)
+        for (auto it = cbd_same_t[i].cbegin(); it != cbd_same_t[i].cend(); ++it)
         {
-          auto it = cbd_same_t[i].cbegin();
-          auto jt = cbd_same_t[j].cbegin();
-          while (it != cbd_same_t[i].cend() && jt != cbd_same_t[j].cend())
+          for (auto jt = cbd_same[it->first].cbegin(); jt != cbd_same[it->first].cend(); ++jt)
           {
-            if (it->first < jt->first)
-              ++it;
-            else if (jt->first < it->first)
-              ++jt;
-            else
-            {
-              double foo = (weights[d+1][jt->first]/weights[d][i])*(it->second * jt->second);
-              lap[i].add(j, (weights[d+1][jt->first]/weights[d][i])*(it->second * jt->second));
-              ++it;
-              ++jt;
-            }
+            lap[i].add(jt->first, (weights[d+1][it->first]/weights[d][jt->first]) * (it->second * jt->second));
           }
         }
 
         // Down part.
-        for (unsigned int j = 0; j < n; ++j)
+        for (auto it = cbd_below[i].cbegin(); it != cbd_below[i].cend(); ++it)
         {
-          auto it = cbd_below[i].cbegin();
-          auto jt = cbd_below[j].cbegin();
-          while (it != cbd_below[i].cend() && jt != cbd_below[j].cend())
+          for (auto jt = cbd_below_t[it->first].cbegin(); jt != cbd_below_t[it->first].cend(); ++jt)
           {
-            if (it->first < jt->first)
-              ++it;
-            else if (jt->first < it->first)
-              ++jt;
-            else
-            {
-              lap[i].add(j, (weights[d][j]/weights[d-1][it->first])*(it->second * jt->second));
-              ++it;
-              ++jt;
-            }
+            lap[i].add(jt->first, (weights[d][jt->first]/weights[d-1][it->first]) * (it->second * jt->second));
           }
         }
 
-- 
cgit v1.2.3