(JEC) 14/10/15 improve the Local Integration strate (quadinteg.h) to compute...

(JEC) 14/10/15 improve the Local Integration strate (quadinteg.h) to compute the laguerre2bessel coeffs.

src/lagsht_testsuite.cc

@@ -124,7 +124,7 @@ void MultiLaguerreTransformSdtAlone() {
     fnOrig[i] = complex<r_8>(rv,iv);
   }
 #if DEBUG >=2 
-  std::copy(fnOrig.begin(), fnOrig.end(), std::ostream_iterator<char>(std::cout, " "));
+  std::copy(fnOrig.begin(), fnOrig.end(), std::ostream_iterator< complex<r_8> >(std::cout, " "));
 #endif
     
   LaguerreTransform trans(N,R);
@@ -136,7 +136,7 @@ void MultiLaguerreTransformSdtAlone() {
   trans.MultiSynthesis(fnOrig,fi,mapSize);
     
 #if DEBUG >=2 
-  std::copy(fi.begin(), fi.end(), std::ostream_iterator<char>(std::cout, " "));
+  std::copy(fi.begin(), fi.end(), std::ostream_iterator< complex<r_8> >(std::cout, " "));
 #endif
 
   vector< complex<r_8> > fn(Ntot);
@@ -221,7 +221,7 @@ void MultiSphericalLaguerreTransform() {
   }//end loop on shell
 
 #if DEBUG >=2 
-  std::copy(flmnOrig.begin(), flmnOrig.end(), std::ostream_iterator<char>(std::cout, " "));
+  std::copy(flmnOrig.begin(), flmnOrig.end(), std::ostream_iterator< complex<r_8> >(std::cout, " "));
 #endif
       
   
@@ -342,7 +342,7 @@ void TestLaguerre2Bessel() {
   }//end loop on shell
 
 #if DEBUG >=2 
-  std::copy(flmnOrig.begin(), flmnOrig.end(), std::ostream_iterator<char>(std::cout, " "));
+  std::copy(flmnOrig.begin(), flmnOrig.end(), std::ostream_iterator< complex<r_8> >(std::cout, " "));
 #endif
         
   cout << "tstack 1 End" <<endl;
@@ -558,7 +558,6 @@ int main(int narg, char *arg[]) {
 
   param.Lmax = Lmax;
   param.N    = N;
-  if(Pmax<N){Pmax = 32*N; cout << "(Warning): Pmax<Nmax => modify the value to "<< Pmax << endl;}
   param.Pmax = Pmax;
   param.R    = R;
   param.alpha = 2;
@@ -617,6 +616,7 @@ int main(int narg, char *arg[]) {
 
     case 5:
       {
+	if(Pmax<N){Pmax = 16*N; cout << "(Warning): Pmax<Nmax => modify the value to "<< Pmax << endl;}
 	tstack_push("TestLaguerre2Bessel");
 	TestLaguerre2Bessel();
  	tstack_pop("TestLaguerre2Bessel");

src/laguerre2bessel.cc

@@ -7,6 +7,8 @@
 
 #include "walltimer.h"    //timing
 
+#include <omp.h> //TEST 
+
 namespace LagSHT {
 
   /* Version 0: test de calcul des J_{lnp} via une quadrature de Laguerre 
@@ -60,50 +62,40 @@ namespace LagSHT {
       facts[n] = facts[n-1]*sqrt(((r_8)n)/((r_8)(n+alpha)) ); 
     }
 
-//     cout << "Rmax, Lmax, Nmax, Pmax: " 
-// 	 << Rmax_ << ", " 
-// 	 << Lmax_ << ", " 
-// 	 << Nmax_ << ", " 
-// 	 << Pmax_ << endl;
 
     jnu_ = new Bessel(Lmax_,Pmax_);
     int Jstride = Lmax_*Pmax_; // nbre of (l,p) couples
     Jlnp_.resize(Nmax_*Jstride);
 
     
-    vector<r_8> vTest = lagTrans_->GetNodes();
-    
     r_8 integLowerBound = 0.;
-    r_8 integUpperBound = 2*Rmax_/tau; //for the integral 
+    r_8 integUpperBound = 10*Rmax_/tau; //for the integral 
     r_8 tol = 1e-6;
 
     tstack_push("J Integ compute");
 
+    Quadrature<r_8,Integrator_t>::values_t loc_val;
     for(int p=0; p<Pmax_; p++){
       for(int l=0; l<Lmax_; l++){
+
 	int ploff7 = p + l*Pmax_;
  	r_8 klp = (*jnu_)(l,p)/Rmax_*tau; //tau-rescaling of the original klp definition
+
 	for (int n=0; n<Nmax_; n++){
+
 	  LaguerreFuncQuad* Ln = new LaguerreFuncQuad(n);
+
 	  IntFunc1D f(l,klp,Ln);
 	  theQuad.SetFuncBounded(f,integLowerBound,integUpperBound);
-
-	  Quadrature<r_8,Integrator_t>::values_t loc_val = theQuad.LocalAdapStrat(tol);
+	  loc_val = theQuad.LocalAdapStrat(tol);
 	  Jlnp_[ploff7 + n*Jstride] = loc_val.first*facts[n];
 
 	  delete Ln;
 	}//n
       }//l
     }//p
-    tstack_pop("J Integ compute");
-
-    /*
-      std::ofstream ofs;
-      ofs.open ("jlnp.txt", std::ofstream::out);
-      Print(ofs,1);
-      ofs.close();
-    */
 
+    tstack_pop("J Integ compute");
 
   }//ComputeJlpnInteg
 

src/laguerreBuilder.h

@@ -33,7 +33,7 @@
 using namespace std;
 
 
-#define DEBUG 1
+#define DEBUG 0
 
 namespace LagSHT {
 
@@ -219,13 +219,13 @@ protected:
     int alphaFact = 1;
     for(int i=1;i<=alpha_; i++) alphaFact *= i;
     
+
+#if DEBUG >= 1  
     r_16 sumRTh = (r_16)(N*(N+alpha_));
     r_16 sumWTh = (r_16)(alphaFact);
 
-#if DEBUG >= 1  
     cout << "Sum roots = " << sumR << " -> diff with theory: " << sumR - sumRTh << endl;
     cout << "Sum weights = " << sumW << " -> diff with theory: " << sumW - sumWTh << endl;
-#endif
     
     if (fabs(sumR-sumRTh)>(0.000001*sumRTh)) {
       cout << " !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n" 
@@ -242,6 +242,7 @@ protected:
 	   << " !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n" 
 	   << endl;
     }
+#endif
     
   }//nodes_and_weights
 

src/quadinteg.h

@@ -24,8 +24,6 @@
  *  main author: J.E Campagne
  */
 
-
-
 // include standard c/c++
 #include <cmath>
 #include <cfloat>
@@ -40,20 +38,22 @@
 #include <string>
 
 #include <limits> //for precision features
+#include <iterator> //ostream_iterator
 
 #include "walltimer.h"    //timing
 
-#define NAMESPACETAG LagSHT
 
-using namespace std;
+#define DUMPLEV 0
 
-namespace NAMESPACETAG {
+using namespace std; //See if better on each routine
 
+namespace LagSHT {
 
-/* 1D function Class (From Sophya)
+
+/* 1D function Class
 */
 
-/*
+/*!
   \class GenericFunction 
   \brief Base class definition for objects representing functions 
 */
@@ -97,7 +97,7 @@ class Function1D : public ClassFunc1D {
 typedef Function1D Function; 
 
 
-/*
+/*!
   \class Quadrature
   \brief quadrature defined in [0,1]
 */
@@ -110,7 +110,7 @@ public:
   typedef pair< bounds_t , values_t > result_t; // { {a,b},{int,err} }
   typedef list< result_t > colresults_t; // collection of result_t
   
-  /*
+  /*!
     Utility  Sort funtion according to error (largest first)  
    */
   struct Comparator {
@@ -121,8 +121,39 @@ public:
     }
   };
 
+  /*! Class for intermediate recursive result
+   */
+  struct ResultClass {
+    T integral;
+    T error;
+    vector<T> scaledNodes;
+
+    ResultClass(T integ, T err, const vector<T>& vec): 
+      integral(integ), error(err), scaledNodes(vec) {}
+
+    ResultClass(const ResultClass& other):
+      integral(other.integral),error(other.error),scaledNodes(other.scaledNodes){}
+    
+    ResultClass& operator=(const ResultClass& other){
+      if(this == &other) return *this;
+      integral = other.integral;
+      error    = other.error;
+      scaledNodes = other.scaledNodes;
+      return *this;
+    }
+    
+
+    void Print() {
+      cout << "Integral = " << integral << " +/- " << error <<endl;
+      cout << "Nodes <" << setprecision(10); 
+      copy(scaledNodes.begin(), scaledNodes.end(), ostream_iterator<T>(cout, ", "));
+      cout << ">" << setprecision(6)<<endl;
+    }
+    
+  };
+
 public:  
-    /*
+    /*!
       Default Ctor
     */
   Quadrature(string name="",size_t npts=0,string fName="", bool init=false):
@@ -147,14 +178,14 @@ public:
 	  fout.close();
 	}//save
 
-      } catch (::exception& ex) {
+      } catch (exception& ex) {
 	string msg = "Init ";
 	msg += name + ": ERROR ";
 	msg += ex.what();
-	::cerr << msg << endl;
+	cerr << msg << endl;
 	exit(EXIT_FAILURE);
       } catch (...) {
-	::cerr << "Init " << name << " unexpected error " << endl;
+	cerr << "Init " << name << " unexpected error " << endl;
 	exit(EXIT_FAILURE);
       }
     } else {
@@ -170,21 +201,21 @@ public:
 	  inFile >> absc_[i]>> absw_[i]  >> errw_[i];
 	}
 	inFile.close();
-      } catch (::ifstream::failure e) {
-	::cerr << name <<" Exception opening/reading/closing file <"<<fName<<">"<<::endl;
+      } catch (ifstream::failure e) {
+	cerr << name <<" Exception opening/reading/closing file <"<<fName<<">"<<endl;
 	exit(EXIT_FAILURE);
       } //try/catch read
 
     }//initialisation
  
   }//Ctor
