From fe4e4276819293c23ba887d2d956fcba22032f1c Mon Sep 17 00:00:00 2001 From: Scott Thornton Date: Tue, 26 Apr 2016 18:24:58 -0400 Subject: [PATCH] got the tree and key sorted out ... things better now --- Matrix.h | 17 +++++++++++--- function1d.h | 59 ++++++++++++++++++++++++++++++++++++++++++------ test_function.cc | 4 +++- tree.h | 2 +- twoscale.h | 25 ++++++++++++++++---- 5 files changed, 90 insertions(+), 17 deletions(-) diff --git a/Matrix.h b/Matrix.h index 01c8a9e..ceaefea 100644 --- a/Matrix.h +++ b/Matrix.h @@ -226,18 +226,24 @@ class VectorT { return v; } - T* begin() { + T* begin() const { return &_p[0]; } - T* end() { - return &_p[size()-1]; + T* end() const { + return &_p[size()]; } }; typedef VectorT real_vector; typedef VectorT > complex_vector; +void print(const real_vector& v) { + printf("["); + for (auto& t : v) printf("%10.5e ", t); + printf("]\n"); +} + template class MatrixT { private: @@ -567,6 +573,11 @@ double inner(const MatrixT& t1, const MatrixT& t2) { return t1.inner(t2); } +template +double inner(const VectorT& t1, const VectorT& t2) { + return t1.inner(t2); +} + double norm2(const real_matrix& t) { return t.norm2(); } diff --git a/function1d.h b/function1d.h index 83f554c..82fbae9 100644 --- a/function1d.h +++ b/function1d.h @@ -7,13 +7,13 @@ using Vector = real_vector; double normf(Vector v) { auto s = 0.0; - for (auto i = 0; i < v.size(); i++) s+=v[i]; - return s; + for (auto t : v) s+=t*t; + return std::sqrt(s); } class Function1D { private: - bool debug = true; + bool debug = false; int k = 8; double thresh = 1e-6; int maxlevel = 30; @@ -95,7 +95,20 @@ class Function1D { s[i+k] = s1[i]; } Vector d = hg*s; - if (debug) printf(" d = hg*s\n"); + if (debug) { + printf("s0: "); + print(s0); + printf("s1: "); + print(s1); + printf("s: "); + print(s); + printf("d: "); + print(d); + printf("d.slice(k,2*k-1): "); + print(Vector(d.slice(k,2*k-1))); + printf("d slice norm: %15.8e\n", normf(Vector(d.slice(k,2*k-1)))); + } + if (debug) printf("\n d = hg*s\n"); if (normf(Vector(d.slice(k,2*k-1))) < thresh || n >= maxlevel-1) { tree[Key(n+1,2*l)] = s0; if (debug) printf("set n+1 2*l coeff (%d %d)\n",n+1,2*l); @@ -107,19 +120,51 @@ class Function1D { } } + double operator()(double x) { + return eval(x, 0, 0); + } + + double eval(double x, int n, int l) { + assert(n < maxlevel); + auto treep = tree.find(Key(n,l)); + if (treep != tree.end()) { + auto p = ScalingFunction::instance()->phi(x, k); + auto t = inner(treep->second,p)*std::sqrt(std::pow(2.0,n)); + return t; + } else { + auto n2 = n + 1; + auto l2 = 2*l; + auto x2 = 2*x; + if (x2 >= 1.0) { + l2 = l2 + 1; + x2 = x2 - 1; + } + return eval(x2, n2, l2); + } + } + void print_coeffs(int n, int l) { auto s = tree[Key(n,l)]; printf("[%d, %d] (", n, l); for (auto v : s) { printf("%8.4f ", v); } - printf(")\n"); + printf(") %15.8e\n",normf(s)); } void print_tree() { for (auto c : tree) { - Key k = c.first; - printf("[%d %d]\n", k.n, k.l); + auto k = c.first; + auto s = c.second; + printf("[%d %d] %15.8e\n", k.n, k.l, normf(s)); } } + + // void summarize() { + // printf("sum coeffs:\n"); + // for (auto c : tree) { + // auto key = c.first; + // auto s = c.second; + // } + // } }; diff --git a/test_function.cc b/test_function.cc index 3520559..46b5d9a 100644 --- a/test_function.cc +++ b/test_function.cc @@ -10,7 +10,9 @@ double func(double x) { int main(int argc, char** argv) { - Function1D f(func, 8, 1e-6, 300); + Function1D f(func, 8, 1e-10, 30, 2); f.print_tree(); + auto x = 0.23111; + printf("x: %15.8e f: %15.8e func: %15.8e error: %15.8e\n", x, f(x), func(x), std::abs(f(x)-func(x))); return 0; } diff --git a/tree.h b/tree.h index 85a226c..93e2467 100644 --- a/tree.h +++ b/tree.h @@ -9,7 +9,7 @@ struct Key { Key(int n, int l) : n(n), l(l) {} bool operator< (const Key &k) const { - return ((1< #include #include #include "Matrix.h" +using Vector = real_vector; + //Return the twoscale coefficients for the multiwavelets of order k. // //Note that a cached value is returned ... if you want to modify it @@ -9,6 +12,7 @@ class TwoScaleCoeffs { private: std::vector coeffs; + // static std::shared_ptr _instance; static TwoScaleCoeffs* _instance; TwoScaleCoeffs() { std::ifstream fil("tscoeffs"); @@ -29,20 +33,25 @@ class TwoScaleCoeffs { } } } - ~TwoScaleCoeffs() {} public: + // static std::shared_ptr instance() { + // if (!_instance) _instance = std::make_shared(); + // return _instance; + // } static TwoScaleCoeffs* instance() { if (!_instance) _instance = new TwoScaleCoeffs(); return _instance; } real_matrix hg(int k) {return coeffs[k-1];} + ~TwoScaleCoeffs() {} }; +// std::shared_ptr TwoScaleCoeffs::_instance = 0; TwoScaleCoeffs* TwoScaleCoeffs::_instance = 0; // Evaluate the Legendre polynomials up to the given order at x // defined on [-1,1]. -std::vector legendre(double x,int order) { - std::vector p(order+1); +Vector legendre(double x,int order) { + Vector p(order+1); p[0] = 1.0; if (order == 0) return p; p[1] = x; @@ -55,14 +64,18 @@ std::vector legendre(double x,int order) { class ScalingFunction { private: double norms[100]; + // static std::shared_ptr _instance; static ScalingFunction* _instance; ScalingFunction() { for (int i = 0; i < 100; i++) { norms[i] = std::sqrt(2*i+1); } } - ~ScalingFunction() {} public: + // static std::shared_ptr instance() { + // if (!_instance) _instance = std::make_shared(); + // return _instance; + // } static ScalingFunction* instance() { if (!_instance) _instance = new ScalingFunction(); return _instance; @@ -73,14 +86,16 @@ class ScalingFunction { // In addition to forming an orthonormal basis on [0,1] we have // phi_j(1/2-x) = (-1)^j phi_j(1/2+x) // (the wavelets are similar with phase (-1)^(j+k)). - std::vector phi(double x,int k) { + Vector phi(double x,int k) { auto order = k-1; auto p = legendre(2.*x-1, order); for (auto n = 0; n < k; n++) p[n] = p[n]*norms[n]; return p; } + ~ScalingFunction() {} }; +// std::shared_ptr ScalingFunction::_instance = 0; ScalingFunction* ScalingFunction::_instance = 0;