From 1f457040a6769af23ca5042da1de9d5562e79460 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?David=20Holzm=C3=BCller?= Date: Sun, 7 Apr 2019 21:54:47 +0200 Subject: [PATCH] moved iterators to separate headers, moved vector index counting to iterator --- src/Interpolation.hpp | 372 ++++-------------------------------------- src/Iterators.hpp | 238 +++++++++++++++++++++++++++ test/main.cpp | 1 + 3 files changed, 269 insertions(+), 342 deletions(-) create mode 100644 src/Iterators.hpp diff --git a/src/Interpolation.hpp b/src/Interpolation.hpp index 7c3014d..5192084 100644 --- a/src/Interpolation.hpp +++ b/src/Interpolation.hpp @@ -22,215 +22,27 @@ limitations under the License. namespace fsi { -size_t binom(size_t n, size_t k) { - if (2 * k > n) { - k = n - k; +inline std::ostream &operator<<(std::ostream &os, boost::numeric::ublas::matrix matrix) { + os << "["; + for (size_t i = 0; i < matrix.size1(); ++i) { + for (size_t j = 0; j < matrix.size2(); ++j) { + os << matrix(i, j) << " "; + } + os << "\n"; } - - size_t prod = 1; - size_t upper_factor = n + 1 - k; - size_t lower_factor = 1; - while (upper_factor <= n) { - prod *= upper_factor; - prod /= lower_factor; - upper_factor += 1; - lower_factor += 1; - } - return prod; + os << " ]"; + return os; } -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; - - 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; - } - } +template +inline std::ostream &operator<<(std::ostream &os, const std::vector &v) { + os << "["; + for (auto ii = v.begin(); ii != v.end(); ++ii) { + os << " " << *ii; } - - 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; - } -}; - -class BoundedSumIterator { - size_t d; - 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; - - public: - BoundedSumIterator(size_t d, size_t bound) - : d(d), 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. - */ - BoundedSumIterator cycle() { - BoundedSumIterator it = *this; - it.reset(); - 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; - } -}; + os << " ]"; + return os; +} class MultiDimVector { public: @@ -265,15 +77,11 @@ void multiply_lower_triangular_inplace(It it, std::vector= data_size) { - second_v_index = 0; - first_v_index += 1; - data_pointer = &v.data[first_v_index][0]; - data_size = v.data[first_v_index].size(); - } - it.next(); + bool first_dimension_changed = it.next(); + if (first_dimension_changed) { + data_pointer = &v.data[it.firstIndex()][0]; + } } v.swap(w); @@ -320,15 +124,11 @@ void multiply_upper_triangular_inplace(It it, std::vector= data_size) { - second_v_index = 0; - first_v_index += 1; - data_pointer = &v.data[first_v_index][0]; - data_size = v.data[first_v_index].size(); - } - it.next(); + bool first_dimension_changed = it.next(); + if (first_dimension_changed) { + data_pointer = &v.data[it.firstIndex()][0]; + } } it = it.cycle(); @@ -517,119 +313,11 @@ MultiDimVector evaluateFunction(It it, Func f, X x) { return v; } -inline std::ostream &operator<<(std::ostream &os, boost::numeric::ublas::matrix matrix) { - os << "["; - for (size_t i = 0; i < matrix.size1(); ++i) { - for (size_t j = 0; j < matrix.size2(); ++j) { - os << matrix(i, j) << " "; - } - os << "\n"; - } - os << " ]"; - return os; -} - -template -inline std::ostream &operator<<(std::ostream &os, const std::vector &v) { - os << "["; - for (auto ii = v.begin(); ii != v.end(); ++ii) { - os << " " << *ii; - } - os << " ]"; - return os; -} - template MultiDimVector interpolate(Func f, It it, Phi phi, X x) { - auto n = it.indexBounds(); - size_t d = it.dim(); - - namespace ublas = boost::numeric::ublas; - typedef ublas::matrix Matrix; - - std::vector Linv, Uinv; - - // create matrices and inverted LU decompositions - for (size_t k = 0; k < d; ++k) { - // std::cout << "Matrix creation loop\n"; - Matrix Mk(n[k], n[k]); - for (size_t i = 0; i < n[k]; ++i) { - for (size_t j = 0; j < n[k]; ++j) { - Mk(i, j) = phi[k](j)(x[k](i)); - } - } - - // std::cout << "Matrices:\n"; - - // std::cout << Mk << "\n"; - - ublas::lu_factorize(Mk); - - // std::cout << Mk << "\n"; - - Matrix Lkinv = ublas::identity_matrix(n[k]); - Matrix Ukinv = ublas::identity_matrix(n[k]); - ublas::inplace_solve(Mk, Lkinv, ublas::unit_lower_tag()); - ublas::inplace_solve(Mk, Ukinv, ublas::upper_tag()); - - // std::cout << Lkinv << "\n"; - - // std::cout << Ukinv << "\n"; - - Linv.push_back(Lkinv); - Uinv.push_back(Ukinv); - } - - MultiDimVector v(n[0]); - - std::cout << "Compute function values\n"; - - // compute function values - it.reset(); - std::vector point(d); - while (not it.done()) { - size_t last_dim_count = it.lastDimensionCount(); - for (size_t dim = 0; dim < d - 1; ++dim) { - point[dim] = x[dim](it.indexAt(dim)); - } - - for (size_t last_dim_idx = 0; last_dim_idx < last_dim_count; ++last_dim_idx) { - point[d - 1] = x[d - 1](last_dim_idx); - - double function_value = f(point); - - v.data[it.firstIndex()].push_back(function_value); - } - - it.next(); - } - - size_t number = 0; - for (size_t dim = 0; dim < n[0]; ++dim) { - number += v.data[dim].size(); - } - std::cout << "number of points: " << number << "\n"; - std::cout << "alternative number of points: " << it.numValues() << "\n"; - - // for (size_t i = 0; i < v.data.size(); ++i) { - // std::cout << v.data[i] << "\n\n"; - // } - - std::cout << "First matrix multiplication\n"; - - // multiply by L^{-1} - - multiply_lower_triangular_inplace(it, Linv, v); - - std::cout << "Second matrix multiplication\n"; - - // multiply by U^{-1} - // the multiplication is based on a cyclic permutation of the indices: the last index of v becomes - // the first index of w - - multiply_upper_triangular_inplace(it, Uinv, v); - - return v; + auto rhs = evaluateFunction(it, f, x); + auto op = createInterpolationOperator(it, phi, x); + return op.solve(rhs); } } /* namespace fsi */ diff --git a/src/Iterators.hpp b/src/Iterators.hpp new file mode 100644 index 0000000..ef18dc2 --- /dev/null +++ b/src/Iterators.hpp @@ -0,0 +1,238 @@ +// 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 + +namespace fsi { + +size_t binom(size_t n, size_t k) { + if (2 * k > n) { + k = n - k; + } + + size_t prod = 1; + size_t upper_factor = n + 1 - k; + size_t lower_factor = 1; + while (upper_factor <= n) { + prod *= upper_factor; + prod /= lower_factor; + upper_factor += 1; + lower_factor += 1; + } + return prod; +} + +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; + size_t middle_dimensions_counter; + bool is_done; + + public: + TemplateBoundedSumIterator(size_t bound) + : bound(bound), + index_head(d - 1, 0), + index_head_sum(0), + middle_dimensions_counter(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; } + + bool next() { + middle_dimensions_counter += lastDimensionCount(); + + 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 > 1) { + index_head_sum -= (index_head[dim] - 1); + index_head[dim] = 0; + index_head[dim - 1] += 1; + } else if (dim == 1) { + index_head_sum -= (index_head[1] - 1); + index_head[1] = 0; + index_head[0] += 1; + middle_dimensions_counter = 0; + return true; + } else if (dim == 0) { + index_head[dim] = 0; + index_head_sum = 0; + is_done = true; + } + } + return false; + } + + 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; } + + size_t getMiddleDimensionsCounter() const { return middle_dimensions_counter; } + + void reset() { + index_head = std::vector(d - 1, 0); + index_head_sum = 0; + middle_dimensions_counter = 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; + } +}; + +class BoundedSumIterator { + size_t d; + 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; + + public: + BoundedSumIterator(size_t d, size_t bound) + : d(d), 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. + */ + BoundedSumIterator cycle() { + BoundedSumIterator it = *this; + it.reset(); + 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 5b139a8..2c9f2d3 100644 --- a/test/main.cpp +++ b/test/main.cpp @@ -14,6 +14,7 @@ limitations under the License. ==============================================================================*/ #include +#include #include #include #include