guilhem_lavaux/cosmotool
0
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Changed interface of newRandomFourier to avoid generating copies.

Browse Source
This commit is contained in:
Guilhem Lavaux 2012-11-11 15:46:03 -05:00
parent ed1fede33d
commit 7ad4911872
5 changed files with 22 additions and 32 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
sample/CMakeLists.txt
View File

@ -47,7 +47,7 @@ target_link_libraries(testBSP ${tolink})
if (FFTW3_FOUND AND FFTW3F_FOUND AND EIGEN3_FOUND) if (FFTW3_FOUND AND FFTW3F_FOUND AND EIGEN3_FOUND)
add_executable(test_fft_calls test_fft_calls.cpp) add_executable(test_fft_calls test_fft_calls.cpp)
target_link_libraries(test_fft_calls ${tolink} ${FFTW3_LIBRARIES} ${FFTW3F_LIBRARIES}) target_link_libraries(test_fft_calls ${tolink} ${FFTW3_LIBRARIES} ${FFTW3F_LIBRARIES})
endif (FFTW3_FOUND AND EIGEN3_FOUND) endif (FFTW3_FOUND AND FFTW3F_FOUND AND EIGEN3_FOUND)
if (SHARP_LIBRARY AND SHARP_INCLUDE_PATH AND EIGEN3_FOUND) if (SHARP_LIBRARY AND SHARP_INCLUDE_PATH AND EIGEN3_FOUND)
include_directories(${SHARP_INCLUDE_PATH}) include_directories(${SHARP_INCLUDE_PATH})

3
sample/test_fft_calls.cpp
View File

@ -37,8 +37,7 @@ void test_2d(int Nx, int Ny)
dft.realSpace().scale(2.0); dft.realSpace().scale(2.0);
dft.fourierSpace().scale(0.2); dft.fourierSpace().scale(0.2);
typename SpectrumFunction<T>::FourierMapPtr m = spectrum.newRandomFourier(rng, dft.fourierSpace()); spectrum.newRandomFourier(rng, dft.fourierSpace());
dft.fourierSpace() = *m.get();
dft.synthesis(); dft.synthesis();
uint32_t dims[2] = { Ny, Nx }; uint32_t dims[2] = { Ny, Nx };

4
src/fourier/base_types.hpp
View File

@ -26,8 +26,8 @@ namespace CosmoTool
typedef FourierMap<std::complex<T> > FourierMapType; typedef FourierMap<std::complex<T> > FourierMapType;
typedef boost::shared_ptr<FourierMapType> FourierMapPtr; typedef boost::shared_ptr<FourierMapType> FourierMapPtr;
virtual FourierMapPtr virtual void
newRandomFourier(gsl_rng *rng, const FourierMapType& like_map) const = 0; newRandomFourier(gsl_rng *rng, FourierMapType& out_map) const = 0;
virtual void mul(FourierMapType& m) const = 0; virtual void mul(FourierMapType& m) const = 0;
virtual void mul_sqrt(FourierMapType& m) const = 0; virtual void mul_sqrt(FourierMapType& m) const = 0;

33
src/fourier/euclidian.hpp
View File

