(JEC) 17/4/15 use ANN library to find theta,phi to spherical index...

(JEC) 17/4/15 use ANN library to find theta,phi to spherical index (lagSphericalTransform) , and a rescalng for radial to index (laguerreTransform)

Makefile

@@ -4,9 +4,9 @@ ifeq ($(MACH), Linux)
 # ===== OpenBLAS ========
 		BLASYES = -DCBLAS
 		BLASDIR = /opt/OpenBLAS
-		BLASINC = $(BLASDIR)/include
+		BLASINC = -I$(BLASDIR)/include
 		BLASLIB = -L$(BLASDIR)/lib
-		BLASLIBN = -lopenblas
+		BLASLIBN = $(BLASLIB) -lopenblas
 	endif
 else 
 	include Darwin_g++_omp_make.inc
@@ -15,14 +15,20 @@ else
 	BLASLIBN = -framework Accelerate
 endif
 
+# ===== ANN  ==========
+ANNDIR = /Users/campagne/Travail/kits/ann_1.1.2/
+ANNLIB = $(ANNDIR)/lib
+ANNINC = -I$(ANNDIR)/include
+ANNLIBN = -L$(ANNLIB) -lANN
+
 
 # ===== Libsharp =====
 #SHARPDIR = <yourdir>/libsharp-code/auto/
 
 SHARPDIR = /Users/campagne/Travail/kits/libsharp-code/auto
 SHARPLIB = $(SHARPDIR)/lib
-SHARPINC = $(SHARPDIR)/include
-SHARPLIBN = -lsharp -lc_utils -lfftpack 
+SHARPINC = -I$(SHARPDIR)/include
+SHARPLIBN = -L$(SHARPLIB) -lsharp -lc_utils -lfftpack 
 
 
 OBJ = ./Objs/
@@ -77,8 +83,8 @@ CXXHDR = lagsht_exceptions.h \
 
 
 
-CPPFLAGS  += $(BLASYES) -I$(SHARPINC) -I$(BLASINC)
-LDFLAGS += -L$(SHARPLIB) $(SHARPLIBN) $(BLASLIB) $(BLASLIBN)  -lm
+CPPFLAGS  += $(BLASYES) $(SHARPINC) $(BLASINC) $(ANNINC)
+LDFLAGS   += $(SHARPLIBN) $(BLASLIBN) $(ANNLIBN)  -lm
 
 
 

lagSphericTransform.cc

