Mercurial > hg > kwantix

changeset 14:e6d2cd3b6e77

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Intermediate whole domain evaluation work
author Jordi Gutiérrez Hermoso <jordigh@gmail.com>
date Sun, 03 Aug 2008 10:59:37 -0500
parents 0531194a431d
children 5144dd3c5468
files Makefile diff_op.cpp include/interpolator.hpp include/linalg.hpp interpolator.cpp linalg.cpp main.cpp
diffstat 7 files changed, 219 insertions(+), 186 deletions(-) [+]
line wrap: on
line diff
--- a/Makefile
+++ b/Makefile
@@ -1,6 +1,6 @@
 CPP = g++
-LINKING = -lgsl -lgslcblas  -pg
-CFLAGS = -pg
+LINKING = -lgsl -lgslcblas  
+CFLAGS = -g
 OPTIONS = -Wall -pedantic  -W -Werror -Wconversion -Wshadow \
 	  -Wpointer-arith -Wcast-qual -Wcast-align -Wwrite-strings  \
 	  -fshort-enums -fno-common -Wfatal-errors 
--- a/diff_op.cpp
+++ b/diff_op.cpp
@@ -47,7 +47,7 @@
     vector grad(dim);
     for(size_t i = 1; i <= dim; i++)
       grad(i) =f.d(p,i);
-    return grad*n/norm(n);
+    return grad%n/norm(n);
   }
 
   double Id_op::at(const realfunc &f, const point &p) const{
--- a/include/interpolator.hpp
+++ b/include/interpolator.hpp
@@ -52,7 +52,7 @@
     void interpolate(const map<point, double>& Xi); 
     void interpolate(shared_ptr<linear_BVP2> bvp);
     
-    /** Interpolate given values on the domain
+    /** \short Interpolate given values on the domain
      * 
      * The given values must match the values returned by other
      * interpolators on compatible domains.
@@ -60,31 +60,34 @@
     void interpolate(const interp_values& values); 
     //@}
     
-    /** @name Evaluations and derivatives
-     * 
-     * Full domain evaluations and derivatives are 
+
+    /** @name Pointwise evaluation and derivatives
      */
     //@{
     ///Point evaluation
     double operator()(const point& p) const;
     ///Point evaluation
     double at(const point& p) const;
-
+    /// Pointwise first derivative
+    double d(const point& p, size_t k) const;
+    /// Pointwise second derivatives
+    double d2(const point &p, size_t k1, size_t k2) const;
+    //@}
+    
+    /** @name Full domain evaluation and derivatives
+     */
+    //@{
     ///Full domain evaluation
     interp_values operator()() const; 
     ///Full domain evaluation
     interp_values at() const; 
-
-    /// Pointwise first derivative
-    double d(const point& p, size_t k) const;
     /// Full domain first erivative
     interp_values d(size_t k) const; 
-    /// Pointwise second derivatives
-    double d2(const point &p, size_t k1, size_t k2) const;
     /// Full domain second derivatives
     interp_values d2(size_t k1, size_t k2) const; 
     //@}
 
+
     /** @name Precomputation of RBFs
      *
      * The following functions will precompute upon request all the
@@ -94,9 +97,9 @@
      */
     //@{
     ///Precompute all the function evaluation of the RBFs.
-    void precompute_at();
+    void precompute_ev();
     ///Precompute all relevant first derivatives.
-    void precompute_d() ;
+    void precompute_d1() ;
     ///Precompute all relevant second derivatives.
     void precompute_d2();
     //@}
@@ -198,8 +201,8 @@
     /// Same, but storing the vectors.
     mutable map<std::vector<size_t>, linalg::vector >
     precomp_values_vec;
-    /// Invert phis\coeffs and return result vector in a map.
-    map<point, double> precompute_values(const matrix& phis) const;
+    /// Convert vector of values to a map
+    map<point, double> valsvec2map(const linalg::vector& values) const;
   };
 
   /// Interpolated values manipulated by the interpolator class.