@ -27,12 +27,12 @@ namespace CosmoTool
{ {
} }
ptr_map newRandomFourier(gsl_rng *rng, const FourierMapType& like_map) const; void newRandomFourier(gsl_rng *rng, FourierMapType& out_map) const;
void mul(FourierMap<std::complex<T> >& m) const; void mul(FourierMapType& m) const;
void mul_sqrt(FourierMap<std::complex<T> >& m) const; void mul_sqrt(FourierMapType& m) const;
void mul_inv(FourierMap<std::complex<T> >& m) const; void mul_inv(FourierMapType& m) const;
void mul_inv_sqrt(FourierMap<std::complex<T> >& m) const; void mul_inv_sqrt(FourierMapType& m) const;
}; };
template<typename T> template<typename T>
@ -388,16 +388,13 @@ namespace CosmoTool
template<typename T> template<typename T>
typename EuclidianSpectrum_1D<T>::ptr_map void EuclidianSpectrum_1D<T>::newRandomFourier(gsl_rng *rng, FourierMapType& out_map) const
EuclidianSpectrum_1D<T>::newRandomFourier(gsl_rng *rng, const FourierMapType& like_map) const
{ {
typedef EuclidianFourierMapComplex<T> MapT; typedef EuclidianFourierMapComplex<T> MapT;
typedef typename EuclidianSpectrum_1D<T>::ptr_map ptr_map; typedef typename EuclidianSpectrum_1D<T>::ptr_map ptr_map;
typedef typename MapT::DimArray DimArray; typedef typename MapT::DimArray DimArray;
const MapT& m_c = dynamic_cast<const MapT&>(like_map); MapT& rand_map = dynamic_cast<MapT&>(out_map);
ptr_map ret_map = ptr_map(m_c.mimick());
MapT& rand_map = dynamic_cast<MapT&>(*ret_map.get());
std::complex<T> *d = rand_map.data(); std::complex<T> *d = rand_map.data();
long idx; long idx;
@ -405,7 +402,7 @@ namespace CosmoTool
long plane_size; long plane_size;
bool alleven = rand_map.allDimensionsEven(); bool alleven = rand_map.allDimensionsEven();
for (long p = 1; p < m_c.size(); p++) for (long p = 1; p < rand_map.size(); p++)
{ {
double A_k = std::sqrt(0.5*f(rand_map.get_K(p))); double A_k = std::sqrt(0.5*f(rand_map.get_K(p)));
d[p] = std::complex<T>(gsl_ran_gaussian(rng, A_k), d[p] = std::complex<T>(gsl_ran_gaussian(rng, A_k),
@ -415,14 +412,14 @@ namespace CosmoTool
d[0] = std::complex<T>(gsl_ran_gaussian(rng, std::sqrt(f(0))), 0); d[0] = std::complex<T>(gsl_ran_gaussian(rng, std::sqrt(f(0))), 0);
if (!rand_map.firstDimensionEven()) if (!rand_map.firstDimensionEven())
return ret_map; return;
// Correct the Nyquist plane // Correct the Nyquist plane
idx = dims[0]-1; // Stick to the last element of the first dimension idx = dims[0]-1; // Stick to the last element of the first dimension
d[idx] = std::complex<T>(d[idx].real() + d[idx].imag(), 0); d[idx] = std::complex<T>(d[idx].real() + d[idx].imag(), 0);
// 1D is special case // 1D is special case
if (dims.size() == 1) if (dims.size() == 1)
return ret_map; return;
plane_size = 1; plane_size = 1;
for (int q = 1; q < dims.size(); q++) for (int q = 1; q < dims.size(); q++)
@ -447,12 +444,10 @@ namespace CosmoTool
q += dims[0]-1; q += dims[0]-1;
d[q] = std::complex<T>(d[q].real()+d[q].imag(),0); d[q] = std::complex<T>(d[q].real()+d[q].imag(),0);
} }
return ret_map;
} }
template<typename T> template<typename T>
void EuclidianSpectrum_1D<T>::mul(FourierMap<std::complex<T> >& m) const void EuclidianSpectrum_1D<T>::mul(FourierMapType& m) const
{ {
EuclidianFourierMapComplex<T>& m_c = dynamic_cast<EuclidianFourierMapComplex<T>&>(m); EuclidianFourierMapComplex<T>& m_c = dynamic_cast<EuclidianFourierMapComplex<T>&>(m);
std::complex<T> *d = m.data(); std::complex<T> *d = m.data();
@ -462,7 +457,7 @@ namespace CosmoTool
} }
template<typename T> template<typename T>
void EuclidianSpectrum_1D<T>::mul_sqrt(FourierMap<std::complex<T> >& m) const void EuclidianSpectrum_1D<T>::mul_sqrt(FourierMapType& m) const
{ {
EuclidianFourierMapComplex<T>& m_c = dynamic_cast<EuclidianFourierMapComplex<T>&>(m); EuclidianFourierMapComplex<T>& m_c = dynamic_cast<EuclidianFourierMapComplex<T>&>(m);
std::complex<T> *d = m.data(); std::complex<T> *d = m.data();
@ -472,7 +467,7 @@ namespace CosmoTool
} }
template<typename T> template<typename T>
void EuclidianSpectrum_1D<T>::mul_inv(FourierMap<std::complex<T> >& m) const void EuclidianSpectrum_1D<T>::mul_inv(FourierMapType& m) const
{ {
EuclidianFourierMapComplex<T>& m_c = dynamic_cast<EuclidianFourierMapComplex<T>&>(m); EuclidianFourierMapComplex<T>& m_c = dynamic_cast<EuclidianFourierMapComplex<T>&>(m);
std::complex<T> *d = m.data(); std::complex<T> *d = m.data();
@ -488,7 +483,7 @@ namespace CosmoTool
} }
template<typename T> template<typename T>
void EuclidianSpectrum_1D<T>::mul_inv_sqrt(FourierMap<std::complex<T> >& m) const void EuclidianSpectrum_1D<T>::mul_inv_sqrt(FourierMapType& m) const
{ {
EuclidianFourierMapComplex<T>& m_c = dynamic_cast<EuclidianFourierMapComplex<T>&>(m); EuclidianFourierMapComplex<T>& m_c = dynamic_cast<EuclidianFourierMapComplex<T>&>(m);
std::complex<T> *d = m.data(); std::complex<T> *d = m.data();

12
src/fourier/healpix.hpp
View File

@ -32,7 +32,7 @@ namespace CosmoTool
const T *data() const { return &cls[0]; } const T *data() const { return &cls[0]; }
long size() const { return cls.size(); } long size() const { return cls.size(); }
ptr_map newRandomFourier(gsl_rng *rng, const FourierMapType& like_map) const; void newRandomFourier(gsl_rng *rng, FourierMapType& like_map) const;
void mul(FourierMapType& m) const; void mul(FourierMapType& m) const;
void mul_sqrt(FourierMapType& m) const; void mul_sqrt(FourierMapType& m) const;
@ -244,14 +244,11 @@ namespace CosmoTool
}; };
template<typename T> template<typename T>
typename HealpixSpectrum<T>::ptr_map void HealpixSpectrum<T>::newRandomFourier(gsl_rng *rng, FourierMapType& out_map) const
HealpixSpectrum<T>::newRandomFourier(gsl_rng *rng, const FourierMapType& like_map) const
{ {
const HealpixFourierALM<T>& alms = dynamic_cast<const HealpixFourierALM<T>&>(like_map); HealpixFourierALM<T>& alms = dynamic_cast<HealpixFourierALM<T>&>(out_map);
HealpixFourierALM<T> *new_alms;
ptr_map r(new_alms = new HealpixFourierALM<T>(alms.Lmax(), alms.Mmax()));
long lmaxGen = std::min(cls.size()-1, alms.Lmax()); long lmaxGen = std::min(cls.size()-1, alms.Lmax());
std::complex<T> *new_data = new_alms->data(); std::complex<T> *new_data = alms.data();
for (long l = 0; l <= lmaxGen; l++) for (long l = 0; l <= lmaxGen; l++)
{ {
@ -266,7 +263,6 @@ namespace CosmoTool
c.imag() = gsl_ran_gaussian(rng, Al); c.imag() = gsl_ran_gaussian(rng, Al);
} }
} }
return r;
} }
template<typename T> template<typename T>