@@ -66,6 +66,10 @@ void BaseLagSphTransform::SetThetaPhiMap(string choice, int nrings, int nphi){
   Npix_ = get_npix(sharpjob_.get_geom_info());    
   Nalm_ = get_nalms(sharpjob_.get_alm_info());
 
+  //fill the ANN-kdTree to prepare the filling of the 2D spheres
+  geom_t vecPix;
+  vecPix = GetPixelPositions();
+  ann.SetkdTree(vecPix);
 
   geomDone_ = true;
     
@@ -76,12 +80,9 @@ void BaseLagSphTransform::SetThetaPhiMap(string choice, int nrings, int nphi){
        << Nalm_ 
        << endl;
 
-
-
-  
 }//SetThetaPhiMap
 
-BaseLagSphTransform::geom_t BaseLagSphTransform::GetPixelPositions() {
+geom_t BaseLagSphTransform::GetPixelPositions() {
   geom_t pixels;
   sharp_geom_info *ginfo; 
   ginfo = sharpjob_.get_geom_info();

lagSphericTransform.h

@@ -10,6 +10,8 @@ using namespace std;
 #include "sharp_almhelpers.h"
 #include "sharp_cxx.h"
 
+//ANN
+#include <ANN/ANN.h>
 
 #include "laguerreTransform.h"
 #include "lagsht_exceptions.h"
@@ -22,17 +24,95 @@ namespace LagSHT {
 //// --------- Class Spherical Harmonic Laguerre Transform -- //
 ////////////////////////////////////////////////////////////////
 
+  typedef pair<r_8,r_8>    pixpos_t;
+  typedef vector<pixpos_t> geom_t;
+
+
+  class ANNServer {
+    
+  public:
+    //! Constructor
+    ANNServer(): nNear_(1), dim_(2), eps_(0), 
+		 npts_(0), dataPts_(NULL), kdTree_(NULL), nnIdx_(NULL), dists_(NULL) {
+      
+    }
+
+    //! Constructor
+    ANNServer(geom_t vec):
+      nNear_(1), dim_(2), eps_(0),npts_(0), dataPts_(NULL), kdTree_(NULL) {
+      nnIdx_ = new ANNidx[nNear_];  // allocate near neigh indices
+      dists_ = new ANNdist[nNear_]; // allocate near neighbor dists      
+    }//Ctor
+
+    //! Set kdtree
+    void SetkdTree(geom_t vec) {
+
+      npts_ = vec.size();
+      dataPts_ = annAllocPts(npts_, dim_);
+      if(dataPts_ == NULL) throw LagSHTError("ANN Failure dataPts alloc");
+
+      for(int i=0;i<npts_;i++){
+	ANNpoint p = dataPts_[i];
+	p[0] = vec[i].first;
+	p[1] = vec[i].second;
+      }
+      kdTree_ = new ANNkd_tree(dataPts_,npts_,dim_);
+      if(kdTree_ == NULL) throw LagSHTError("ANN Failure kdTree alloc");
+      
+    }//SetkdTree
+
+    //! Destructor
+    virtual ~ANNServer() {
+      if (dataPts_) annDeallocPts(dataPts_);
+      if (nnIdx_) delete [] nnIdx_;	
+      if (dists_) delete [] dists_;
+      if (kdTree_) delete kdTree_;      
+      annClose();
+    }//Dtor
+
+    int GetClosestId(r_8 theta, r_8 phi) const {
+
+      ANNpoint queryPt;
+      queryPt = annAllocPt(dim_);
+      if(queryPt == NULL) throw LagSHTError("Get Failure queryPt alloc");
+      queryPt[0] = theta;
+      queryPt[1] = phi;
+      
+      if(kdTree_ == NULL) throw LagSHTError("Get Failure kdTree alloc");
+      if(nnIdx_ == NULL) throw LagSHTError("Get Failure nnIdx alloc");
+      if(dists_ == NULL) throw LagSHTError("Get Failure dists alloc");
+      
+      kdTree_->annkSearch(queryPt,nNear_,nnIdx_,dists_,eps_);
+      annDeallocPt(queryPt);
+      
+      return nnIdx_[0]; //first closest
+    }//GetClosestId
+
+
+  protected:
+    int nNear_; //nomber of solution to the closest point
+    int dim_;   //dimsion of the point space
+    r_8 eps_;   // error bound
+    int npts_;  //number of points to fill the kdTree
+
+    ANNkd_tree* kdTree_;  //the tree
+    ANNidxArray  nnIdx_; // near neighbor indices
+    ANNdistArray dists_; // near neighbor distances
+    ANNpointArray dataPts_; //the initial grid points
+
+  };//ANNServer
+
 class BaseLagSphTransform {
  public:
 
-  typedef pair<r_8,r_8>    pixpos_t;
-  typedef vector<pixpos_t> geom_t;
 
   //! Creator
   BaseLagSphTransform(int Lmax, int Mmax, int Nmax, int alpha=2):
     L_(Lmax), M_(Mmax), N_(Nmax), alpha_(alpha),
     Nrings_(0), Nphi_(0), Npix_(0), Nalm_(0),geomDone_(false) {
+
      if (Mmax<0) M_=Lmax;
+              
   }
   //! Destructor
   virtual ~BaseLagSphTransform() {}
@@ -48,7 +128,7 @@ class BaseLagSphTransform {
   //! Total number of Alm coefficients
   int Nalm()   { return Nalm_; }
   
-  //! List of pixel positions
+  //! Sorted List of pixel center (theta, phi) positions
   geom_t GetPixelPositions();
     
   /*! Choose the 2D map representation & initialization of libsharp
@@ -58,6 +138,12 @@ class BaseLagSphTransform {
    */
   void SetThetaPhiMap(string choice="Fejer1", int nrings=-1, int nphi=-1);
 
+  /*! Get Pixel index in the 2D map
+    \input theta in [0, pi]
+    \input phi   in [0, 2pi]
+  */
+  int PixIndexSph(r_8 theta, r_8 phi){ return ann.GetClosestId(theta, phi); }
+
 
 protected:
 
@@ -93,7 +179,9 @@ protected:
 
   sharp_cxxjob<r_8> sharpjob_; //!< libsharp job // voir si r_8 -> r_16 extend the validity of the computation
 
-};
+  ANNServer ann; //!< utility to fill the spheric shells
+
+}; //BaseLagSphTransform
 
 class LaguerreSphericalTransform : public BaseLagSphTransform {
  public:

laguerreTransform.cc

@@ -73,6 +73,15 @@ LaguerreTransform::LaguerreTransform(int N, r_8 R, int alpha) :
 
 }//Ctor
 
+int LaguerreTransform::RadIndex(r_8 r) const {
+  using namespace std;
+  r_8 rscaled = r*nodes_[N_-1]/R_;
+  
+  vector<r_8> dist(N_);
+  transform(nodes_.begin(),nodes_.end(),dist.begin(), distance(rscaled));
+  return min_element(dist.begin(),dist.end()) - dist.begin();
+}//RadIndex
+
 
 void LaguerreTransform::MultiAnalysis(const vector< complex<r_8> >& fi, 
 				      vector< complex<r_8> >& fn, int stride) {

laguerreTransform.h

@@ -33,15 +33,17 @@ protected:
 class LaguerreTransform : public BaseLaguerreTransform {
  public:
   //! Creator
-  LaguerreTransform(int n, r_8 R, int alpha=2);
+  LaguerreTransform(int N, r_8 R, int alpha=2);
   //! Destructor
   virtual ~LaguerreTransform() {}
 
   //! Getters
-/*   vector<r_16> GetNodes()   const { return nodes_; } //const ? */
-/*   vector<r_16> GetWeights() const { return weights_; } //const ? */
-  vector<r_8> GetNodes()   const { return nodes_; } //const ?
-  vector<r_8> GetWeights() const { return weights_; } //const ?
+  vector<r_8> GetNodes()   const { return nodes_; }
+  vector<r_8> GetWeights() const { return weights_; }
+
+
+  //! Find the closest nodes
+  int RadIndex(r_8 r) const;
 
 
   //! Multi Analysis at once (Values -> Coefficients)
@@ -61,8 +63,6 @@ class LaguerreTransform : public BaseLaguerreTransform {
   */
   void MultiSynthesis(const vector< complex<r_8> >& fn, vector< complex<r_8> >& fi, int stride);
 
-
-
   //! Single Analysis (Values -> Coefs)
   /*!
     \param fi: values f(r_i) with r_i zero L^(alpha)_N
@@ -80,13 +80,20 @@ class LaguerreTransform : public BaseLaguerreTransform {
 protected:
 
   r_8 alphaFact_; //!< sqrt(alpha!)
-  r_8 R_; //!< maximal radius (not used yet)
+  r_8 R_; //!< maximal radius
 
-/*   vector<r_16> nodes_;    //!< nodes   of the Gauss-Laguerre quadrature used */
-/*   vector<r_16> weights_;  //!< weights of the Gauss-Laguerre quadrature used */
   vector<r_8> nodes_;    //!< nodes   of the Gauss-Laguerre quadrature used
   vector<r_8> weights_;  //!< weights of the Gauss-Laguerre quadrature used
 
+ public:
+  class distance {
+  public: 
+    distance(r_8 r0): r0_(r0) {}
+    r_8 operator()(r_8 x){ r_8 xn=x-r0_; return xn*xn;}
+  private:
+    r_8 r0_;
+  };//distance
+
 };//LaguerreTransform