From 246b04540647e6523edc9d97b37705ac742b39ea Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?David=20Holzm=C3=BCller?= Date: Thu, 11 Apr 2019 18:50:10 +0200 Subject: [PATCH] MultiDimVector now provides a begin() end() interface with a corresponding iterator, also added some operations for vector arithmetic and reconstruction error check --- src/Interpolation.hpp | 109 ++++++++++++++++++++++++++++--- src/Iterators.hpp | 148 ++++++------------------------------------ test/main.cpp | 10 ++- 3 files changed, 130 insertions(+), 137 deletions(-) diff --git a/src/Interpolation.hpp b/src/Interpolation.hpp index 1b50378..15594d4 100644 --- a/src/Interpolation.hpp +++ b/src/Interpolation.hpp @@ -19,6 +19,7 @@ limitations under the License. #include #include #include +#include "Iterators.hpp" namespace fsi { @@ -48,6 +49,47 @@ template class MultiDimVector { It it; + class Iterator { + MultiDimVector &v; + StepIterator stepIt; + + public: + Iterator(MultiDimVector &v, It it) : v(v), stepIt(it){}; + + typedef Iterator self_type; + typedef double value_type; + typedef double &reference; + typedef double *pointer; + typedef std::forward_iterator_tag iterator_category; + typedef int difference_type; + + std::vector index() { return stepIt.index(); } + + // post-increment + self_type operator++() { + self_type i = *this; + stepIt.next(); + return i; + } + + // pre-increment + self_type operator++(int junk) { + stepIt.next(); + return *this; + } + reference operator*() { return v.data[stepIt.firstIndex()][stepIt.tailDimsCounter()]; } + pointer operator->() { return &(*(*this)); } + bool operator==(const self_type &rhs) { + if (stepIt.valid() != rhs.stepIt.valid()) { + return false; + } + + return stepIt.valid() or (stepIt.firstIndex() == rhs.stepIt.firstIndex() && + stepIt.tailDimsCounter() == rhs.stepIt.tailDimsCounter()); + } + bool operator!=(const self_type &rhs) { return !(*this == rhs); } + }; + public: // put the first dimension into an outer vector for processing reasons std::vector> data; @@ -59,7 +101,10 @@ class MultiDimVector { } }; - void swap(MultiDimVector &other) { data.swap(other.data); } + void swap(MultiDimVector &other) { + data.swap(other.data); + std::swap(it, other.it); + } void clear() { for (size_t dim = 0; dim < data.size(); ++dim) { @@ -68,8 +113,58 @@ class MultiDimVector { } It getJumpIterator() const { return it; } + + Iterator begin() { return Iterator(*this, it); } + + Iterator end() { + It end_it = it; + end_it.goToEnd(); + return Iterator(*this, end_it); + } + + MultiDimVector &operator+=(MultiDimVector const &other) { + for (size_t i = 0; i < data.size(); ++i) { + for (size_t j = 0; j < data[i].size(); ++j) { + data[i][j] += other.data[i][j]; + } + } + return *this; + } + + MultiDimVector &operator-=(MultiDimVector const &other) { + for (size_t i = 0; i < data.size(); ++i) { + for (size_t j = 0; j < data[i].size(); ++j) { + data[i][j] -= other.data[i][j]; + } + } + return *this; + } }; +template +MultiDimVector operator+(MultiDimVector first, MultiDimVector const &second) { + first += second; + return first; +} + +template +MultiDimVector operator-(MultiDimVector first, MultiDimVector const &second) { + first -= second; + return first; +} + +template +double squared_l2_norm(MultiDimVector const &v) { + double sum = 0.0; + for (size_t i = 0; i < v.data.size(); ++i) { + for (size_t j = 0; j < v.data[i].size(); ++j) { + double value = v.data[i][j]; + sum += value * value; + } + } + return sum; +} + template void multiply_lower_triangular_inplace(std::vector> L, MultiDimVector &v) { @@ -91,7 +186,7 @@ void multiply_lower_triangular_inplace(std::vector= data_size) { @@ -144,7 +238,7 @@ void multiply_upper_triangular_inplace(std::vector= data_size) { @@ -317,7 +410,7 @@ MultiDimVector evaluateFunction(It it, Func f, X x) { it.reset(); std::vector point(d); - while (not it.done()) { + while (it.valid()) { size_t last_dim_count = it.lastDimensionCount(); for (size_t dim = 0; dim < d - 1; ++dim) { point[dim] = x[dim](it.indexAt(dim)); diff --git a/src/Iterators.hpp b/src/Iterators.hpp index 4f7fc20..241cb70 100644 --- a/src/Iterators.hpp +++ b/src/Iterators.hpp @@ -56,7 +56,7 @@ class StepIterator { } } - bool done() { return jump_it.done(); } + bool valid() const { return jump_it.valid(); } void reset() { jump_it.reset(); @@ -65,75 +65,14 @@ class StepIterator { last_dim_count = jump_it.lastDimensionCount(); tail_dims_counter = 0; } -}; -template -class TemplateBoundedSumIterator { - size_t bound; - std::vector index_head; // contains all entries of the index except the last one - size_t index_head_sum; - bool is_done; + size_t firstIndex() const { return first_dim_value; } + size_t tailDimsCounter() const { return tail_dims_counter; } - public: - TemplateBoundedSumIterator(size_t bound) - : bound(bound), index_head(d - 1, 0), index_head_sum(0), is_done(false){}; - - /** - * At the current multi-index (i_1, ..., i_{d-1}, 0), return how many multi-indices starting with - * (i_1, ..., i_{d-1}) are contained in the multi-index set, then advance to the next multi-index - * that ends with a zero. - */ - size_t lastDimensionCount() { return bound - index_head_sum + 1; } - - void next() { - if (bound > index_head_sum) { - index_head_sum += 1; - index_head[d - 2] += 1; - } else { - int dim = d - 2; - - for (; dim >= 0 && index_head[dim] == 0; --dim) { - // reduce dimension until entry is nonzero - } - - if (dim > 0) { - index_head_sum -= (index_head[dim] - 1); - index_head[dim] = 0; - index_head[dim - 1] += 1; - } else if (dim == 0) { - index_head[dim] = 0; - index_head_sum = 0; - is_done = true; - } - } - } - - size_t firstIndex() const { return index_head[0]; } - - size_t indexAt(size_t dim) const { return index_head[dim]; } - - bool done() const { return is_done; } - - void reset() { - index_head = std::vector(d - 1, 0); - index_head_sum = 0; - is_done = false; - } - - size_t dim() const { return d; } - - std::vector indexBounds() const { return std::vector(d, bound + 1); } - - size_t numValues() const { return binom(bound + d, d); } - - /** - * Returns an iterator where the last index moves to the front. For an index set defined by a sum - * bound, nothing changes. - */ - TemplateBoundedSumIterator cycle() const { - TemplateBoundedSumIterator it = *this; - it.reset(); - return it; + std::vector index() const { + std::vector result = jump_it.getIndexHead(); + result.push_back(last_dim_value); + return result; } }; @@ -142,11 +81,11 @@ class BoundedSumIterator { size_t bound; std::vector index_head; // contains all entries of the index except the last one size_t index_head_sum; - bool is_done; + bool is_valid; public: BoundedSumIterator(size_t d, size_t bound) - : d(d), bound(bound), index_head(d - 1, 0), index_head_sum(0), is_done(false){}; + : d(d), bound(bound), index_head(d - 1, 0), index_head_sum(0), is_valid(true){}; /** * At the current multi-index (i_1, ..., i_{d-1}, 0), return how many multi-indices starting with @@ -171,9 +110,9 @@ class BoundedSumIterator { index_head[dim] = 0; index_head[dim - 1] += 1; } else if (dim == 0) { - index_head[dim] = 0; + index_head[0] = 0; index_head_sum = 0; - is_done = true; + is_valid = false; } } } @@ -182,12 +121,14 @@ class BoundedSumIterator { size_t indexAt(size_t dim) const { return index_head[dim]; } - bool done() const { return is_done; } + std::vector getIndexHead() const { return index_head; } + + bool valid() const { return is_valid; } void reset() { index_head = std::vector(d - 1, 0); index_head_sum = 0; - is_done = false; + is_valid = true; } size_t dim() const { return d; } @@ -206,6 +147,12 @@ class BoundedSumIterator { return result; } + void goToEnd() { + index_head = std::vector(d - 1, 0); + index_head_sum = 0; + is_valid = false; + } + /** * Returns an iterator where the last index moves to the front. For an index set defined by a sum * bound, nothing changes. @@ -216,57 +163,4 @@ class BoundedSumIterator { return it; } }; - -class StandardBoundedSumIterator { - size_t d; - size_t bound; - std::vector index; // contains all entries of the index except the last one - size_t index_sum; - bool is_done; - - public: - StandardBoundedSumIterator(size_t d, size_t bound) - : d(d), bound(bound), index(d, 0), index_sum(0), is_done(false){}; - - /** - * At the current multi-index (i_1, ..., i_{d-1}, 0), return how many multi-indices starting with - * (i_1, ..., i_{d-1}) are contained in the multi-index set, then advance to the next multi-index - * that ends with a zero. - */ - bool next() { - if (bound > index_sum) { - index_sum += 1; - index[d - 1] += 1; - } else { - int dim = d - 1; - - for (; dim >= 0 && index[dim] == 0; --dim) { - // reduce dimension until entry is nonzero - } - - if (dim > 0) { - index_sum -= index[dim]; - index[dim] = 0; - index[dim - 1] += 1; - } else if (dim == 0) { - index[dim] = 0; - index_sum = 0; - is_done = true; - } - } - return is_done; - } - - size_t firstIndex() { return index[0]; } - - size_t indexSum() { return index_sum; } - - bool done() { return is_done; } - - void reset() { - index = std::vector(d - 1, 0); - index_sum = 0; - is_done = false; - } -}; } diff --git a/test/main.cpp b/test/main.cpp index 3458dec..1a2bd12 100644 --- a/test/main.cpp +++ b/test/main.cpp @@ -26,8 +26,8 @@ limitations under the License. */ void runFunctions() { - constexpr size_t d = 8; - size_t bound = 24; + constexpr size_t d = 30; + size_t bound = 8; // fsi::TemplateBoundedSumIterator it(bound); fsi::BoundedSumIterator it(d, bound); std::vector phi(d); @@ -56,6 +56,12 @@ void runFunctions() { std::cout << "Number of points: " << it.numValues() << "\n"; + std::cout << "Reconstruction error (L2 norm): " << sqrt(fsi::squared_l2_norm(b - rhs)) << "\n"; + + // for (auto it = c.begin(); it != c.end(); ++it) { + // std::cout << "Value at index " << it.index() << ": " << *it << "\n"; + // } + // for (size_t i = 0; i < result.data.size(); ++i) { // std::cout << result.data[i] << "\n\n"; // }