From 999de0cbf791e279e0dabcfa7582d442d0b6d743 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?David=20Holzm=C3=BCller?= Date: Sun, 7 Apr 2019 20:42:24 +0200 Subject: [PATCH] new interface, new timing --- CMakeLists.txt | 17 +- header.txt | 2 +- src/Interpolation.cpp | 8 - src/Interpolation.hpp | 426 ++++++++++++++++++++---------------------- test/main.cpp | 81 +++++++- 5 files changed, 275 insertions(+), 259 deletions(-) delete mode 100644 src/Interpolation.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 8e62184..ec690aa 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -16,21 +16,14 @@ file(GLOB INC_TEST "test/*.hpp") option(BUILD_TEST ON) ## Build application -add_library(fast_sparse_interpolation ${SRC}) +# add_library(fast_sparse_interpolation ${SRC}) include_directories("${PROJECT_SOURCE_DIR}/src") -if(BUILD_TEST) - if(BUILD_TEST STREQUAL "ON") - file(GLOB SRC_TEST "test/*.cpp") - add_executable(fsi-test ${SRC_TEST}) - include_directories("${PROJECT_SOURCE_DIR}/src") - target_link_libraries(fast_sparse_interpolation) - endif() -endif() +file(GLOB SRC_TEST "test/*.cpp") +add_executable(fsi-test ${SRC_TEST}) +include_directories("${PROJECT_SOURCE_DIR}/src") - - -install(TARGETS fast_sparse_interpolation DESTINATION lib) +# install(TARGETS fast_sparse_interpolation DESTINATION lib) install(FILES ${INC} DESTINATION include/sfcpp) diff --git a/header.txt b/header.txt index 4cb4aea..6a4b637 100644 --- a/header.txt +++ b/header.txt @@ -1,4 +1,4 @@ -/* Copyright 2017 The sfcpp Authors. All Rights Reserved. +/* 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. diff --git a/src/Interpolation.cpp b/src/Interpolation.cpp deleted file mode 100644 index 736cb9b..0000000 --- a/src/Interpolation.cpp +++ /dev/null @@ -1,8 +0,0 @@ -// 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 - -#include "Interpolation.hpp" - -namespace fsi {} /* namespace fsi */ diff --git a/src/Interpolation.hpp b/src/Interpolation.hpp index 1f794a6..7c3014d 100644 --- a/src/Interpolation.hpp +++ b/src/Interpolation.hpp @@ -1,7 +1,17 @@ -// 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 +/* 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. +==============================================================================*/ #pragma once @@ -29,18 +39,6 @@ size_t binom(size_t n, size_t k) { return prod; } -template -void iterateRecursive(size_t remaining, std::function callback) { - for (size_t val = 0; val <= remaining; ++val) { - iterateRecursive(remaining - val, callback); - } -} - -template <> -void iterateRecursive<0>(size_t remaining, std::function callback) { - callback(remaining); -} - template class TemplateBoundedSumIterator { size_t bound; @@ -52,14 +50,6 @@ class TemplateBoundedSumIterator { TemplateBoundedSumIterator(size_t bound) : bound(bound), index_head(d - 1, 0), index_head_sum(0), is_done(false){}; - void iterate(std::function callback, - std::function outerLoopCallback) { - for (size_t val = 0; val <= bound; ++val) { - outerLoopCallback(val); - iterateRecursive(bound - val, callback); - } - } - /** * 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 @@ -106,6 +96,8 @@ class TemplateBoundedSumIterator { 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. @@ -240,12 +232,6 @@ class StandardBoundedSumIterator { } }; -// class FastSparseInterpolation { -// public: -// FastSparseInterpolation(std::function)> const &f, -// BoundedSumIterator it, std::vector> phi); -//}; - class MultiDimVector { public: // put the first dimension into an outer vector for processing reasons @@ -261,28 +247,6 @@ class MultiDimVector { } }; -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 void multiply_lower_triangular_inplace(It it, std::vector> L, MultiDimVector &v) { @@ -393,6 +357,188 @@ void multiply_upper_triangular_inplace(It it, std::vector +class SparseTPOperator { + It it; + std::vector> M; + std::vector> LU; + std::vector> L; + std::vector> U; + std::vector> Linv; + std::vector> Uinv; + size_t d; + + public: + SparseTPOperator(It it, std::vector> matrices) + : it(it), M(matrices), d(matrices.size()){}; + + void prepareCommon() { + if (LU.size() > 0) { + return; // already prepared + } + + namespace ublas = boost::numeric::ublas; + typedef ublas::matrix Matrix; + + for (size_t k = 0; k < d; ++k) { + // std::cout << "Matrix creation loop\n"; + Matrix Mk = M[k]; + + ublas::lu_factorize(Mk); + + LU.push_back(Mk); + } + } + + void prepareApply() { + if (L.size() > 0) { + return; // already prepared + } + + prepareCommon(); + + namespace ublas = boost::numeric::ublas; + typedef ublas::matrix Matrix; + + for (size_t k = 0; k < d; ++k) { + size_t nk = M[k].size1(); + Matrix &LUk = LU[k]; + Matrix Lk(nk, nk); + Matrix Uk(nk, nk); + + for (size_t i = 0; i < nk; ++i) { + for (size_t j = 0; j < i; ++j) { + Lk(i, j) = LUk(i, j); + } + + Lk(i, i) = 1.0; + + for (size_t j = i; j < nk; ++j) { + Uk(i, j) = LUk(i, j); + } + } + + L.push_back(Lk); + U.push_back(Uk); + } + } + + void prepareSolve() { + if (Linv.size() > 0) { + return; // already prepared + } + + prepareCommon(); + + namespace ublas = boost::numeric::ublas; + typedef ublas::matrix Matrix; + + for (size_t k = 0; k < d; ++k) { + Matrix Lkinv = ublas::identity_matrix(M[k].size1()); + Matrix Ukinv = ublas::identity_matrix(M[k].size1()); + + ublas::inplace_solve(LU[k], Lkinv, ublas::unit_lower_tag()); + ublas::inplace_solve(LU[k], Ukinv, ublas::upper_tag()); + + Linv.push_back(Lkinv); + Uinv.push_back(Ukinv); + } + } + + MultiDimVector apply(MultiDimVector input) { + prepareApply(); + multiply_upper_triangular_inplace(it, U, input); + multiply_lower_triangular_inplace(it, L, input); + return input; + } + + MultiDimVector solve(MultiDimVector rhs) { + prepareSolve(); + multiply_lower_triangular_inplace(it, Linv, rhs); + multiply_upper_triangular_inplace(it, Uinv, rhs); + return rhs; + } +}; + +template +SparseTPOperator createInterpolationOperator(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 matrices; + + // 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)); + } + } + + matrices.push_back(Mk); + } + + return SparseTPOperator(it, matrices); +} + +template +MultiDimVector evaluateFunction(It it, Func f, X x) { + size_t d = it.dim(); + auto n = it.indexBounds(); + MultiDimVector v(n[0]); + + 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(); + } + + 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(); @@ -486,180 +632,4 @@ MultiDimVector interpolate(Func f, It it, Phi phi, X x) { return v; } -template -MultiDimVector interpolate2(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"; - - // for (size_t i = 0; i < v.data.size(); ++i) { - // std::cout << v.data[i] << "\n\n"; - // } - - size_t num_sum_op = 0; - - std::cout << "First matrix multiplication\n"; - - // multiply by L^{-1} - // the multiplication is based on a cyclic permutation of the indices: the last index of v becomes - // the first index of w - - for (int k = d - 1; k >= 0; --k) { - MultiDimVector w(n[k]); - for (size_t idx = 0; idx < n[k]; ++idx) { - // TODO: make compatible with other iterators - w.data[idx].reserve(v.data[idx].size()); - } - - auto &Lkinv = Linv[k]; - it.reset(); - - size_t second_v_index = 0; - - double *data_pointer = &v.data[0][0]; - - it.iterate( - [&](size_t last_dim_count) { - double *offset_data_pointer = data_pointer + second_v_index; - for (size_t i = 0; i < last_dim_count; ++i) { - double sum = 0.0; - for (size_t j = 0; j <= i; ++j) { - sum += Lkinv(i, j) * (*(offset_data_pointer + j)); - ++num_sum_op; - } - w.data[i].push_back(sum); - } - second_v_index += last_dim_count; - }, - [&](size_t first_dim_index) { - second_v_index = 0; - data_pointer = &v.data[first_dim_index][0]; - }); - - v.swap(w); - - it = it.cycle(); - - // for (size_t i = 0; i < v.data.size(); ++i) { - // std::cout << v.data[i] << "\n\n"; - // } - } - - 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 - - for (int k = d - 1; k >= 0; --k) { - MultiDimVector w(n[k]); - for (size_t idx = 0; idx < n[k]; ++idx) { - // TODO: make compatible with other iterators - w.data[idx].reserve(v.data[idx].size()); - } - - auto &Ukinv = Uinv[k]; - it.reset(); - - size_t second_v_index = 0; - - double *data_pointer = &v.data[0][0]; - - it.iterate( - [&](size_t last_dim_count) { - double *offset_data_pointer = data_pointer + second_v_index; - for (size_t i = 0; i < last_dim_count; ++i) { - double sum = 0.0; - for (size_t j = i; j < last_dim_count; ++j) { - sum += Ukinv(i, j) * (*(offset_data_pointer + j)); - ++num_sum_op; - } - w.data[i].push_back(sum); - } - second_v_index += last_dim_count; - }, - [&](size_t first_dim_index) { - second_v_index = 0; - data_pointer = &v.data[first_dim_index][0]; - }); - - v.swap(w); - - it = it.cycle(); - - // for (size_t i = 0; i < v.data.size(); ++i) { - // std::cout << v.data[i] << "\n\n"; - // } - } - - std::cout << "Total number of summation operations: " << num_sum_op << "\n"; - - return v; -} - } /* namespace fsi */ diff --git a/test/main.cpp b/test/main.cpp index 1a7e092..5b139a8 100644 --- a/test/main.cpp +++ b/test/main.cpp @@ -1,3 +1,17 @@ +/* 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 #include @@ -52,6 +66,14 @@ inline std::ostream &operator<<(std::ostream &os, const std::vector &v) { 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 @@ -66,19 +88,58 @@ inline std::ostream &operator<<(std::ostream &os, const std::vector &v) { * - Computation of derivatives? */ -// using namespace fsi; -int main() { - constexpr size_t d = 5; - size_t bound = 57; - // fsi::TemplateBoundedSumIterator it(bound); - fsi::BoundedSumIterator it(d, bound); - std::vector phi(d); - std::vector x(d); +double measure_execution_time(std::function f) { auto start = std::chrono::high_resolution_clock::now(); - auto result = fsi::interpolate(f, it, phi, x); + f(); auto finish = std::chrono::high_resolution_clock::now(); std::chrono::duration elapsed = finish - start; - std::cout << "Elapsed time: " << elapsed.count() << " s\n"; + 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 elapsed = finish - start; + return elapsed.count(); + } +}; + +// using namespace fsi; +int main() { + constexpr size_t d = 8; + size_t bound = 24; + fsi::TemplateBoundedSumIterator it(bound); + // fsi::BoundedSumIterator it(d, bound); + std::vector phi(d); + std::vector x(d); + + auto rhs = evaluateFunction(it, f, x); + auto op = createInterpolationOperator(it, phi, x); + + // std::cout << rhs << "\n"; + + Timer timer; + op.prepareSolve(); + // auto result = fsi::interpolate(f, it, phi, x); + std::cout << "Time for prepareSolve(): " << timer.elapsed() << " s\n"; + + timer.reset(); + auto c = op.solve(rhs); + // auto c = fsi::interpolate(f, it, phi, x); + std::cout << "Time for solve(): " << timer.elapsed() << " s\n"; + // std::cout << c << "\n"; + + timer.reset(); + auto b = op.apply(c); + std::cout << "Time for apply(): " << timer.elapsed() << " s\n"; + // std::cout << b << "\n"; + + std::cout << "Number of points: " << it.numValues() << "\n"; // for (size_t i = 0; i < result.data.size(); ++i) { // std::cout << result.data[i] << "\n\n";