@@ -217,6 +220,8 @@
     interp_values operator-(const interp_values& w) const;
     interp_values operator*(const interp_values& w) const;
     interp_values operator/(const interp_values& w) const;
+    interp_values operator*(const double a) const;
+    interp_values operator/(const double a) const;
     //@}
   private:
     
@@ -234,14 +239,13 @@
     void different_rbfs(int line, string file) const;
 
     template<typename RBF> friend class interpolator;
+    friend interp_values operator*(const interp_values& v, const
+				   interp_values& w); 
   };
 
   ///For comfortable syntax
   template <typename RBF>
-  interpolator<RBF> operator*(double a, const interpolator<RBF>& u)
-  {
-    return u*a;
-  }
+  interpolator<RBF> operator*(double a, const interpolator<RBF>& u);
 }
 
 #endif
--- a/include/linalg.hpp
+++ b/include/linalg.hpp
@@ -261,7 +261,9 @@
 
     //Arithmetic operations:
     ///Scale the vector.
-    vector operator*(const double a)  const; 
+    vector operator*(const double a)  const;
+    ///Scale the vector.
+    vector operator/(const double a)  const;
     ///Vector addition.
     vector operator+(const vector& w) const;
     ///Vector subtration.
--- a/interpolator.cpp
+++ b/interpolator.cpp
@@ -159,6 +159,9 @@
 	map<std::vector<size_t>, matrix>::const_iterator I;
       for(I = precomp_rbfs.begin(); I != precomp_rbfs.end(); I++)
 	precomp_values_vec[I -> first] = (I -> second)*coeffs;
+      
+      //Pointwise map evaluations we will not compute again.
+      precomp_values.clear();
     }
   }
 
@@ -239,37 +242,12 @@
     return seed;
   }
   
-
-  //***** Evaluations, derivatives *************
-
-  template<typename RBF>
-  double interpolator<RBF>::operator()(const point& p) const
-  {
-    return at(p);
-  }
-
-  template<typename RBF>
-  map<point, double> interpolator<RBF>::precompute_values(const matrix& phis) 
-    const
-  {
-    linalg::vector values = phis*coeffs;
-    map<point, double> out;
-    {
-      set<point>::const_iterator I;
-      size_t i;
-      for(I = thebvp -> get_domain() -> get_interior().begin(), i = 1;
-	  i <= n; I++, i++)
-	out[*I] = values(i);
-      
-      for(I = thebvp -> get_domain() -> get_boundary().begin(), i = n+1;
-	  i <= n+m; I++, i++)
-	out[*I] = values(i);
-    }
-    return out;
-  }
+  // ******************** Evaluations, derivatives ********************
+ 
+  //  ---------- Precomputation ---------------------------
   
   template<typename RBF>
-  void interpolator<RBF>::precompute_at()
+  void interpolator<RBF>::precompute_ev()
   {
     if(!initted){
       not_initted(__LINE__, __FILE__);
@@ -293,65 +271,14 @@
 	  phis(i,j+1) = rbfs[j].at(*I);
       
       precomp_rbfs[alpha] = phis;
-      precomp_values[alpha] = precompute_values(phis);
-    }
-
-    ev_precomp = true;
-  }
-
-  template<typename RBF>
-  double interpolator<RBF>::at(const point& p) const
-  {
-    if(!initted){
-      not_initted(__LINE__, __FILE__);
+      computecoeffs();
+      
+      ev_precomp = true;
     }
-    std::vector<size_t> alpha; //empty vector
-
-    //Are we evaluating at a precomputed point in the domain?
-    if(ev_precomp and 
-       precomp_values[alpha].find(p) != precomp_values[alpha].end())
-      return precomp_values[alpha][p];
-     
-    
-    //Else, must compute the value
-    double result = 0;
-    for(size_t i = 1; i <= coeffs.size(); i++)
-      result += coeffs(i)*rbfs[i-1].at(p);
-        
-    return result;
   }
 
   template<typename RBF>
