added files for performance measurement and plotting

.gitignore

@@ -13,6 +13,7 @@
 *.auxlock
 *.dpth
 *.md5
+*.csv
 
 build/*
 TexCode/*

src/Iterators.hpp

@@ -9,7 +9,7 @@
 
 namespace fsi {
 
-size_t binom(size_t n, size_t k) {
+inline size_t binom(size_t n, size_t k) {
   if (2 * k > n) {
     k = n - k;
   }
@@ -26,6 +26,47 @@ size_t binom(size_t n, size_t k) {
   return prod;
 }
 
+template <typename JumpIt>
+class StepIterator {
+  JumpIt jump_it;
+  size_t first_dim_value;
+  size_t last_dim_value;
+  size_t last_dim_count;
+  size_t tail_dims_counter;
+
+ public:
+  StepIterator(JumpIt jump_it)
+      : jump_it(jump_it),
+        first_dim_value(0),
+        last_dim_value(0),
+        last_dim_count(jump_it.lastDimensionCount()),
+        tail_dims_counter(0){};
+
+  void next() {
+    tail_dims_counter += 1;
+    last_dim_value += 1;
+    if (last_dim_value >= last_dim_count) {
+      last_dim_value = 0;
+      jump_it.next();
+      if (jump_it.firstIndex() != first_dim_value) {
+        first_dim_value = jump_it.firstIndex();
+        tail_dims_counter = 0;
+      }
+      last_dim_count = jump_it.lastDimensionCount();
+    }
+  }
+
+  bool done() { return jump_it.done(); }
+
+  void reset() {
+    jump_it.reset();
+    first_dim_value = 0;
+    last_dim_value = 0;
+    last_dim_count = jump_it.lastDimensionCount();
+    tail_dims_counter = 0;
+  }
+};
+
 template <size_t d>
 class TemplateBoundedSumIterator {
   size_t bound;

test/common.hpp

@@ -0,0 +1,106 @@
+// Copyright (C) 2008-today The SG++ project
+// This file is part of the SG++ project. For conditions of distribution and
+// use, please see the copyright notice provided with SG++ or at
+// sgpp.sparsegrids.org
+
+#pragma once
+
+#include <Interpolation.hpp>
+#include <Iterators.hpp>
+#include <chrono>
+#include <cmath>
+#include <iostream>
+#include <vector>
+
+inline double f(std::vector<double> point) {
+  double prod = 1.0;
+  for (size_t dim = 0; dim < point.size(); ++dim) {
+    prod *= point[dim];
+  }
+  return prod;
+}
+
+class GoldenPointDistribution {
+  static constexpr double golden_ratio = 0.5 * (1.0 + sqrt(5));
+
+ public:
+  double operator()(size_t idx) {
+    double value = (idx + 1) * golden_ratio;
+    return value - int(value);
+  }
+};
+
+class SimplePointDistribution {
+ public:
+  double operator()(size_t idx) {
+    if (idx == 0) {
+      return 0.0;
+    } else if (idx == 1) {
+      return 1.0;
+    } else {
+      return 0.5;
+    }
+  }
+};
+
+class MonomialFunctions {
+ public:
+  std::function<double(double)> operator()(size_t idx) {
+    return [=](double x) { return pow(x, idx); };
+  }
+};
+
+template <class T>
+inline std::ostream &operator<<(std::ostream &os, const std::vector<T> &v) {
+  os << "[";
+  for (auto ii = v.begin(); ii != v.end(); ++ii) {
+    os << " " << *ii;
+  }
+  os << " ]";
+  return os;
+}
+
+inline std::ostream &operator<<(std::ostream &os, fsi::MultiDimVector const &v) {
+  for (size_t i = 0; i < v.data.size(); ++i) {
+    std::cout << v.data[i] << "\n\n";
+  }
+
+  return os;
+}
+
+// TODO: refactoring:
+/**
+ * - Interface for MultiDimVector
+ * - Separate Functions for L- and U-Multiplication, Creation of LU decomposition, computation of
+ * function values
+ * - Pass a callback function to iterator instead of calling it.next() - this might improve the
+ * vector functions (could be made recursive or even loops for fixed dimension implementation)
+ * - Typed interface?
+ * - Tests?
+ * - More point distributions / Basis functions?
+ * - Forward evaluation?
+ * - Computation of derivatives?
+ */
+
+inline double measure_execution_time(std::function<void()> f) {
+  auto start = std::chrono::high_resolution_clock::now();
+  f();
+  auto finish = std::chrono::high_resolution_clock::now();
+  std::chrono::duration<double> elapsed = finish - start;
+  return elapsed.count();
+}
+
+class Timer {
+  std::chrono::system_clock::time_point start;
+
+ public:
+  Timer() : start(std::chrono::high_resolution_clock::now()){};
+  void reset() { start = std::chrono::high_resolution_clock::now(); }
+  double elapsed() {
+    auto finish = std::chrono::high_resolution_clock::now();
+    std::chrono::duration<double> elapsed = finish - start;
+    return elapsed.count();
+  }
+};
+
+void measurePerformance();

