ok inserted empty vector in the interior nodes in the scaling functio…

…n basis

function1d.h

@@ -35,6 +35,7 @@ class Function1D {
   real_matrix quad_phiw;
 
 public:
+  // I know that friends are evil, but I decided to have them anyway
   friend Function1D operator*(const double& s, const Function1D& f);
   friend Function1D operator*(const Function1D& f, const double& s);
   friend Function1D compress(const Function1D& f);
@@ -55,6 +56,12 @@ class Function1D {
     init_quadrature(k);
     int ntrans = std::pow(2, initiallevel);
     for (auto l = 0; l < ntrans; l++) refine(f, initiallevel, l);
+    for (auto n = 0; n <= initiallevel; n++) {
+      ntrans = std::pow(2,n);
+      for (auto l = 0; l < ntrans; l++) {
+        stree[Key(n,l)] = Vector();
+      }
+    }
   }
 
   ~Function1D() {}
@@ -112,6 +119,7 @@ class Function1D {
       stree[Key(n+1,2*l)] = s0;
       stree[Key(n+1,2*l+1)] = s1;
     } else {
+      stree[Key(n,l)] = Vector();
       refine(f, n+1, 2*l);
       refine(f, n+1, 2*l+1);
     }
@@ -135,6 +143,7 @@ class Function1D {
       }
       reconstruct_spawn(stree_r, s0, n+1, 2*l);
       reconstruct_spawn(stree_r, s1, n+1, 2*l+1);
+      stree_r[Key(n,l)] = Vector();
     } else {
       stree_r[Key(n,l)] = ss;
     }
@@ -143,8 +152,18 @@ class Function1D {
   Vector compress_spawn(CoeffTree& dtree_r, int n, int l) const {
     auto s0p = stree.find(Key(n+1,2*l));
     auto s1p = stree.find(Key(n+1,2*l+1));
-    Vector s0 = (s0p == stree.end()) ? compress_spawn(dtree_r, n+1, 2*l)   : s0p->second;
-    Vector s1 = (s1p == stree.end()) ? compress_spawn(dtree_r, n+1, 2*l+1) : s1p->second;
+    Vector s0;
+    Vector s1;
+    if (s0p != stree.end()) {
+      s0 = ((s0p->second).size() == 0) ? compress_spawn(dtree_r, n+1, 2*l) : s0p->second;
+    } else {
+      s0 = compress_spawn(dtree_r, n+1, 2*l);
+    }
+    if (s1p != stree.end()) {
+      s1 = ((s1p->second).size() == 0) ? compress_spawn(dtree_r, n+1, 2*l+1) : s1p->second;
+    } else {
+      s1 = compress_spawn(dtree_r, n+1, 2*l+1);
+    }
     Vector s(2*k);
     for (auto i = 0; i < k; i++) {
       s[i]   = s0[i];
@@ -222,7 +241,7 @@ class Function1D {
   double eval(double x, int n, int l) {
     assert(n < maxlevel);
     auto treep = stree.find(Key(n,l));
-    if (treep != stree.end()) {
+    if (treep != stree.end() && (treep->second).size()) {
       auto p = ScalingFunction::instance()->phi(x, k);
       auto t = inner(treep->second,p)*std::sqrt(std::pow(2.0,n));
       return t;
@@ -276,10 +295,12 @@ Function1D operator*(const double& s, const Function1D& f) {
   r.form = f.form;
   for (auto cp : f.stree) {
     auto key = cp.first;
-    auto c = cp.second; 
-    auto c2 = copy(c);
-    c2.scale(s);
-    r.stree[key] = c2;
+    if ((cp.second).size()) {
+      auto c = cp.second; 
+      auto c2 = copy(c);
+      c2.scale(s);
+      r.stree[key] = c2;
+    }
   }
   for (auto cp : f.dtree) {
     auto key = cp.first;

test_function.cc

@@ -1,8 +1,9 @@
 #include "function1d.h"
 
 const double PI = 3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825;
-const int k = 16;
+const int k = 4;
 const double thresh = 1e-8;
+const int initiallevel = 2;
 
 double func1(double x) {
   auto c = 1.0;
@@ -20,37 +21,43 @@ double sum_f12(double x) {
   return func1(x) + func2(x);
 }
 
+void test_function_point() {
+   Function1D f(func1, k, thresh, 30, initiallevel);
+   f.print_tree();
+   auto x = 0.23111;
+   printf("x: %15.8e f: %15.8e func: %15.8e error: %15.8e\n", x, f(x), func1(x), std::abs(f(x)-func1(x)));
+}
+
+void test_function_compress() {
+   Function1D f(func1, k, thresh, 30, initiallevel);
+   printf("f : \n");
+   f.print_tree();
+   Function1D g = reconstruct(compress(f));
+   //Function1D g = compress(f);
+   printf("\ng : \n");
+   g.print_tree();
+}
+
 void test_add() {
-  auto fr = Function1D(func1, k, thresh, 30, 2);
+  auto fr = Function1D(func1, k, thresh, 30, initiallevel);
   auto fc = compress(fr);
   printf("\nfunction fc:\n");
-  fc.print_tree();
-  auto gr = Function1D(func2, k, thresh, 30, 2);
+  fr.print_tree();
+  auto gr = Function1D(func2, k, thresh, 30, initiallevel);
   auto gc = compress(gr);
   printf("\nfunction gc:\n");
-  gc.print_tree();
+  gr.print_tree();
   auto fgc = fc + gc;
   printf("\nfunction f+g(c):\n");
   fgc.print_tree();
-  auto f12r = Function1D(sum_f12, k, thresh, 30, 2);
+  auto f12r = Function1D(sum_f12, k, thresh, 30, initiallevel);
   auto f12c = compress(f12r);
   printf("\nfunction f12c:\n");
   f12c.print_tree();
 }
 
 int main(int argc, char** argv) {
-
-//   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)));
-//   auto g = compress(f);
-//   auto g2 = g + g;
-//   g.print_tree();
-//   g2.print_tree();
-//   auto h = 2.0*reconstruct(g);
-//   h.print_tree();
-
+  //test_function_compress();
   test_add();
   return 0;
 }