-  void interpolator<RBF>::not_precomputed(int line, string file) const
-  {
-    badArgument exc;
-    exc.reason = "Cannot do whole domain evaluations without "
-      "precomputed RBFs.";
-    exc.line = line;
-    exc.file = file;
-    throw exc;
-  }
-
-  template<typename RBF>
-  interp_values interpolator<RBF>::operator()() const
-  {
-    return at();
-  }
-  
-  template<typename RBF>
-  interp_values interpolator<RBF>::at() const
-  {
-    if(!ev_precomp)
-      not_precomputed(__LINE__, __FILE__);
-     
-    std::vector<size_t> alpha;
-    return interp_values(precomp_values_vec[alpha],rbfs_hash);
-  }
-  
-  
-
-  template<typename RBF>
-  void interpolator<RBF>::precompute_d()
+  void interpolator<RBF>::precompute_d1()
   {
     if(!initted){
       not_initted(__LINE__, __FILE__);
@@ -380,13 +307,74 @@
 	}
 
 	precomp_rbfs[alpha] = phis;
-	precomp_values[alpha] = precompute_values(phis);
+	computecoeffs();
       }
     }
     d1_precomp = true;
   }
 
   template<typename RBF>
+  void interpolator<RBF>::precompute_d2()
+  {
+    if(!initted)
+      not_initted(__LINE__, __FILE__);
+    if(!d2_precomp){
+      size_t dim = thebvp -> get_domain() -> get_dimension();
+      //precompute all second derivatives
+      for(size_t k1 = 1; k1 <= dim; k1++){
+	for(size_t k2 = k1; k2 <= dim; k2++){
+	  std::vector<size_t> alpha(k1>k2?k1:k2); alpha[k1-1]++; alpha[k2-1]++;
+	    matrix phis(n+m,n+m);
+	    set<point>::iterator I;
+	    size_t i;
+	    for(I = thebvp -> get_domain () -> get_interior().begin(), i = 1; 
+		i <= n; i++, I++)
+	      for(size_t j = 0; j < n+m; j++)
+		phis(i,j+1) = rbfs[j].d2(*I,k1,k2);
+      
+	    for(I = thebvp -> get_domain() -> get_boundary().begin(), i=n+1; 
+		i <= n+m; i++, I++)
+	      for(size_t j = 0; j < n+m; j++)
+		phis(i,j+1) = rbfs[j].d2(*I,k1,k2);
+      
+	    precomp_rbfs[alpha] = phis;
+	    computecoeffs();	  
+	}
+      }
+    }
+    d2_precomp = true;
+  }
+
+  // ----------------- Point evaluation ---------------------------
+
+  template<typename RBF>
+  double interpolator<RBF>::operator()(const point& p) const
+  {
+    return at(p);
+  }
+
+  template<typename RBF>
+  double interpolator<RBF>::at(const point& p) const
+  {
+    if(!initted){
+      not_initted(__LINE__, __FILE__);
+    }
+    std::vector<size_t> alpha; //empty vector
+
+    //Are we evaluating at a precomputed point in the domain?
+    if(ev_precomp and 
+       precomp_values[alpha].find(p) != precomp_values[alpha].end())
+      return precomp_values[alpha][p];
+     
+    //Else, must compute the value
+    double result = 0;
+    for(size_t i = 1; i <= coeffs.size(); i++)
+      result += coeffs(i)*rbfs[i-1].at(p);
+        
+    return result;
+  }
+
+  template<typename RBF>
   double interpolator<RBF>::d(const point& p, size_t k) const{
     if(!initted){
       not_initted(__LINE__, __FILE__);
@@ -407,49 +395,6 @@
   }
 
   template<typename RBF>
-  interp_values interpolator<RBF>::d(size_t k) const
-  {
-    if(!d1_precomp)
-      not_precomputed(__LINE__, __FILE__);
-    
-    std::vector<size_t> alpha(k); alpha[k-1]++;
-    return interp_values(precomp_values_vec[alpha],rbfs_hash);
-  }
-
-  template<typename RBF>
-  void interpolator<RBF>::precompute_d2()
-  {
-    if(!initted)
-      not_initted(__LINE__, __FILE__);
-    if(!d2_precomp){
-      size_t dim = thebvp -> get_domain() -> get_dimension();
-      //precompute all second derivatives of t
-      for(size_t k1 = 1; k1 <= dim; k1++){
-	for(size_t k2 = k1; k2 <= dim; k2++){
-	  std::vector<size_t> alpha(k1>k2?k1:k2); alpha[k1-1]++; alpha[k2-1]++;
-	    matrix phis(n+m,n+m);
-	    set<point>::iterator I;
-	    size_t i;
-	    for(I = thebvp -> get_domain () -> get_interior().begin(), i = 1; 
-		i <= n; i++, I++)
-	      for(size_t j = 0; j < n+m; j++)
-		phis(i,j+1) = rbfs[j].d2(*I,k1,k2);
-      
-	    for(I = thebvp -> get_domain() -> get_boundary().begin(), i=n+1; 
-		i <= n+m; i++, I++)
-	      for(size_t j = 0; j < n+m; j++)
-		phis(i,j+1) = rbfs[j].d2(*I,k1,k2);
-      
-	    precomp_rbfs[alpha] = phis;
-	    precomp_values[alpha] = precompute_values(phis);
-	  
-	}
-      }
-    }
-    d2_precomp = true;
-  }
-
-  template<typename RBF>
   double interpolator<RBF>::d2(const point& p, size_t k1, size_t k2) const{
     if(!initted){
       not_initted(__LINE__, __FILE__);
@@ -469,6 +414,45 @@
     return result;
   }
 
+  // ------------------- Whole domain evaluations ----------------
+
+  template<typename RBF>
+  void interpolator<RBF>::not_precomputed(int line, string file) const
+  {
+    badArgument exc;
+    exc.reason = "Cannot do whole domain evaluations without "
+      "precomputed RBFs.";
+    exc.line = line;
+    exc.file = file;
+    throw exc;
+  }
+
+  template<typename RBF>
+  interp_values interpolator<RBF>::operator()() const
+  {
+    return at();
+  }
+  
+  template<typename RBF>
+  interp_values interpolator<RBF>::at() const
+  {
+    if(!ev_precomp)
+      not_precomputed(__LINE__, __FILE__);
+     
+    std::vector<size_t> alpha;
+    return interp_values(precomp_values_vec[alpha],rbfs_hash);
+  }
+   
+  template<typename RBF>
+  interp_values interpolator<RBF>::d(size_t k) const
+  {
+    if(!d1_precomp)
+      not_precomputed(__LINE__, __FILE__);
+    
+    std::vector<size_t> alpha(k); alpha[k-1]++;
+    return interp_values(precomp_values_vec[alpha],rbfs_hash);
+  }
+
   template<typename RBF>
   interp_values interpolator<RBF>::d2(size_t k1, size_t k2) const
   {
@@ -479,7 +463,7 @@
     return interp_values(precomp_values_vec[alpha],rbfs_hash);
   }
 
-  //********** linear operations ***************
+  //***************** linear operations ***************************
 
   template<typename RBF>
   interpolator<RBF> interpolator<RBF>::
@@ -535,7 +519,7 @@
     return u;
   }
   
-  //********** data manipulation ***************
+  //******************* data manipulation **************************
 
   template<typename RBF>
   void interpolator<RBF>::set_f(const realfunc& f){
@@ -611,14 +595,36 @@
     coeffs = M.inv(rhs);
     
     //Precompute the values for the derivatives and evaluations that
-    //have already  been computed.
+    //have already been computed.
     {
       typename map<std::vector<size_t>, matrix>::const_iterator I;
       
-      for(I = precomp_rbfs.begin(); I != precomp_rbfs.end(); I++)
-	precomp_values[I -> first] = precompute_values(I->second);
+      for(I = precomp_rbfs.begin(); I != precomp_rbfs.end(); I++){
+	linalg::vector vals = (I->second)*coeffs;
+	precomp_values_vec[I -> first] = vals;
+	precomp_values[I -> first] = valsvec2map(vals);
+      }
     }
   }
+
+  template<typename RBF>
+  map<point, double> interpolator<RBF>::
+  valsvec2map(const linalg::vector& values) const
+  {
+    map<point, double> out;
+    {
+      set<point>::const_iterator I;
+      size_t i;
+      for(I = thebvp -> get_domain() -> get_interior().begin(), i = 1;
+	  i <= n; I++, i++)
+	out[*I] = values(i);
+      
+      for(I = thebvp -> get_domain() -> get_boundary().begin(), i = n+1;
+	  i <= n+m; I++, i++)
+	out[*I] = values(i);
+    }
+    return out;
+  }
   
   // ************************ interp_values stuff *****************
   void interp_values::different_rbfs(int line, string file) const
@@ -666,6 +672,35 @@
     return interp_values(v / w.v,rbfs_hash);
   }
 
