jess
/
fast_sparse_interpolation
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

MultiDimVector now provides a begin() end() interface with a corresponding iterator, also added some operations for vector arithmetic and reconstruction error check

This commit is contained in:
David Holzmüller 2019-04-11 18:50:10 +02:00
parent 7d02f0213a
commit 246b045406
3 changed files with 130 additions and 137 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

109
src/Interpolation.hpp
View File

@ -19,6 +19,7 @@ limitations under the License.
#include <boost/numeric/ublas/matrix.hpp>
#include <iostream>
#include <vector>
#include "Iterators.hpp"
namespace fsi {
@ -48,6 +49,47 @@ template <typename It>
class MultiDimVector {
It it;
class Iterator {
MultiDimVector<It> &v;
StepIterator<It> stepIt;
public:
Iterator(MultiDimVector<It> &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<size_t> 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<std::vector<double>> data;
@ -59,7 +101,10 @@ class MultiDimVector {
}
};
void swap(MultiDimVector &other) { data.swap(other.data); }
void swap(MultiDimVector<It> &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<It> &operator+=(MultiDimVector<It> 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<It> &operator-=(MultiDimVector<It> 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 <typename It>
MultiDimVector<It> operator+(MultiDimVector<It> first, MultiDimVector<It> const &second) {
first += second;
return first;
}
template <typename It>
MultiDimVector<It> operator-(MultiDimVector<It> first, MultiDimVector<It> const &second) {
first -= second;
return first;
}
template <typename It>
double squared_l2_norm(MultiDimVector<It> 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 <typename It>
void multiply_lower_triangular_inplace(std::vector<boost::numeric::ublas::matrix<double>> L,
MultiDimVector<It> &v) {
@ -91,7 +186,7 @@ void multiply_lower_triangular_inplace(std::vector<boost::numeric::ublas::matrix
double *data_pointer = &v.data[0][0];
size_t data_size = v.data[0].size();
while (not it.done()) {
while (it.valid()) {
size_t last_dim_count = it.lastDimensionCount();
double *offset_data_pointer = data_pointer + second_v_index;
for (size_t i = 0; i < last_dim_count; ++i) {
@ -99,8 +194,7 @@ void multiply_lower_triangular_inplace(std::vector<boost::numeric::ublas::matrix
for (size_t j = 0; j <= i; ++j) {
sum += Lk(i, j) * (*(offset_data_pointer + j));
}
w.data[i][indexes[i]] = sum;
++indexes[i];
w.data[i][indexes[i]++] = sum;
}
second_v_index += last_dim_count;
if (second_v_index >= data_size) {
@ -144,7 +238,7 @@ void multiply_upper_triangular_inplace(std::vector<boost::numeric::ublas::matrix
double *data_pointer = &v.data[0][0];
size_t data_size = v.data[0].size();
while (not it.done()) {
while (it.valid()) {
size_t last_dim_count = it.lastDimensionCount();
double *offset_data_pointer = data_pointer + second_v_index;
for (size_t i = 0; i < last_dim_count; ++i) {
@ -152,8 +246,7 @@ void multiply_upper_triangular_inplace(std::vector<boost::numeric::ublas::matrix
for (size_t j = i; j < last_dim_count; ++j) {
sum += Uk(i, j) * (*(offset_data_pointer + j));
}
w.data[i][indexes[i]] = sum;
++indexes[i];
w.data[i][indexes[i]++] = sum;
}
second_v_index += last_dim_count;
if (second_v_index >= data_size) {
@ -317,7 +410,7 @@ MultiDimVector<It> evaluateFunction(It it, Func f, X x) {
it.reset();
std::vector<double> 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));

148
src/Iterators.hpp
View File

@ -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 <size_t d>
class TemplateBoundedSumIterator {
size_t bound;
std::vector<size_t> 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<size_t>(d - 1, 0);
index_head_sum = 0;
is_done = false;
}
size_t dim() const { return d; }
std::vector<size_t> indexBounds() const { return std::vector<size_t>(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<d> cycle() const {
TemplateBoundedSumIterator<d> it = *this;
it.reset();
return it;
std::vector<size_t> index() const {
std::vector<size_t> result = jump_it.getIndexHead();
result.push_back(last_dim_value);
return result;
}
};
@ -142,11 +81,11 @@ class BoundedSumIterator {
size_t bound;
std::vector<size_t> 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<size_t> getIndexHead() const { return index_head; }
bool valid() const { return is_valid; }
void reset() {
index_head = std::vector<size_t>(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<size_t>(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<size_t> 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<size_t>(d - 1, 0);
index_sum = 0;
is_done = false;
}
};
}

10
test/main.cpp
View File

@ -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<d> it(bound);
fsi::BoundedSumIterator it(d, bound);
std::vector<MonomialFunctions> 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";
// }