got the tree and key sorted out ... things better now

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<double> real_vector;
 typedef VectorT<std::complex<double> > complex_vector;
 
+void print(const real_vector& v) {
+  printf("[");
+  for (auto& t : v) printf("%10.5e  ", t);
+  printf("]\n");
+}
+
 template <typename T>
 class MatrixT {
 private:
@@ -567,6 +573,11 @@ double inner(const MatrixT<Q>& t1, const MatrixT<Q>& t2) {
   return t1.inner(t2);
 }
 
+template <typename Q>
+double inner(const VectorT<Q>& t1, const VectorT<Q>& t2) {
+  return t1.inner(t2);
+}
+
 double norm2(const real_matrix& t) {
   return t.norm2();
 }

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;
+  //   }
+  // }
 };

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;
 }

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<<n)+l < (k.n<<n)+k.l);
+    return ((1<<n)+l < (1<<k.n)+k.l);
   }
 };
 

twoscale.h

@@ -1,14 +1,18 @@
+#include <memory>
 #include <string>
 #include <fstream>
 #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
 //take a copy first
 class TwoScaleCoeffs {
 private:
   std::vector<real_matrix> coeffs;
+  // static std::shared_ptr<TwoScaleCoeffs> _instance;
   static TwoScaleCoeffs* _instance;
   TwoScaleCoeffs() {
     std::ifstream fil("tscoeffs");
@@ -29,20 +33,25 @@ class TwoScaleCoeffs {
       }
     }
   }
-  ~TwoScaleCoeffs() {}
 public:
+  // static std::shared_ptr<TwoScaleCoeffs> instance() {
+  //   if (!_instance) _instance = std::make_shared<TwoScaleCoeffs>();
+  //   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> TwoScaleCoeffs::_instance = 0;
 TwoScaleCoeffs* TwoScaleCoeffs::_instance = 0;
 
 // Evaluate the Legendre polynomials up to the given order at x
 // defined on [-1,1].
-std::vector<double> legendre(double x,int order) {
-  std::vector<double> 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<double> legendre(double x,int order) {
 class ScalingFunction {
 private:
   double norms[100];
+  // static std::shared_ptr<ScalingFunction> _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<ScalingFunction> instance() {
+  //   if (!_instance) _instance = std::make_shared<ScalingFunction>();
+  //   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<double> 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> ScalingFunction::_instance = 0;
 ScalingFunction* ScalingFunction::_instance = 0;