+  interp_values
+  interp_values::operator*(const double a) const
+  {
+    return interp_values(v*a,rbfs_hash);
+  }
+  
+  interp_values
+  interp_values::operator/(const double a) const
+  {
+    return interp_values(v/a, rbfs_hash);
+  }
+
+  //*************** Non-member global friend functions ******************
+
+  template <typename RBF>
+  interpolator<RBF> operator*(double a, const interpolator<RBF>& u)
+  {
+    return u*a;
+  }  
+
+  interp_values operator*(const interp_values& v, const interp_values&
+  w)
+  {
+    if(v.rbfs_hash != w.rbfs_hash)
+      v.different_rbfs(__LINE__, __FILE__);
+    
+    return interp_values(v.v / w.v,v.rbfs_hash);
+  } 
+
   //Instantiations
   using namespace rbf;
   template class interpolator<piecewise_polynomial>;
--- a/linalg.cpp
+++ b/linalg.cpp
@@ -492,6 +492,12 @@
     gsl_vector_scale(v.x, a);
     return v;
   }
+  
+  vector vector::operator/(const double a) const{
+    vector v = *this;
+    gsl_vector_scale(v.x, 1.0/a);
+    return v;
+  }
 
   vector vector::operator+(const vector& w) const{
     if(x -> size != w.x->size){
--- a/main.cpp
+++ b/main.cpp
@@ -36,7 +36,7 @@
 const double g = 9.8; // m/s^2
 
 template<typename RBF>
-class Fgen : public realfunc{
+class Fgen {
 public:
   void set_interps(const interpolator<RBF>& u0, 
 		   const interpolator<RBF>& v0,
@@ -59,7 +59,7 @@
   using Fgen<RBF>::h;
   using Fgen<RBF>::dt;
 public:
-  double at(const point& p) const;
+  bvp::interp_values at() const;
 };
 
 template<typename RBF>
@@ -69,7 +69,7 @@
   using Fgen<RBF>::h;
   using Fgen<RBF>::dt;
 public:
-  double at(const point& p) const;
+  bvp::interp_values at() const;
 };
 
 template<typename RBF>
@@ -79,9 +79,10 @@
   using Fgen<RBF>::h;
   using Fgen<RBF>::dt;
 public:
-  double at(const point& p) const;
+  bvp::interp_values at() const;
 };
 
+//FIXME: This iteration must work on interpolators without boundary
 class iter_neumann : public bdry_diff_op{
 public:
   //If this is u's, set 1, for v, set 2.
@@ -164,15 +165,17 @@
       h1, h0(Omega, h_init);
 
     //Precompute some RBFs
-    u0.precompute_at(); v0.precompute_at(), h0.precompute_at();
-    u0.precompute_d(); v0.precompute_d(), h0.precompute_d();
+    u0.precompute_ev(); v0.precompute_ev(), h0.precompute_ev();
+    u0.precompute_d1(); v0.precompute_d1(), h0.precompute_d1();
 
 
     u1 = u0; v1 = v0; h1 = h0;
 
     //Intermediate interpolators for RK4
     std::vector<interpolator<RBF_TYPE> >
-      k1(3,h0), k2(3,h0), k3(3,h0), k4(3,h0);
+      // FIXME: make the copy ctor for interps make different copies
+      // of the bvp.
+      k1(3,h0), k2(3,h0), k3(3,h0), k4(3,h0); 
 
     Fu<RBF_TYPE> f_u;
     Fv<RBF_TYPE> f_v;
@@ -182,8 +185,6 @@
     f_v.set_dt(dt);
     f_h.set_dt(dt);
 
-
-
     //main loop
     size_t maxiter = 1000;
     for(size_t i = 1; i <= maxiter; i++){ 
@@ -199,9 +200,6 @@
       f_v.set_interps(u0,v0,h0);
       f_h.set_interps(u0,v0,h0);
       
-      k1[0].set_f(f_u);
-      k1[1].set_f(f_v);
-      k1[2].set_f(f_h);
       cout << "k1 done!" << endl;
 
       bdry_iter(k1[0],k1[1],Omega);
@@ -211,10 +209,6 @@
       f_v.set_interps(u0+(k1[0]/2), v0+(k1[1]/2), h0+(k1[2]/2));
       f_h.set_interps(u0+(k1[0]/2), v0+(k1[1]/2), h0+(k1[2]/2));
  
-      k2[0].set_f(f_u);
-      k2[1].set_f(f_v);
-      k2[2].set_f(f_h);
-
       bdry_iter(k2[0],k2[1],Omega);
 
       cout << "k2 done!" << endl;
@@ -224,10 +218,6 @@
       f_v.set_interps(u0+(k2[0]/2), v0+(k2[1]/2), h0+(k2[2]/2));
       f_h.set_interps(u0+(k2[0]/2), v0+(k2[1]/2), h0+(k2[2]/2));
       
-      k3[0].set_f(f_u);
-      k3[1].set_f(f_v);
-      k3[2].set_f(f_h);
-
       bdry_iter(k3[0],k3[1],Omega);
       
       cout << "k3 done!" << endl;
@@ -237,10 +227,6 @@
       f_v.set_interps(u0+k3[0], v0+k3[1], h0+k3[2]);
       f_h.set_interps(u0+k3[0], v0+k3[1], h0+k3[2]);
       
-      k4[0].set_f(f_u);
-      k4[1].set_f(f_v);
-      k4[2].set_f(f_h);
-
       bdry_iter(k4[0],k4[1],Omega);
 
       cout << "k4 done!" << endl;
@@ -262,7 +248,7 @@
   }
   catch(error exc){
     utils::show_exception(exc);
-    return 1;
+    return ;
   }
 
 }
@@ -304,21 +290,21 @@
 
 
 template<typename RBF>
-double Fu<RBF>::at(const point& p) const
+interp_values Fu<RBF>::at() const
 {
-  return   dt*(u(p)*u.d(p,1) +  v(p)*u.d(p,2) +  g*h.d(p,1));
+  return   dt*(u()*u.d(1) +  v()*u.d(2) +  g*h.d(1));
 }
 
 template<typename RBF>
-double Fv<RBF>::at(const point& p) const
+interp_values Fv<RBF>::at() const
 {
-  return   dt*(u(p)*v.d(p,1) +  v(p)*v.d(p,2) +  g*h.d(p,2));
+  return   dt*(u()*v.d(1) +  v()*v.d(2) +  g*h.d(2));
 }
 
 template<typename RBF>
-double Fh<RBF>::at(const point& p) const
+interp_values Fh<RBF>::at() const
 {
-  return   dt*(u(p)*h.d(p,1) + h(p)*u.d(p,1) +   v(p)*h.d(p,2) + h(p)*v.d(p,2));
+  return   dt*(u()*h.d(1) + h()*u.d(1) +   v()*h.d(2) + h()*v.d(2));
 }
 
 iter_neumann::iter_neumann(size_t u_or_v_in)