-  /*
+  /*!
     Dtor: care cannot be virtual because of template
    */
   ~Quadrature() {}
     
    
-  /*
+  /*!
     Rescale abscissa and weights from [xInmin, xInmax] to [xOutmin, xOutmax]
     Typical use-case the abs & weights are computed for (-1,1) range
     and we need to translate them to (0,1) or vice-versa
@@ -202,15 +233,15 @@ public:
     }
   }//RescaleAbsWeights
   
-   /*
+   /*!
      Set function
    */
   void SetFunction(ClassFunc1D& aFunc) { func_ = &aFunc;}
-  /*
+  /*!
     Set the integration bounds
   */
   void SetBounds(T xmin, T xmax) { initBounds_ = make_pair(xmin,xmax); }
-  /*
+  /*!
     Set both Function & bounds
   */
   void SetFuncBounded(ClassFunc1D& aFunc, T xmin, T xmax) {
@@ -218,31 +249,31 @@ public:
     SetBounds(xmin,xmax);
   }
   
-  /*
+  /*!
     Get abcissa, abscissa weights, error weights 
   */
   vector<T> GetAbscissa()  const {return absc_;}
   vector<T> GetAbscissaW() const {return absw_;}
   vector<T> GetErrorW()    const {return errw_;}
 
-  /* 
+  /*!
      Get number of points used
    */
   size_t GetNPoints() const {return npts_;}
   size_t GetOrder()   const {return 0;} //should be reimplemented
 
   