test/main.cpp

@@ -13,105 +13,9 @@ See the License for the specific language governing permissions and
 limitations under the License.
 ==============================================================================*/
 
-#include <Interpolation.hpp>
-#include <Iterators.hpp>
-#include <chrono>
-#include <cmath>
-#include <iostream>
+#include "common.hpp"
 
-double f(std::vector<double> point) {
-  double prod = 1.0;
-  for (size_t dim = 0; dim < point.size(); ++dim) {
-    prod *= point[dim];
-  }
-  return prod;
-}
-
-class GoldenPointDistribution {
-  static constexpr double golden_ratio = 0.5 * (1.0 + sqrt(5));
-
- public:
-  double operator()(size_t idx) {
-    double value = (idx + 1) * golden_ratio;
-    return value - int(value);
-  }
-};
-
-class SimplePointDistribution {
- public:
-  double operator()(size_t idx) {
-    if (idx == 0) {
-      return 0.0;
-    } else if (idx == 1) {
-      return 1.0;
-    } else {
-      return 0.5;
-    }
-  }
-};
-
-class MonomialFunctions {
- public:
-  std::function<double(double)> operator()(size_t idx) {
-    return [=](double x) { return pow(x, idx); };
-  }
-};
-
-template <class T>
-inline std::ostream &operator<<(std::ostream &os, const std::vector<T> &v) {
-  os << "[";
-  for (auto ii = v.begin(); ii != v.end(); ++ii) {
-    os << " " << *ii;
-  }
-  os << " ]";
-  return os;
-}
-
-std::ostream &operator<<(std::ostream &os, fsi::MultiDimVector const &v) {
-  for (size_t i = 0; i < v.data.size(); ++i) {
-    std::cout << v.data[i] << "\n\n";
-  }
-
-  return os;
-}
-
-// TODO: refactoring:
-/**
- * - Interface for MultiDimVector
- * - Separate Functions for L- and U-Multiplication, Creation of LU decomposition, computation of
- * function values
- * - Pass a callback function to iterator instead of calling it.next() - this might improve the
- * vector functions (could be made recursive or even loops for fixed dimension implementation)
- * - Typed interface?
- * - Tests?
- * - More point distributions / Basis functions?
- * - Forward evaluation?
- * - Computation of derivatives?
- */
-
-double measure_execution_time(std::function<void()> f) {
-  auto start = std::chrono::high_resolution_clock::now();
-  f();
-  auto finish = std::chrono::high_resolution_clock::now();
-  std::chrono::duration<double> elapsed = finish - start;
-  return elapsed.count();
-}
-
-class Timer {
-  std::chrono::system_clock::time_point start;
-
- public:
-  Timer() : start(std::chrono::high_resolution_clock::now()){};
-  void reset() { start = std::chrono::high_resolution_clock::now(); }
-  double elapsed() {
-    auto finish = std::chrono::high_resolution_clock::now();
-    std::chrono::duration<double> elapsed = finish - start;
-    return elapsed.count();
-  }
-};
-
-// using namespace fsi;
-int main() {
+void runFunctions() {
   constexpr size_t d = 8;
   size_t bound = 24;
   fsi::TemplateBoundedSumIterator<d> it(bound);
@@ -146,3 +50,9 @@ int main() {
   //   std::cout << result.data[i] << "\n\n";
   // }
 }
+
+// using namespace fsi;
+int main() {
+  // runFunctions();
+  measurePerformance();
+}

test/performance.cpp

@@ -0,0 +1,81 @@
+/* Copyright 2019 The fast_sparse_interpolation Authors. All Rights Reserved.
+
+Licensed under the Apache License, Version 2.0 (the "License");
+you may not use this file except in compliance with the License.
+You may obtain a copy of the License at
+
+    http://www.apache.org/licenses/LICENSE-2.0
+
+Unless required by applicable law or agreed to in writing, software
+distributed under the License is distributed on an "AS IS" BASIS,
+WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+See the License for the specific language governing permissions and
+limitations under the License.
+==============================================================================*/
+
+#include "common.hpp"
+
+#include <fstream>
+
+double measureRuntime(std::function<double()> f, size_t k = 3) {
+  size_t n = 2 * k + 1;
+
+  std::vector<double> singleResults(n);
+  for (size_t i = 0; i < n; ++i) {
+    singleResults[i] = f();
+  }
+  std::sort(singleResults.begin(), singleResults.end());
+  // return (singleResults[k-1] + singleResults[k] +
+  // singleResults[k+1]) / 3.0;
+  return singleResults[k];
+}
+
+void measurePerformance() {
+  // std::vector < size_t > dimensions = {2, 3, 4, 6, 8, 12, 16, 24, 32, 48, 64 };
+  // std::vector<size_t> dimensions = {5, 10};
+
+  std::vector<size_t> dimensions = {2, 4, 8, 16, 32};
+  // size_t maxNumPoints = 20000000;
+  size_t maxNumPoints = 500000;
+
+  std::ostringstream stream;
+
+  size_t approxNumSteps = 8;
+  size_t k = 1;
+
+  for (size_t d : dimensions) {
+    std::cout << "Dimension: " << d << "\n";
+    size_t maxBound = 0;
+
+    while (fsi::binom(maxBound + d, d) <= maxNumPoints) {
+      maxBound += 1;
+    }
+
+    size_t stepsize = maxBound / approxNumSteps + 1;
+
+    for (size_t bound = 2; bound < maxBound; bound += stepsize) {
+      std::cout << "Bound: " << bound << "\n";
+      double runtime = measureRuntime(
+          [&]() {
+            fsi::BoundedSumIterator it(d, bound);
+            std::vector<MonomialFunctions> phi(d);
+            std::vector<GoldenPointDistribution> x(d);
+
+            auto rhs = evaluateFunction(it, f, x);
+            auto op = createInterpolationOperator(it, phi, x);
+
+            Timer timer;
+            op.solve(rhs);
+            double time = timer.elapsed();
+            std::cout << "Time for solve(): " << time << " s\n";
+            return time;
+          },
+          k);
+
+      stream << d << ", " << bound << ", " << fsi::binom(bound + d, d) << ", " << runtime << "\n";
+    }
+  }
+
+  std::ofstream file("performance_data.csv");
+  file << stream.str();
+}

test/plot.py

@@ -0,0 +1,32 @@
+import matplotlib.pyplot as plt
+from mpl_toolkits.mplot3d import Axes3D
+import numpy as np
+import os.path
+import math
+
+
+def readFromFile(filename):
+    if not os.path.isfile(filename):
+        return ''
+
+    file = open(filename, 'r')
+    result = file.read()
+    file.close()
+    return result
+
+datastr = readFromFile("../performance_data.csv")
+rows = datastr.split("\n")[:-1]
+values = np.array([[float(s) for s in r.split(', ')] for r in rows])
+# print(values)
+
+dimensions = values[:, 0]
+bounds = values[:, 1]
+points = values[:, 2]
+times = values[:, 3]
+
+fig = plt.figure()
+ax = fig.add_subplot(1,1,1, projection='3d')
+surf = ax.plot_trisurf(np.log2(dimensions), np.log2(points), np.log10(times), cmap='viridis', edgecolor='none')
+plt.show()
+
+