Changeset 3510 in Sophya for trunk

Timestamp:

Aug 8, 2008, 3:11:45 PM (17 years ago)

Author:

ansari

Message:

Nettoyage alm.h .cc , utilisation RandomGenerator par Alm<T> et SphericalTransformServer , Reza 08/08/2008

Location:

trunk/SophyaLib/Samba

Files:

• alm.cc (modified) (6 diffs)
• alm.h (modified) (2 diffs)
• sphericaltransformserver.cc (modified) (4 diffs)
• sphericaltransformserver.h (modified) (4 diffs)

trunk/SophyaLib/Samba/alm.cc

@@ -1 +1 @@
 #include "sopnamsp.h"
 #include "alm.h"
+
+
+/*!
+  \class Alm 
+  \ingroup Samba 
+  Class for manipulating the coefficients \f$a_{lm}\f$ of the development
+  in spherical harmonics of  a function efined on a sphere.
+*/
+ 
+/*!
+  fwhm specifies the gaussian beam half witdh in arc.minutes
+*/
 template <class T>
 Alm<T>::Alm(const TVector<T>& clin, const r_8 fwhm) 
 
 {
+  int_4 nlmax= clin.NElts()-1;
+
+  //alm.ReSize(nlmax);
+  this->ReSizeRow(nlmax+1);
+  RandomGenerator rg(1, false);
+  GenFromCl(clin, fwhm, rg);
+}
+
+/*!
+  fwhm specifies the gaussian beam half witdh in arc.minutes
+*/
+template <class T>
+Alm<T>::Alm(const TVector<T>& clin, const r_8 fwhm, RandomGenerator & rg)
+{
+  int_4 nlmax= clin.NElts()-1;
+
+  //alm.ReSize(nlmax);
+  this->ReSizeRow(nlmax+1);
+  GenFromCl(clin, fwhm, rg);  
+}
 
 
+
+template <class T>
+void Alm<T>::GenFromCl(const TVector<T> & clin, const r_8 fwhm, RandomGenerator & rg)
+{
   /*=======================================================================
      creates the a_lm from the power spectrum, 
@@ -21 +57 @@
   int_4 nlmax= clin.NElts()-1;
 
-  //alm.ReSize(nlmax);
-  this->ReSizeRow(nlmax+1);
-
   r_8 sig_smooth = fwhm/sqrt(8.*log(2.))/(60.*180.)* M_PI;
   int_4 n_l = nlmax+1;
@@ -29 +62 @@
 
   //    --- smoothes the initial power spectrum ---
-  TVector<T> cl=clin;
+  TVector<T> cl(clin, false);
   int l;
   for (l=0;l<n_l;l++)
@@ -45 +78 @@
       T rms=sqrt(cl(l));
       //        ------ m = 0 ------
-      complex<T> zeta1(NorRand());
+      complex<T> zeta1((T)rg.Gaussian() );
       (*this)(l,0)   = zeta1 * rms;
 
@@ -52 +85 @@
         {
           complex<T> aux1(hsqrt2);
-          complex<T> aux2(NorRand(),NorRand());
+          complex<T> aux2((T)rg.Gaussian() , (T)rg.Gaussian() );
           zeta1=aux1*aux2;
           (*this)(l,m)=rms*zeta1;
@@ -79 +112 @@
   return powsp;
 }
+
+/*!
+  \class Bm
+  \ingroup Samba 
+  Class for a vector with an index running from \f$-m_{max}\f$ to \f$+m_{max}\f$ 
+  (then the size of the vector will be actually \f$2m_{max}+1)\f$.
+  This class is used by the spherical harmonics transform server 
+*/
+
  
 #ifdef __CXX_PRAGMA_TEMPLATES__

trunk/SophyaLib/Samba/alm.h

@@ -6 +6 @@
 #define ALM_SEEN
 
-#include "srandgen.h"
+#include "stsrand.h"
 #include "nbmath.h"
 #include "triangmtx.h"
@@ -13 +13 @@
 namespace SOPHYA {  
 
-/*! class for the coefficients \f$a_{lm}\f$ of the development of a function defined on a sphere, in spherical harmonics */
+// ----- Classe Alm 
 template <class T>
-class Alm : public TriangularMatrix<complex<T> >
-  {
-    public :
+class Alm : public TriangularMatrix< complex<T> > {
+public :
  
-Alm() : TriangularMatrix<complex<T> >()  {;};
-    /* instanciate the class from the maximum value of the index \f$l\f$ in the sequence of the \f$a_{lm}\f$'s */
-Alm(sa_size_t nlmax) : TriangularMatrix<complex<T> >(nlmax+1) {;}
-Alm(const Alm<T>& a,  bool share=false)  : TriangularMatrix<complex<T> >(a, share)  {;}
+  Alm() : TriangularMatrix<complex<T> >()  {}
+//! instanciate the class from the maximum value of the index \f$l\f$ in the sequence of the \f$a_{lm}\f$'s 
+  Alm(sa_size_t nlmax) : TriangularMatrix<complex<T> >(nlmax+1) {}
+//! copy constructor 
+  Alm(const Alm<T>& a,  bool share=false)  : TriangularMatrix<complex<T> >(a, share)  {}
+//! Constructor, generates alm from a cl vector
+  Alm(const TVector<T> & clin, const r_8 fwhm) ;
+//! Constructor, generates alm from a cl vector, using the specified random generator object
+  Alm(const TVector<T> & clin, const r_8 fwhm, RandomGenerator & rg);
 
-Alm(const TVector<T>& clin, const r_8 fwhm) ;
-/*! resize with a new lmax */
-inline void ReSizeToLmax(sa_size_t nlmax) {this->ReSizeRow(nlmax+1);}
-inline sa_size_t Lmax() const {return this->rowNumber()-1;}
-TVector<T> powerSpectrum() const;
+//! Resize with a new lmax 
+  inline void ReSizeToLmax(sa_size_t nlmax) {this->ReSizeRow(nlmax+1);}
+  inline sa_size_t Lmax() const {return this->rowNumber()-1;}
 
-inline  Alm<T>& SetComplexT(complex<T> a)
-   {
-     this->SetT(a);
-     return *this;
-   }
- inline Alm<T>& operator = (complex<T> a)
- {
-   return SetComplexT(a);
- } 
+//! Computes the corresponding pwer spectrum (Cl vector)
+  TVector<T> powerSpectrum() const;
 
+  inline  Alm<T>& SetComplexT(complex<T> a)
+    {  this->SetT(a);  return *this;  }
+  inline Alm<T>& operator = (complex<T> a)
+    {  return SetComplexT(a);  } 
+
+private:
+  void GenFromCl(const TVector<T> & clin, const r_8 fwhm, RandomGenerator& rg);
 };
-/*! class for a vector with an index running from \f$-m_{max}\f$ to \f$+m_{max}\f$ (then the size of the vector will be actually \f$2m_{max}+1)\f$ */
+
+
+//------------ Classe Bm , vecteur avec index  -mMax <= index <= +mMax 
 template <class T>
-class Bm
-  {
-    public :
-Bm(): nmmax_(0) {;};
-    /* instanciate from the maximum absolute value of the index */
-Bm(sa_size_t mmax) : nmmax_(mmax)   { bm_.ReSize( (2*mmax+1) );}
- Bm(const Bm<T>& b,  bool share=false) : nmmax_(b.nmmax_), bm_(b.bm_, share) {;}
- /*! resize with a new value of mmax */
-inline void ReSizeToMmax(sa_size_t mmax) 
-  {
-    nmmax_= mmax;
-    bm_.ReSize(2* nmmax_+1);
-  }
-//inline sa_size_t Size() const {return 2*nmmax_+1;};
-inline T& operator()(sa_size_t m) {return bm_( adr_i(m));};
-inline T const& operator()(sa_size_t m) const {return *(bm_.Begin()+ adr_i(m));};
-/*! return the current value of the maximum absolute value of the index */
-inline sa_size_t Mmax() const 
-   {
-     return nmmax_;
-   }
-inline Bm<T>& operator = (const Bm<T>& b)
- {
-   if (this != &b)
-     {
+class Bm  {
+public :
+  Bm(): nmmax_(0) {;};
+//! Constructor with the specification of the maximum absolute value of the index : |mMax|
+//! Copy contrsuctor 
+  Bm(sa_size_t mmax) : nmmax_(mmax)   { bm_.ReSize( (2*mmax+1) );}
+  Bm(const Bm<T>& b,  bool share=false) : nmmax_(b.nmmax_), bm_(b.bm_, share) { }
+ //! resize with a new value of mmax 
+  inline void ReSizeToMmax(sa_size_t mmax) 
+    { nmmax_= mmax; bm_.ReSize(2* nmmax_+1); }
+
+//! Acces to element with index m  ( -mMax <= m <= +mMax  ) 
+  inline T operator()(sa_size_t m) const { return bm_( m+nmmax_); }
+//! Acces to element with index m  ( -mMax <= m <= +mMax  ) 
+  inline T& operator()(sa_size_t m) { return bm_( m+nmmax_); }
+
+//! Return the current value of the maximum absolute value of the index 
+  inline sa_size_t Mmax() const   { return nmmax_; }
+
+  inline Bm<T>& operator = (const Bm<T>& b)
+    {
+    if (this != &b) {
        nmmax_= b.nmmax_;
        bm_=  b.bm_;
-     }
-   return *this;
- }
-
-  private:
-/*! return the address in the array representing the vector */
-inline sa_size_t adr_i(sa_size_t i) const 
-{
-  return (i+nmmax_);
-}
+       }
+    return *this;
+    }
 
 
-
- sa_size_t nmmax_;
- NDataBlock<T> bm_;
-
+ private:
+  sa_size_t nmmax_;
+  NDataBlock<T> bm_;
   };
 

trunk/SophyaLib/Samba/sphericaltransformserver.cc

@@ -6 +6 @@
 #include "sphericaltransformserver.h"
 #include "tvector.h"
-#include "srandgen.h"
+#include "stsrand.h"
 #include "nbmath.h"
 #include "timing.h"
@@ -152 +152 @@
 
  */
+
+//!  Default constructor - Creates a non thread-safe RandomGenerator to be used by GenerateFromCl 
+template<class T>
+SphericalTransformServer<T>::SphericalTransformServer()
+: rg_(1, false) 
+{
+  fftIntfPtr_=new FFTPackServer(true); // preserveinput = true
+  fftIntfPtr_->setNormalize(false);
+}
+
+//!  Constructor with the specification of a RandomGenerator object to be used by GenerateFromCl 
+template<class T>
+SphericalTransformServer<T>::SphericalTransformServer(RandomGenerator const & rg)
+: rg_(rg) 
+{
+  fftIntfPtr_=new FFTPackServer(true); // preserveinput = true
+  fftIntfPtr_->setNormalize(false);
+}
+
+template<class T>
+SphericalTransformServer<T>::~SphericalTransformServer()
+{
+  if (fftIntfPtr_!=NULL) delete fftIntfPtr_; 
+}
+
+/*!
+ Set a fft server. The constructor sets a default fft server (fft-pack). 
+ So it is not necessary to call this method for a standard use.
+ \warning The FFTServerInterface object should NOT overwrite the input arrays
+*/
+template<class T>
+void SphericalTransformServer<T>::SetFFTServer(FFTServerInterface* srv) 
+{
+  if (fftIntfPtr_!=NULL) delete fftIntfPtr_;
+  fftIntfPtr_=srv;
+  fftIntfPtr_->setNormalize(false);
+}
+
 
  /*! \fn void SOPHYA::SphericalTransformServer::GenerateFromAlm( SphericalMap<T>& map, int_4 pixelSizeIndex, const Alm<T>& alm) const
@@ -1275 +1313 @@
   //  Alm<T> a2lme = almFromCl(Cle, fwhm);
   // Alm<T> a2lmb = almFromCl(Clb, fwhm);
-  Alm<T> a2lme(Cle, fwhm);
-  Alm<T> a2lmb(Clb, fwhm);
+  Alm<T> a2lme(Cle, fwhm, rg_);
+  Alm<T> a2lmb(Clb, fwhm, rg_);
 
   GenerateFromAlm(sphq,sphu,pixelSizeIndex,a2lme,a2lmb); 
@@ -1293 +1331 @@
 {
 
-  Alm<T> alm(Cl, fwhm);
+  Alm<T> alm(Cl, fwhm, rg_);
   GenerateFromAlm(sph,pixelSizeIndex, alm ); 
 }

trunk/SophyaLib/Samba/sphericaltransformserver.h

@@ -5 +5 @@
 #include "fftservintf.h"
 #include "fftpserver.h"
+#include "stsrand.h"
 #include "alm.h"
 #include "lambdaBuilder.h"
@@ -18 +19 @@
  public:
 
-  SphericalTransformServer()  
-    {
-      fftIntfPtr_=new FFTPackServer(true); // preserveinput = true
-      fftIntfPtr_->setNormalize(false);
-    };
-  ~SphericalTransformServer(){ if (fftIntfPtr_!=NULL) delete fftIntfPtr_;};
+  SphericalTransformServer();  
+  SphericalTransformServer(RandomGenerator const & rg);  
+  virtual ~SphericalTransformServer(); 
+  void SetFFTServer(FFTServerInterface* srv=NULL);
   
- /*!
- Set a fft server. The constructor sets a default fft server (fft-pack). So it is not necessary to call this method for a standard use.
- \warning The FFTServerInterface object should NOT overwrite the input arrays
-*/
-  void SetFFTServer(FFTServerInterface* srv=NULL) 
-    {
-      if (fftIntfPtr_!=NULL) delete fftIntfPtr_;
-      fftIntfPtr_=srv;
-      fftIntfPtr_->setNormalize(false);
-    }
 
   void GenerateFromAlm( SphericalMap<T>& map, int_4 pixelSizeIndex, const Alm<T>& alm) const; 
@@ -44 +33 @@
                      int_4 pixelSizeIndex,
                      const TVector<T>& Cle, const TVector<T>& Clb, 
-                    const r_8 fwhm) const; 
+                     const r_8 fwhm) const; 
 
 
@@ -118 +107 @@
 
  FFTServerInterface* fftIntfPtr_;
+ mutable RandomGenerator rg_;
 
 };