-  /*
+  /*!
     Debug
   */
-  void ToOut(::ostream& out, const result_t& aResult) {
+   void ToOut(ostream& out, const result_t& aResult) {
      out << "{" << aResult.first.first << ",," 
 	  <<  aResult.first.second << "}, "
 	  <<  aResult.second.first << ", "
 	  << aResult.second.second << "}";
    }//ToOut
   
-  void ToOut(::ostream& out, const colresults_t& aColRes) {
+  void ToOut(ostream& out, const colresults_t& aColRes) {
     typename colresults_t::const_iterator iRes,iEnd;
     iEnd = aColRes.end();
     out << "{";
@@ -252,80 +283,89 @@ public:
     out << "}" << endl;
   }//ToOut
   
-   /*
+   /*!
      Compute func integral between the current bounds
      The abscissa/weights of the Quadratures are defined for (0,1) interval.
    */
-  values_t IntEval(const bounds_t& curBounds) {
-
+  ResultClass IntEval(const bounds_t& curBounds) {
     T integral=0.;
     T error   = 0.;
     T boundsDist = curBounds.second - curBounds.first;
+    vector<T>xi(npts_);
+    
     for (size_t i=0;i<npts_;i++){
-      T xi = curBounds.first + absc_[i]*boundsDist;
-      T fi = (*func_)(xi);
+      xi[i] = curBounds.first + absc_[i]*boundsDist;
+      T fi = (*func_)(xi[i]);
       integral += absw_[i] * fi;
       error    += errw_[i] * fi;
     }
-
-    return make_pair(integral*boundsDist,fabs((double)error*boundsDist));
+    
+    return ResultClass(integral*boundsDist,fabs((double)error*boundsDist),xi);
   }//IntEval
   
 
-  /* Recursive computation of integral on a finite range
+  /*! Recursive computation of integral on a finite range
     \input curBounds = {x0,x1}
     \input cut: threshold to be compared to the integral error to stop or not
     \output : {integral, error}
+  */
+  ResultClass LocalAdapStratRecur(const bounds_t& curBounds, T cut, int maxdepth, int depth=0) {
 
-    Todo: 1) improve the partitionning (eg. Mathematica uses the quadrature abscissa rescaled) 
-          2) revisit the condition
-   */
-  values_t LocalAdapStratRecur(const bounds_t& curBounds, T cut){
+#if DUMPLEV > 2
+    cout << "depth ["<<depth<<"]";
+    for(int i=0;i<depth;i++) cout << "...";
+#endif
 
-    
-    values_t curValues = IntEval(curBounds); //{int,err}
-    //Not used    T integral = curValues.first;
-    T error    = curValues.second;
-    if (error <= cut) return curValues;
-    
-    //split the initial range in two partition and decrease by 2 the cut on each part 
-    T xmiddle = (curBounds.first+curBounds.second)*0.5; //xmid
-    bounds_t lowerB = make_pair(curBounds.first,xmiddle); //{x0,xmid}
-    values_t valLow = LocalAdapStratRecur(lowerB,cut/2.0); 
-    
-    bounds_t upperB = make_pair(xmiddle,curBounds.second); //{xmid,x1}
-    values_t valUp = LocalAdapStratRecur(upperB,cut/2.0);
+    ResultClass curResult = IntEval(curBounds); //{int, err, nodes}
 
-    //total integral= sum of sub integral, idem for the error
-    return make_pair(valLow.first+valUp.first,valLow.second+valUp.second);
-    
+    T error    = curResult.error;
+  
+    if (error <= cut || depth >= maxdepth){
+#if DUMPLEV > 2
+      cout << "Int[" <<curBounds.first << ", " <<curBounds.second << "]: "
+	   << curResult.integral << " +/- " << curResult.error << ".... Ok" << endl;
+#endif
+      if(depth == maxdepth) cout << "Warning, max recursive reached" <<endl;
+      return curResult;
+    }  
+    int nnodes = curResult.scaledNodes.size();
+    int lastNodes=nnodes-1;
+    T integral = (T)0.;
+    error      = (T)0.;
+    for(int i=0;i<lastNodes;i++){
+      ResultClass locres = 
+	LocalAdapStratRecur(make_pair(curResult.scaledNodes[i],curResult.scaledNodes[i+1]),cut,maxdepth,depth+1);
+      integral += locres.integral;
+      error    += locres.error;
+    }
 
+    return ResultClass(integral, error, curResult.scaledNodes);
+        
   }//LocalAdapStratRecur
 
-
-  /*
+  /*!
     Local Adaptative Strategy
    */
   values_t LocalAdapStrat(T tol = (T)1.0e-6, 
-			  size_t MaxNumRecurssion = 1000, // not used 
-			  size_t MinNumRecurssion = 5    // not used
+			  size_t MaxNumRecurssion = 100
 			  ) {
-
+    
     bounds_t curBounds = initBounds_; // {a,b}
-    values_t curValues = IntEval(curBounds); //{int,err}