diff --git a/CMakeLists.txt b/CMakeLists.txt index 77bd8b8..f9df817 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -55,6 +55,7 @@ include(FindHDF5) include(FindPkgConfig) pkg_check_modules(FFTW3 fftw3>=3.3) +pkg_check_modules(FFTW3F fftw3f>=3.3) pkg_check_modules(EIGEN3 eigen3) include(FindPackageHandleStandardArgs) diff --git a/build_tools/gather_sources.py b/build_tools/gather_sources.py new file mode 100644 index 0000000..ca8c3d1 --- /dev/null +++ b/build_tools/gather_sources.py @@ -0,0 +1,161 @@ +#+ +# This is CosmoTool (./build_tools/gather_sources.py) -- Copyright (C) Guilhem Lavaux (2007-2013) +# +# guilhem.lavaux@gmail.com +# +# This software is a computer program whose purpose is to provide a toolbox for cosmological +# data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...) +# +# This software is governed by the CeCILL license under French law and +# abiding by the rules of distribution of free software. You can use, +# modify and/ or redistribute the software under the terms of the CeCILL +# license as circulated by CEA, CNRS and INRIA at the following URL +# "http://www.cecill.info". +# +# As a counterpart to the access to the source code and rights to copy, +# modify and redistribute granted by the license, users are provided only +# with a limited warranty and the software's author, the holder of the +# economic rights, and the successive licensors have only limited +# liability. +# +# In this respect, the user's attention is drawn to the risks associated +# with loading, using, modifying and/or developing or reproducing the +# software by the user in light of its specific status of free software, +# that may mean that it is complicated to manipulate, and that also +# therefore means that it is reserved for developers and experienced +# professionals having in-depth computer knowledge. Users are therefore +# encouraged to load and test the software's suitability as regards their +# requirements in conditions enabling the security of their systems and/or +# data to be ensured and, more generally, to use and operate it in the +# same conditions as regards security. +# +# The fact that you are presently reading this means that you have had +# knowledge of the CeCILL license and that you accept its terms. +#+ +import shutil +import tempfile +import re +from git import Repo,Tree,Blob + +def apply_license(license, relimit, filename): + header = re.sub(r'@FILENAME@', filename, license) + + f = file(filename) + lines = f.read() + f.close() + + lines = re.sub(relimit, '', lines) + lines = header + lines + + with tempfile.NamedTemporaryFile(delete=False) as tmp_sources: + tmp_sources.write(lines) + + shutil.move(tmp_sources.name, filename) + + +def apply_python_license(filename): + license="""#+ +# This is CosmoTool (@FILENAME@) -- Copyright (C) Guilhem Lavaux (2007-2013) +# +# guilhem.lavaux@gmail.com +# +# This software is a computer program whose purpose is to provide a toolbox for cosmological +# data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...) +# +# This software is governed by the CeCILL license under French law and +# abiding by the rules of distribution of free software. You can use, +# modify and/ or redistribute the software under the terms of the CeCILL +# license as circulated by CEA, CNRS and INRIA at the following URL +# "http://www.cecill.info". +# +# As a counterpart to the access to the source code and rights to copy, +# modify and redistribute granted by the license, users are provided only +# with a limited warranty and the software's author, the holder of the +# economic rights, and the successive licensors have only limited +# liability. +# +# In this respect, the user's attention is drawn to the risks associated +# with loading, using, modifying and/or developing or reproducing the +# software by the user in light of its specific status of free software, +# that may mean that it is complicated to manipulate, and that also +# therefore means that it is reserved for developers and experienced +# professionals having in-depth computer knowledge. Users are therefore +# encouraged to load and test the software's suitability as regards their +# requirements in conditions enabling the security of their systems and/or +# data to be ensured and, more generally, to use and operate it in the +# same conditions as regards security. +# +# The fact that you are presently reading this means that you have had +# knowledge of the CeCILL license and that you accept its terms. +#+ +""" + + print("Shell/Python file: %s" % filename) + relimit=r'^#\+\n(#.*\n)*#\+\n' + apply_license(license, relimit, filename) + + +def apply_cpp_license(filename): + license="""/*+ +This is CosmoTool (@FILENAME@) -- Copyright (C) Guilhem Lavaux (2007-2013) + +guilhem.lavaux@gmail.com + +This software is a computer program whose purpose is to provide a toolbox for cosmological +data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...) + +This software is governed by the CeCILL license under French law and +abiding by the rules of distribution of free software. You can use, +modify and/ or redistribute the software under the terms of the CeCILL +license as circulated by CEA, CNRS and INRIA at the following URL +"http://www.cecill.info". + +As a counterpart to the access to the source code and rights to copy, +modify and redistribute granted by the license, users are provided only +with a limited warranty and the software's author, the holder of the +economic rights, and the successive licensors have only limited +liability. + +In this respect, the user's attention is drawn to the risks associated +with loading, using, modifying and/or developing or reproducing the +software by the user in light of its specific status of free software, +that may mean that it is complicated to manipulate, and that also +therefore means that it is reserved for developers and experienced +professionals having in-depth computer knowledge. Users are therefore +encouraged to load and test the software's suitability as regards their +requirements in conditions enabling the security of their systems and/or +data to be ensured and, more generally, to use and operate it in the +same conditions as regards security. + +The fact that you are presently reading this means that you have had +knowledge of the CeCILL license and that you accept its terms. ++*/ +""" + relimit = r'^(?s)/\*\+.*\+\*/' + print("C++ file: %s" % filename) + apply_license(license, relimit, filename) + + +def analyze_tree(prefix, t): + for ename,entry in t.items(): + if ename == 'external' or ename == 'zobov': + continue + if type(entry) == Tree: + analyze_tree(prefix + "/" + ename, entry) + elif type(entry) == Blob: + if ename == './src/hdf5_flash.h' or ename == './src/h5_readFlash.cpp' or ename == './src/h5_readFlash.hpp': + continue + + if re.match(".*\.(sh|py|pyx)$",ename) != None: + fname=prefix+"/"+ename + apply_python_license(fname) + if re.match('.*\.(cpp|hpp|h)$', ename) != None: + fname=prefix+"/"+ename + apply_cpp_license(fname) + + +if __name__=="__main__": + repo = Repo(".") + assert repo.bare == False + t = repo.tree() + analyze_tree(".", t) diff --git a/sample/CMakeLists.txt b/sample/CMakeLists.txt index 816d7bb..d158ae2 100644 --- a/sample/CMakeLists.txt +++ b/sample/CMakeLists.txt @@ -20,6 +20,9 @@ target_link_libraries(testkd ${tolink}) add_executable(testkd2 testkd2.cpp) target_link_libraries(testkd2 ${tolink}) +add_executable(testkd3 testkd3.cpp) +target_link_libraries(testkd3 ${tolink}) + add_executable(testDelaunay testDelaunay.cpp) target_link_libraries(testDelaunay ${tolink}) @@ -44,10 +47,10 @@ target_link_libraries(testAlgo ${tolink}) add_executable(testBSP testBSP.cpp) target_link_libraries(testBSP ${tolink}) -if (FFTW3_FOUND AND EIGEN3_FOUND) +if (FFTW3_FOUND AND FFTW3F_FOUND AND EIGEN3_FOUND) add_executable(test_fft_calls test_fft_calls.cpp) - target_link_libraries(test_fft_calls ${tolink} ${FFTW3_LIBRARIES}) -endif (FFTW3_FOUND AND EIGEN3_FOUND) + target_link_libraries(test_fft_calls ${tolink} ${FFTW3_LIBRARIES} ${FFTW3F_LIBRARIES}) +endif (FFTW3_FOUND AND FFTW3F_FOUND AND EIGEN3_FOUND) if (SHARP_LIBRARY AND SHARP_INCLUDE_PATH AND EIGEN3_FOUND) include_directories(${SHARP_INCLUDE_PATH}) @@ -55,3 +58,7 @@ if (SHARP_LIBRARY AND SHARP_INCLUDE_PATH AND EIGEN3_FOUND) target_link_libraries(test_healpix_calls ${tolink} ${SHARP_LIBRARIES}) set_target_properties(test_healpix_calls PROPERTIES COMPILE_FLAGS ${OpenMP_CXX_FLAGS} LINK_FLAGS ${OpenMP_CXX_FLAGS}) endif (SHARP_LIBRARY AND SHARP_INCLUDE_PATH AND EIGEN3_FOUND) + +add_executable(test_cosmopower test_cosmopower.cpp) +target_link_libraries(test_cosmopower ${tolink}) + diff --git a/sample/testAlgo.cpp b/sample/testAlgo.cpp index 65614f5..cdf7971 100644 --- a/sample/testAlgo.cpp +++ b/sample/testAlgo.cpp @@ -1,3 +1,38 @@ +/*+ +This is CosmoTool (./sample/testAlgo.cpp) -- Copyright (C) Guilhem Lavaux (2007-2013) + +guilhem.lavaux@gmail.com + +This software is a computer program whose purpose is to provide a toolbox for cosmological +data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...) + +This software is governed by the CeCILL license under French law and +abiding by the rules of distribution of free software. You can use, +modify and/ or redistribute the software under the terms of the CeCILL +license as circulated by CEA, CNRS and INRIA at the following URL +"http://www.cecill.info". + +As a counterpart to the access to the source code and rights to copy, +modify and redistribute granted by the license, users are provided only +with a limited warranty and the software's author, the holder of the +economic rights, and the successive licensors have only limited +liability. + +In this respect, the user's attention is drawn to the risks associated +with loading, using, modifying and/or developing or reproducing the +software by the user in light of its specific status of free software, +that may mean that it is complicated to manipulate, and that also +therefore means that it is reserved for developers and experienced +professionals having in-depth computer knowledge. Users are therefore +encouraged to load and test the software's suitability as regards their +requirements in conditions enabling the security of their systems and/or +data to be ensured and, more generally, to use and operate it in the +same conditions as regards security. + +The fact that you are presently reading this means that you have had +knowledge of the CeCILL license and that you accept its terms. ++*/ + #include #include #include "algo.hpp" diff --git a/sample/testBQueue.cpp b/sample/testBQueue.cpp index 822cd6b..08beb6d 100644 --- a/sample/testBQueue.cpp +++ b/sample/testBQueue.cpp @@ -1,3 +1,38 @@ +/*+ +This is CosmoTool (./sample/testBQueue.cpp) -- Copyright (C) Guilhem Lavaux (2007-2013) + +guilhem.lavaux@gmail.com + +This software is a computer program whose purpose is to provide a toolbox for cosmological +data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...) + +This software is governed by the CeCILL license under French law and +abiding by the rules of distribution of free software. You can use, +modify and/ or redistribute the software under the terms of the CeCILL +license as circulated by CEA, CNRS and INRIA at the following URL +"http://www.cecill.info". + +As a counterpart to the access to the source code and rights to copy, +modify and redistribute granted by the license, users are provided only +with a limited warranty and the software's author, the holder of the +economic rights, and the successive licensors have only limited +liability. + +In this respect, the user's attention is drawn to the risks associated +with loading, using, modifying and/or developing or reproducing the +software by the user in light of its specific status of free software, +that may mean that it is complicated to manipulate, and that also +therefore means that it is reserved for developers and experienced +professionals having in-depth computer knowledge. Users are therefore +encouraged to load and test the software's suitability as regards their +requirements in conditions enabling the security of their systems and/or +data to be ensured and, more generally, to use and operate it in the +same conditions as regards security. + +The fact that you are presently reading this means that you have had +knowledge of the CeCILL license and that you accept its terms. ++*/ + #include #include "bqueue.hpp" diff --git a/sample/testBSP.cpp b/sample/testBSP.cpp index f142fc2..def0717 100644 --- a/sample/testBSP.cpp +++ b/sample/testBSP.cpp @@ -1,3 +1,38 @@ +/*+ +This is CosmoTool (./sample/testBSP.cpp) -- Copyright (C) Guilhem Lavaux (2007-2013) + +guilhem.lavaux@gmail.com + +This software is a computer program whose purpose is to provide a toolbox for cosmological +data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...) + +This software is governed by the CeCILL license under French law and +abiding by the rules of distribution of free software. You can use, +modify and/ or redistribute the software under the terms of the CeCILL +license as circulated by CEA, CNRS and INRIA at the following URL +"http://www.cecill.info". + +As a counterpart to the access to the source code and rights to copy, +modify and redistribute granted by the license, users are provided only +with a limited warranty and the software's author, the holder of the +economic rights, and the successive licensors have only limited +liability. + +In this respect, the user's attention is drawn to the risks associated +with loading, using, modifying and/or developing or reproducing the +software by the user in light of its specific status of free software, +that may mean that it is complicated to manipulate, and that also +therefore means that it is reserved for developers and experienced +professionals having in-depth computer knowledge. Users are therefore +encouraged to load and test the software's suitability as regards their +requirements in conditions enabling the security of their systems and/or +data to be ensured and, more generally, to use and operate it in the +same conditions as regards security. + +The fact that you are presently reading this means that you have had +knowledge of the CeCILL license and that you accept its terms. ++*/ + #if 0 #include "bsp_simple.hpp" diff --git a/sample/testDelaunay.cpp b/sample/testDelaunay.cpp index fe14553..8d2d9d4 100644 --- a/sample/testDelaunay.cpp +++ b/sample/testDelaunay.cpp @@ -1,3 +1,38 @@ +/*+ +This is CosmoTool (./sample/testDelaunay.cpp) -- Copyright (C) Guilhem Lavaux (2007-2013) + +guilhem.lavaux@gmail.com + +This software is a computer program whose purpose is to provide a toolbox for cosmological +data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...) + +This software is governed by the CeCILL license under French law and +abiding by the rules of distribution of free software. You can use, +modify and/ or redistribute the software under the terms of the CeCILL +license as circulated by CEA, CNRS and INRIA at the following URL +"http://www.cecill.info". + +As a counterpart to the access to the source code and rights to copy, +modify and redistribute granted by the license, users are provided only +with a limited warranty and the software's author, the holder of the +economic rights, and the successive licensors have only limited +liability. + +In this respect, the user's attention is drawn to the risks associated +with loading, using, modifying and/or developing or reproducing the +software by the user in light of its specific status of free software, +that may mean that it is complicated to manipulate, and that also +therefore means that it is reserved for developers and experienced +professionals having in-depth computer knowledge. Users are therefore +encouraged to load and test the software's suitability as regards their +requirements in conditions enabling the security of their systems and/or +data to be ensured and, more generally, to use and operate it in the +same conditions as regards security. + +The fact that you are presently reading this means that you have had +knowledge of the CeCILL license and that you accept its terms. ++*/ + #include #include #include "dinterpolate.hpp" diff --git a/sample/testEskow.cpp b/sample/testEskow.cpp index 60c1821..aa48d2b 100644 --- a/sample/testEskow.cpp +++ b/sample/testEskow.cpp @@ -1,3 +1,38 @@ +/*+ +This is CosmoTool (./sample/testEskow.cpp) -- Copyright (C) Guilhem Lavaux (2007-2013) + +guilhem.lavaux@gmail.com + +This software is a computer program whose purpose is to provide a toolbox for cosmological +data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...) + +This software is governed by the CeCILL license under French law and +abiding by the rules of distribution of free software. You can use, +modify and/ or redistribute the software under the terms of the CeCILL +license as circulated by CEA, CNRS and INRIA at the following URL +"http://www.cecill.info". + +As a counterpart to the access to the source code and rights to copy, +modify and redistribute granted by the license, users are provided only +with a limited warranty and the software's author, the holder of the +economic rights, and the successive licensors have only limited +liability. + +In this respect, the user's attention is drawn to the risks associated +with loading, using, modifying and/or developing or reproducing the +software by the user in light of its specific status of free software, +that may mean that it is complicated to manipulate, and that also +therefore means that it is reserved for developers and experienced +professionals having in-depth computer knowledge. Users are therefore +encouraged to load and test the software's suitability as regards their +requirements in conditions enabling the security of their systems and/or +data to be ensured and, more generally, to use and operate it in the +same conditions as regards security. + +The fact that you are presently reading this means that you have had +knowledge of the CeCILL license and that you accept its terms. ++*/ + #include #include #include diff --git a/sample/testInterpolate.cpp b/sample/testInterpolate.cpp index df146b2..b863170 100644 --- a/sample/testInterpolate.cpp +++ b/sample/testInterpolate.cpp @@ -1,3 +1,38 @@ +/*+ +This is CosmoTool (./sample/testInterpolate.cpp) -- Copyright (C) Guilhem Lavaux (2007-2013) + +guilhem.lavaux@gmail.com + +This software is a computer program whose purpose is to provide a toolbox for cosmological +data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...) + +This software is governed by the CeCILL license under French law and +abiding by the rules of distribution of free software. You can use, +modify and/ or redistribute the software under the terms of the CeCILL +license as circulated by CEA, CNRS and INRIA at the following URL +"http://www.cecill.info". + +As a counterpart to the access to the source code and rights to copy, +modify and redistribute granted by the license, users are provided only +with a limited warranty and the software's author, the holder of the +economic rights, and the successive licensors have only limited +liability. + +In this respect, the user's attention is drawn to the risks associated +with loading, using, modifying and/or developing or reproducing the +software by the user in light of its specific status of free software, +that may mean that it is complicated to manipulate, and that also +therefore means that it is reserved for developers and experienced +professionals having in-depth computer knowledge. Users are therefore +encouraged to load and test the software's suitability as regards their +requirements in conditions enabling the security of their systems and/or +data to be ensured and, more generally, to use and operate it in the +same conditions as regards security. + +The fact that you are presently reading this means that you have had +knowledge of the CeCILL license and that you accept its terms. ++*/ + #include #include "interpolate3d.hpp" diff --git a/sample/testNewton.cpp b/sample/testNewton.cpp index 27415b8..5677bbc 100644 --- a/sample/testNewton.cpp +++ b/sample/testNewton.cpp @@ -1,3 +1,38 @@ +/*+ +This is CosmoTool (./sample/testNewton.cpp) -- Copyright (C) Guilhem Lavaux (2007-2013) + +guilhem.lavaux@gmail.com + +This software is a computer program whose purpose is to provide a toolbox for cosmological +data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...) + +This software is governed by the CeCILL license under French law and +abiding by the rules of distribution of free software. You can use, +modify and/ or redistribute the software under the terms of the CeCILL +license as circulated by CEA, CNRS and INRIA at the following URL +"http://www.cecill.info". + +As a counterpart to the access to the source code and rights to copy, +modify and redistribute granted by the license, users are provided only +with a limited warranty and the software's author, the holder of the +economic rights, and the successive licensors have only limited +liability. + +In this respect, the user's attention is drawn to the risks associated +with loading, using, modifying and/or developing or reproducing the +software by the user in light of its specific status of free software, +that may mean that it is complicated to manipulate, and that also +therefore means that it is reserved for developers and experienced +professionals having in-depth computer knowledge. Users are therefore +encouraged to load and test the software's suitability as regards their +requirements in conditions enabling the security of their systems and/or +data to be ensured and, more generally, to use and operate it in the +same conditions as regards security. + +The fact that you are presently reading this means that you have had +knowledge of the CeCILL license and that you accept its terms. ++*/ + #include #include diff --git a/sample/testPool.cpp b/sample/testPool.cpp index ee01468..e4ff8c1 100644 --- a/sample/testPool.cpp +++ b/sample/testPool.cpp @@ -1,3 +1,38 @@ +/*+ +This is CosmoTool (./sample/testPool.cpp) -- Copyright (C) Guilhem Lavaux (2007-2013) + +guilhem.lavaux@gmail.com + +This software is a computer program whose purpose is to provide a toolbox for cosmological +data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...) + +This software is governed by the CeCILL license under French law and +abiding by the rules of distribution of free software. You can use, +modify and/ or redistribute the software under the terms of the CeCILL +license as circulated by CEA, CNRS and INRIA at the following URL +"http://www.cecill.info". + +As a counterpart to the access to the source code and rights to copy, +modify and redistribute granted by the license, users are provided only +with a limited warranty and the software's author, the holder of the +economic rights, and the successive licensors have only limited +liability. + +In this respect, the user's attention is drawn to the risks associated +with loading, using, modifying and/or developing or reproducing the +software by the user in light of its specific status of free software, +that may mean that it is complicated to manipulate, and that also +therefore means that it is reserved for developers and experienced +professionals having in-depth computer knowledge. Users are therefore +encouraged to load and test the software's suitability as regards their +requirements in conditions enabling the security of their systems and/or +data to be ensured and, more generally, to use and operate it in the +same conditions as regards security. + +The fact that you are presently reading this means that you have had +knowledge of the CeCILL license and that you accept its terms. ++*/ + #include "pool.hpp" using namespace CosmoTool; diff --git a/sample/testReadFlash.cpp b/sample/testReadFlash.cpp index 61d038e..93a5b64 100644 --- a/sample/testReadFlash.cpp +++ b/sample/testReadFlash.cpp @@ -1,3 +1,38 @@ +/*+ +This is CosmoTool (./sample/testReadFlash.cpp) -- Copyright (C) Guilhem Lavaux (2007-2013) + +guilhem.lavaux@gmail.com + +This software is a computer program whose purpose is to provide a toolbox for cosmological +data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...) + +This software is governed by the CeCILL license under French law and +abiding by the rules of distribution of free software. You can use, +modify and/ or redistribute the software under the terms of the CeCILL +license as circulated by CEA, CNRS and INRIA at the following URL +"http://www.cecill.info". + +As a counterpart to the access to the source code and rights to copy, +modify and redistribute granted by the license, users are provided only +with a limited warranty and the software's author, the holder of the +economic rights, and the successive licensors have only limited +liability. + +In this respect, the user's attention is drawn to the risks associated +with loading, using, modifying and/or developing or reproducing the +software by the user in light of its specific status of free software, +that may mean that it is complicated to manipulate, and that also +therefore means that it is reserved for developers and experienced +professionals having in-depth computer knowledge. Users are therefore +encouraged to load and test the software's suitability as regards their +requirements in conditions enabling the security of their systems and/or +data to be ensured and, more generally, to use and operate it in the +same conditions as regards security. + +The fact that you are presently reading this means that you have had +knowledge of the CeCILL license and that you accept its terms. ++*/ + #include #include #include diff --git a/sample/testSmooth.cpp b/sample/testSmooth.cpp index d1c94a9..4f3011d 100644 --- a/sample/testSmooth.cpp +++ b/sample/testSmooth.cpp @@ -1,3 +1,38 @@ +/*+ +This is CosmoTool (./sample/testSmooth.cpp) -- Copyright (C) Guilhem Lavaux (2007-2013) + +guilhem.lavaux@gmail.com + +This software is a computer program whose purpose is to provide a toolbox for cosmological +data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...) + +This software is governed by the CeCILL license under French law and +abiding by the rules of distribution of free software. You can use, +modify and/ or redistribute the software under the terms of the CeCILL +license as circulated by CEA, CNRS and INRIA at the following URL +"http://www.cecill.info". + +As a counterpart to the access to the source code and rights to copy, +modify and redistribute granted by the license, users are provided only +with a limited warranty and the software's author, the holder of the +economic rights, and the successive licensors have only limited +liability. + +In this respect, the user's attention is drawn to the risks associated +with loading, using, modifying and/or developing or reproducing the +software by the user in light of its specific status of free software, +that may mean that it is complicated to manipulate, and that also +therefore means that it is reserved for developers and experienced +professionals having in-depth computer knowledge. Users are therefore +encouraged to load and test the software's suitability as regards their +requirements in conditions enabling the security of their systems and/or +data to be ensured and, more generally, to use and operate it in the +same conditions as regards security. + +The fact that you are presently reading this means that you have had +knowledge of the CeCILL license and that you accept its terms. ++*/ + #include #include "sphSmooth.hpp" #include "yorick.hpp" diff --git a/sample/test_cosmopower.cpp b/sample/test_cosmopower.cpp new file mode 100644 index 0000000..1bf1ec8 --- /dev/null +++ b/sample/test_cosmopower.cpp @@ -0,0 +1,63 @@ +/*+ +This is CosmoTool (./sample/test_cosmopower.cpp) -- Copyright (C) Guilhem Lavaux (2007-2013) + +guilhem.lavaux@gmail.com + +This software is a computer program whose purpose is to provide a toolbox for cosmological +data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...) + +This software is governed by the CeCILL license under French law and +abiding by the rules of distribution of free software. You can use, +modify and/ or redistribute the software under the terms of the CeCILL +license as circulated by CEA, CNRS and INRIA at the following URL +"http://www.cecill.info". + +As a counterpart to the access to the source code and rights to copy, +modify and redistribute granted by the license, users are provided only +with a limited warranty and the software's author, the holder of the +economic rights, and the successive licensors have only limited +liability. + +In this respect, the user's attention is drawn to the risks associated +with loading, using, modifying and/or developing or reproducing the +software by the user in light of its specific status of free software, +that may mean that it is complicated to manipulate, and that also +therefore means that it is reserved for developers and experienced +professionals having in-depth computer knowledge. Users are therefore +encouraged to load and test the software's suitability as regards their +requirements in conditions enabling the security of their systems and/or +data to be ensured and, more generally, to use and operate it in the +same conditions as regards security. + +The fact that you are presently reading this means that you have had +knowledge of the CeCILL license and that you accept its terms. ++*/ +#include +#include +#include "cosmopower.hpp" + +using namespace std; +using namespace CosmoTool; + +int main() +{ + CosmoPower cp; + CosmoPower::CosmoFunction f[] = { CosmoPower::POWER_EFSTATHIOU, CosmoPower::HU_WIGGLES, CosmoPower::HU_BARYON, CosmoPower::POWER_SUGIYAMA }; + int num_F = sizeof(f)/sizeof(f[0]); + + cp.setFunction(f[0]); + cp.normalize(); + for (int ik = 0; ik < 100; ik++) + { + double k = pow(10.0, 4*ik/100.-3); + + cout << k << " "; + for (int q = 0; q < num_F; q++) + { + cp.setFunction(f[q]); + cout << cp.power(k) << " "; + } + cout << endl; + } + return 0; +} diff --git a/sample/test_fft_calls.cpp b/sample/test_fft_calls.cpp index 06c1c9a..b98fcea 100644 --- a/sample/test_fft_calls.cpp +++ b/sample/test_fft_calls.cpp @@ -1,20 +1,93 @@ +/*+ +This is CosmoTool (./sample/test_fft_calls.cpp) -- Copyright (C) Guilhem Lavaux (2007-2013) + +guilhem.lavaux@gmail.com + +This software is a computer program whose purpose is to provide a toolbox for cosmological +data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...) + +This software is governed by the CeCILL license under French law and +abiding by the rules of distribution of free software. You can use, +modify and/ or redistribute the software under the terms of the CeCILL +license as circulated by CEA, CNRS and INRIA at the following URL +"http://www.cecill.info". + +As a counterpart to the access to the source code and rights to copy, +modify and redistribute granted by the license, users are provided only +with a limited warranty and the software's author, the holder of the +economic rights, and the successive licensors have only limited +liability. + +In this respect, the user's attention is drawn to the risks associated +with loading, using, modifying and/or developing or reproducing the +software by the user in light of its specific status of free software, +that may mean that it is complicated to manipulate, and that also +therefore means that it is reserved for developers and experienced +professionals having in-depth computer knowledge. Users are therefore +encouraged to load and test the software's suitability as regards their +requirements in conditions enabling the security of their systems and/or +data to be ensured and, more generally, to use and operate it in the +same conditions as regards security. + +The fact that you are presently reading this means that you have had +knowledge of the CeCILL license and that you accept its terms. ++*/ + +#include +#include "yorick.hpp" +#include #include +#include #include "fourier/euclidian.hpp" using namespace CosmoTool; +using boost::format; +using boost::str; using namespace std; -int main() +double spectrum_generator(double k) { - EuclidianFourierTransform_2d dft(128,128,1.0,1.0); - double volume = 128*128; + if (k==0) + return 0; + return 1/(0.01+pow(k, 3.0)); +} + +template +void test_2d(int Nx, int Ny) +{ + EuclidianFourierTransform_2d dft(Nx,Ny,1.0,1.0); + EuclidianSpectrum_1D spectrum(spectrum_generator); + double volume = Nx*Ny; + gsl_rng *rng = gsl_rng_alloc(gsl_rng_default); dft.realSpace().eigen().setRandom(); + dft.fourierSpace().scale(std::complex(0,0)); dft.analysis(); + cout << format("Testing (%d,%d)") % Nx % Ny << endl; cout << "Map dot-product = " << dft.realSpace().dot_product(dft.realSpace()) << endl; cout << "Fourier dot-product = " << dft.fourierSpace().dot_product(dft.fourierSpace()).real()*volume << endl; dft.synthesis(); cout << "Resynthesis dot-product = " << dft.realSpace().dot_product(dft.realSpace()) << endl; - return 0; + dft.realSpace().scale(2.0); + dft.fourierSpace().scale(0.2); + + spectrum.newRandomFourier(rng, dft.fourierSpace()); + dft.synthesis(); + + uint32_t dims[2] = { Ny, Nx }; + CosmoTool::saveArray(str(format("generated_map_%d_%d.nc") %Nx % Ny) , dft.realSpace().data(), dims, 2); + + gsl_rng_free(rng); +} + +int main(int argc, char **argv) +{ + test_2d(128,128); + test_2d(131,128); + test_2d(130,128); + test_2d(128,128); + test_2d(128,131); + test_2d(128,130); + return 0; } diff --git a/sample/test_healpix_calls.cpp b/sample/test_healpix_calls.cpp index ae1844e..75473d8 100644 --- a/sample/test_healpix_calls.cpp +++ b/sample/test_healpix_calls.cpp @@ -1,3 +1,38 @@ +/*+ +This is CosmoTool (./sample/test_healpix_calls.cpp) -- Copyright (C) Guilhem Lavaux (2007-2013) + +guilhem.lavaux@gmail.com + +This software is a computer program whose purpose is to provide a toolbox for cosmological +data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...) + +This software is governed by the CeCILL license under French law and +abiding by the rules of distribution of free software. You can use, +modify and/ or redistribute the software under the terms of the CeCILL +license as circulated by CEA, CNRS and INRIA at the following URL +"http://www.cecill.info". + +As a counterpart to the access to the source code and rights to copy, +modify and redistribute granted by the license, users are provided only +with a limited warranty and the software's author, the holder of the +economic rights, and the successive licensors have only limited +liability. + +In this respect, the user's attention is drawn to the risks associated +with loading, using, modifying and/or developing or reproducing the +software by the user in light of its specific status of free software, +that may mean that it is complicated to manipulate, and that also +therefore means that it is reserved for developers and experienced +professionals having in-depth computer knowledge. Users are therefore +encouraged to load and test the software's suitability as regards their +requirements in conditions enabling the security of their systems and/or +data to be ensured and, more generally, to use and operate it in the +same conditions as regards security. + +The fact that you are presently reading this means that you have had +knowledge of the CeCILL license and that you accept its terms. ++*/ + #include #include "fourier/healpix.hpp" @@ -6,15 +41,20 @@ using namespace std; int main() { - HealpixFourierTransform dft(128,3*128,3*128, 40); + HealpixFourierTransform dft(128,2*128,2*128, 40); + HealpixUtilityFunction utils; long Npix = dft.realSpace().size(); dft.realSpace().eigen().setRandom(); dft.analysis(); cout << "Map dot-product = " << dft.realSpace().dot_product(dft.realSpace()) << endl; cout << "Fourier dot-product = " << dft.fourierSpace().dot_product(dft.fourierSpace()).real()*Npix/(4*M_PI) << endl; + cout << "Spectrum analysis" << endl; + HealpixUtilityFunction::Spectrum_ptr s = utils.estimateSpectrumFromMap(dft.fourierSpace()); + dft.synthesis(); cout << "Resynthesis dot-product = " << dft.realSpace().dot_product(dft.realSpace()) << endl; + return 0; } diff --git a/sample/testkd.cpp b/sample/testkd.cpp index eabcfb2..f810571 100644 --- a/sample/testkd.cpp +++ b/sample/testkd.cpp @@ -1,3 +1,38 @@ +/*+ +This is CosmoTool (./sample/testkd.cpp) -- Copyright (C) Guilhem Lavaux (2007-2013) + +guilhem.lavaux@gmail.com + +This software is a computer program whose purpose is to provide a toolbox for cosmological +data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...) + +This software is governed by the CeCILL license under French law and +abiding by the rules of distribution of free software. You can use, +modify and/ or redistribute the software under the terms of the CeCILL +license as circulated by CEA, CNRS and INRIA at the following URL +"http://www.cecill.info". + +As a counterpart to the access to the source code and rights to copy, +modify and redistribute granted by the license, users are provided only +with a limited warranty and the software's author, the holder of the +economic rights, and the successive licensors have only limited +liability. + +In this respect, the user's attention is drawn to the risks associated +with loading, using, modifying and/or developing or reproducing the +software by the user in light of its specific status of free software, +that may mean that it is complicated to manipulate, and that also +therefore means that it is reserved for developers and experienced +professionals having in-depth computer knowledge. Users are therefore +encouraged to load and test the software's suitability as regards their +requirements in conditions enabling the security of their systems and/or +data to be ensured and, more generally, to use and operate it in the +same conditions as regards security. + +The fact that you are presently reading this means that you have had +knowledge of the CeCILL license and that you accept its terms. ++*/ + #define __KD_TREE_SAVE_ON_DISK #include #include diff --git a/sample/testkd2.cpp b/sample/testkd2.cpp index bc65899..7c29b1e 100644 --- a/sample/testkd2.cpp +++ b/sample/testkd2.cpp @@ -1,3 +1,38 @@ +/*+ +This is CosmoTool (./sample/testkd2.cpp) -- Copyright (C) Guilhem Lavaux (2007-2013) + +guilhem.lavaux@gmail.com + +This software is a computer program whose purpose is to provide a toolbox for cosmological +data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...) + +This software is governed by the CeCILL license under French law and +abiding by the rules of distribution of free software. You can use, +modify and/ or redistribute the software under the terms of the CeCILL +license as circulated by CEA, CNRS and INRIA at the following URL +"http://www.cecill.info". + +As a counterpart to the access to the source code and rights to copy, +modify and redistribute granted by the license, users are provided only +with a limited warranty and the software's author, the holder of the +economic rights, and the successive licensors have only limited +liability. + +In this respect, the user's attention is drawn to the risks associated +with loading, using, modifying and/or developing or reproducing the +software by the user in light of its specific status of free software, +that may mean that it is complicated to manipulate, and that also +therefore means that it is reserved for developers and experienced +professionals having in-depth computer knowledge. Users are therefore +encouraged to load and test the software's suitability as regards their +requirements in conditions enabling the security of their systems and/or +data to be ensured and, more generally, to use and operate it in the +same conditions as regards security. + +The fact that you are presently reading this means that you have had +knowledge of the CeCILL license and that you accept its terms. ++*/ + #include #include #include diff --git a/sample/testkd3.cpp b/sample/testkd3.cpp new file mode 100644 index 0000000..50eba07 --- /dev/null +++ b/sample/testkd3.cpp @@ -0,0 +1,187 @@ +/*+ +This is CosmoTool (./sample/testkd3.cpp) -- Copyright (C) Guilhem Lavaux (2007-2013) + +guilhem.lavaux@gmail.com + +This software is a computer program whose purpose is to provide a toolbox for cosmological +data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...) + +This software is governed by the CeCILL license under French law and +abiding by the rules of distribution of free software. You can use, +modify and/ or redistribute the software under the terms of the CeCILL +license as circulated by CEA, CNRS and INRIA at the following URL +"http://www.cecill.info". + +As a counterpart to the access to the source code and rights to copy, +modify and redistribute granted by the license, users are provided only +with a limited warranty and the software's author, the holder of the +economic rights, and the successive licensors have only limited +liability. + +In this respect, the user's attention is drawn to the risks associated +with loading, using, modifying and/or developing or reproducing the +software by the user in light of its specific status of free software, +that may mean that it is complicated to manipulate, and that also +therefore means that it is reserved for developers and experienced +professionals having in-depth computer knowledge. Users are therefore +encouraged to load and test the software's suitability as regards their +requirements in conditions enabling the security of their systems and/or +data to be ensured and, more generally, to use and operate it in the +same conditions as regards security. + +The fact that you are presently reading this means that you have had +knowledge of the CeCILL license and that you accept its terms. ++*/ + +#include +#include +#include +#include +#include +#define __KD_TREE_NUMNODES +#include "mykdtree.hpp" +#include "kdtree_splitters.hpp" +#include + +#define NTRY 30 +#define NGB 24 +#define ND 3 + +using namespace std; +using namespace CosmoTool; + +typedef KDTree > MyTree; +//typedef KDTree MyTree; +typedef KDCell MyCell; + +static double periodic(double a) +{ + while (a < -0.5) + a+=1; + while (a > 0.5) + a -= 1; + return a; +} + +MyCell *findNearest(MyTree::coords& xc, MyCell *cells, uint32_t Ncells) +{ + MyCell *near2 = 0; + double R2 = INFINITY; + for (int i = 0; i < Ncells; i++) + { + double d2 = 0; + for (int j = 0; j < ND; j++) + { + double delta = periodic(xc[j]-cells[i].coord[j]); + d2 += delta*delta; + } + if (d2 < R2) + { + near2 = &cells[i]; + R2 = d2; + } + } + return near2; +} + +int main() +{ + uint32_t Ncells = 3000000; + MyCell *cells = new MyCell[Ncells]; + + for (int i = 0; i < Ncells; i++) + { + cells[i].active = true; + for (int l = 0; l < ND; l++) + cells[i].coord[l] = drand48(); + } + + // Check timing + clock_t startTimer = clock(); + MyTree tree(cells, Ncells); + clock_t endTimer = clock(); + + tree.setPeriodic(true, 1.0); + + clock_t delta = endTimer-startTimer; + double myTime = delta*1.0/CLOCKS_PER_SEC * 1.0; + + cout << "KDTree build = " << myTime << " s" << endl; + + MyTree::coords *xc = new MyTree::coords[NTRY]; + + cout << "Generating seeds..." << endl; + for (int k = 0; k < NTRY; k++) + { + for (int l = 0; l < ND; l++) + xc[k][l] = drand48(); + } + + // Check consistency + cout << "Check consistency..." << endl; + MyCell **ngb = new MyCell *[NGB]; + double *distances = new double[NGB]; + + ofstream fngb("nearest.txt"); + for (int k = 0; k < NTRY; k++) { + cout << "Seed = " << xc[k][0] << " " << xc[k][1] << " " << xc[k][2] << endl; + tree.getNearestNeighbours(xc[k], NGB, ngb, distances); + int last = -1; + + for (uint32_t i = 0; i < NGB; i++) + { + if (ngb[i] == 0) + continue; + + last = i; + + double d2 = 0; + for (int l = 0; l < 3; l++) + d2 += ({double delta = periodic(xc[k][l] - ngb[i]->coord[l]); delta*delta;}); + fngb << ngb[i]->coord[0] << " " << ngb[i]->coord[1] << " " << ngb[i]->coord[2] << " " << sqrt(d2) << " " << sqrt(distances[i]) << endl; + } + fngb << endl << endl; + double farther_dist = distances[last]; + for (uint32_t i = 0; i < Ncells; i++) + { + bool found = false; + // If the points is not in the list, it means it is farther than the farthest point + for (int j =0; j < NGB; j++) + { + if (&cells[i] == ngb[j]) { + found = true; + break; + } + } + double dist_to_seed = 0; + for (int l = 0; l < 3; l++) + { + double delta = periodic(xc[k][l]-cells[i].coord[l]); + dist_to_seed += delta*delta; + } + double delta_epsilon = fabs(farther_dist/dist_to_seed-1); + if (!found) + { + if (dist_to_seed <= farther_dist && + delta_epsilon > 1000*sqrt(numeric_limits::epsilon())) + { + cout << "SHOULD BE IN (d=" << dist_to_seed << "): Failed iteration " << k << " particle " << i << " (farthest=" << sqrt(farther_dist) << ")" << endl; + cout << "delta_epsilon = " << delta_epsilon << endl; + abort(); + } + } + else + { + if (dist_to_seed > farther_dist && + delta_epsilon > 1000*sqrt(numeric_limits::epsilon())) + { + cout << "SHOULD BE OUT (d=" << sqrt(dist_to_seed) << "): Failed iteration " << k << " particle " << i << " (farthest=" << sqrt(farther_dist) << ")" << endl; + cout << "delta_epsilon = " << delta_epsilon << " epsilon = " << 100*sqrt(numeric_limits::epsilon()) << endl; + abort(); + } + } + } + } + + return 0; +} diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index c99918e..8901027 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -8,6 +8,7 @@ SET(CosmoTool_SRCS powerSpectrum.cpp miniargs.cpp growthFactor.cpp + cosmopower.cpp ) IF(FOUND_NETCDF3) diff --git a/src/algo.hpp b/src/algo.hpp index 485db75..2dce3eb 100644 --- a/src/algo.hpp +++ b/src/algo.hpp @@ -1,3 +1,38 @@ +/*+ +This is CosmoTool (./src/algo.hpp) -- Copyright (C) Guilhem Lavaux (2007-2013) + +guilhem.lavaux@gmail.com + +This software is a computer program whose purpose is to provide a toolbox for cosmological +data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...) + +This software is governed by the CeCILL license under French law and +abiding by the rules of distribution of free software. You can use, +modify and/ or redistribute the software under the terms of the CeCILL +license as circulated by CEA, CNRS and INRIA at the following URL +"http://www.cecill.info". + +As a counterpart to the access to the source code and rights to copy, +modify and redistribute granted by the license, users are provided only +with a limited warranty and the software's author, the holder of the +economic rights, and the successive licensors have only limited +liability. + +In this respect, the user's attention is drawn to the risks associated +with loading, using, modifying and/or developing or reproducing the +software by the user in light of its specific status of free software, +that may mean that it is complicated to manipulate, and that also +therefore means that it is reserved for developers and experienced +professionals having in-depth computer knowledge. Users are therefore +encouraged to load and test the software's suitability as regards their +requirements in conditions enabling the security of their systems and/or +data to be ensured and, more generally, to use and operate it in the +same conditions as regards security. + +The fact that you are presently reading this means that you have had +knowledge of the CeCILL license and that you accept its terms. ++*/ + #ifndef __COSMOTOOL_ALGO_HPP #define __COSMOTOOL_ALGO_HPP diff --git a/src/bqueue.hpp b/src/bqueue.hpp index 792d309..872830e 100644 --- a/src/bqueue.hpp +++ b/src/bqueue.hpp @@ -1,3 +1,38 @@ +/*+ +This is CosmoTool (./src/bqueue.hpp) -- Copyright (C) Guilhem Lavaux (2007-2013) + +guilhem.lavaux@gmail.com + +This software is a computer program whose purpose is to provide a toolbox for cosmological +data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...) + +This software is governed by the CeCILL license under French law and +abiding by the rules of distribution of free software. You can use, +modify and/ or redistribute the software under the terms of the CeCILL +license as circulated by CEA, CNRS and INRIA at the following URL +"http://www.cecill.info". + +As a counterpart to the access to the source code and rights to copy, +modify and redistribute granted by the license, users are provided only +with a limited warranty and the software's author, the holder of the +economic rights, and the successive licensors have only limited +liability. + +In this respect, the user's attention is drawn to the risks associated +with loading, using, modifying and/or developing or reproducing the +software by the user in light of its specific status of free software, +that may mean that it is complicated to manipulate, and that also +therefore means that it is reserved for developers and experienced +professionals having in-depth computer knowledge. Users are therefore +encouraged to load and test the software's suitability as regards their +requirements in conditions enabling the security of their systems and/or +data to be ensured and, more generally, to use and operate it in the +same conditions as regards security. + +The fact that you are presently reading this means that you have had +knowledge of the CeCILL license and that you accept its terms. ++*/ + #ifndef __COSMO_QUEUE_HPP #define __COSMO_QUEUE_HPP diff --git a/src/bsp_simple.hpp b/src/bsp_simple.hpp index b61bcf3..298256e 100644 --- a/src/bsp_simple.hpp +++ b/src/bsp_simple.hpp @@ -1,3 +1,38 @@ +/*+ +This is CosmoTool (./src/bsp_simple.hpp) -- Copyright (C) Guilhem Lavaux (2007-2013) + +guilhem.lavaux@gmail.com + +This software is a computer program whose purpose is to provide a toolbox for cosmological +data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...) + +This software is governed by the CeCILL license under French law and +abiding by the rules of distribution of free software. You can use, +modify and/ or redistribute the software under the terms of the CeCILL +license as circulated by CEA, CNRS and INRIA at the following URL +"http://www.cecill.info". + +As a counterpart to the access to the source code and rights to copy, +modify and redistribute granted by the license, users are provided only +with a limited warranty and the software's author, the holder of the +economic rights, and the successive licensors have only limited +liability. + +In this respect, the user's attention is drawn to the risks associated +with loading, using, modifying and/or developing or reproducing the +software by the user in light of its specific status of free software, +that may mean that it is complicated to manipulate, and that also +therefore means that it is reserved for developers and experienced +professionals having in-depth computer knowledge. Users are therefore +encouraged to load and test the software's suitability as regards their +requirements in conditions enabling the security of their systems and/or +data to be ensured and, more generally, to use and operate it in the +same conditions as regards security. + +The fact that you are presently reading this means that you have had +knowledge of the CeCILL license and that you accept its terms. ++*/ + #ifndef __COSMOTOOL_SIMPLE_BSP_HPP #define __COSMOTOOL_SIMPLE_BSP_HPP diff --git a/src/cic.cpp b/src/cic.cpp index 71a106d..8a24597 100644 --- a/src/cic.cpp +++ b/src/cic.cpp @@ -1,3 +1,38 @@ +/*+ +This is CosmoTool (./src/cic.cpp) -- Copyright (C) Guilhem Lavaux (2007-2013) + +guilhem.lavaux@gmail.com + +This software is a computer program whose purpose is to provide a toolbox for cosmological +data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...) + +This software is governed by the CeCILL license under French law and +abiding by the rules of distribution of free software. You can use, +modify and/ or redistribute the software under the terms of the CeCILL +license as circulated by CEA, CNRS and INRIA at the following URL +"http://www.cecill.info". + +As a counterpart to the access to the source code and rights to copy, +modify and redistribute granted by the license, users are provided only +with a limited warranty and the software's author, the holder of the +economic rights, and the successive licensors have only limited +liability. + +In this respect, the user's attention is drawn to the risks associated +with loading, using, modifying and/or developing or reproducing the +software by the user in light of its specific status of free software, +that may mean that it is complicated to manipulate, and that also +therefore means that it is reserved for developers and experienced +professionals having in-depth computer knowledge. Users are therefore +encouraged to load and test the software's suitability as regards their +requirements in conditions enabling the security of their systems and/or +data to be ensured and, more generally, to use and operate it in the +same conditions as regards security. + +The fact that you are presently reading this means that you have had +knowledge of the CeCILL license and that you accept its terms. ++*/ + #include #include #include diff --git a/src/cic.hpp b/src/cic.hpp index d522d00..31ba069 100644 --- a/src/cic.hpp +++ b/src/cic.hpp @@ -1,3 +1,38 @@ +/*+ +This is CosmoTool (./src/cic.hpp) -- Copyright (C) Guilhem Lavaux (2007-2013) + +guilhem.lavaux@gmail.com + +This software is a computer program whose purpose is to provide a toolbox for cosmological +data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...) + +This software is governed by the CeCILL license under French law and +abiding by the rules of distribution of free software. You can use, +modify and/ or redistribute the software under the terms of the CeCILL +license as circulated by CEA, CNRS and INRIA at the following URL +"http://www.cecill.info". + +As a counterpart to the access to the source code and rights to copy, +modify and redistribute granted by the license, users are provided only +with a limited warranty and the software's author, the holder of the +economic rights, and the successive licensors have only limited +liability. + +In this respect, the user's attention is drawn to the risks associated +with loading, using, modifying and/or developing or reproducing the +software by the user in light of its specific status of free software, +that may mean that it is complicated to manipulate, and that also +therefore means that it is reserved for developers and experienced +professionals having in-depth computer knowledge. Users are therefore +encouraged to load and test the software's suitability as regards their +requirements in conditions enabling the security of their systems and/or +data to be ensured and, more generally, to use and operate it in the +same conditions as regards security. + +The fact that you are presently reading this means that you have had +knowledge of the CeCILL license and that you accept its terms. ++*/ + #ifndef __CICFILTER_HPP #define __CICFILTER_HPP diff --git a/src/config.hpp b/src/config.hpp index b046e2f..8fd7b2d 100644 --- a/src/config.hpp +++ b/src/config.hpp @@ -1,3 +1,38 @@ +/*+ +This is CosmoTool (./src/config.hpp) -- Copyright (C) Guilhem Lavaux (2007-2013) + +guilhem.lavaux@gmail.com + +This software is a computer program whose purpose is to provide a toolbox for cosmological +data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...) + +This software is governed by the CeCILL license under French law and +abiding by the rules of distribution of free software. You can use, +modify and/ or redistribute the software under the terms of the CeCILL +license as circulated by CEA, CNRS and INRIA at the following URL +"http://www.cecill.info". + +As a counterpart to the access to the source code and rights to copy, +modify and redistribute granted by the license, users are provided only +with a limited warranty and the software's author, the holder of the +economic rights, and the successive licensors have only limited +liability. + +In this respect, the user's attention is drawn to the risks associated +with loading, using, modifying and/or developing or reproducing the +software by the user in light of its specific status of free software, +that may mean that it is complicated to manipulate, and that also +therefore means that it is reserved for developers and experienced +professionals having in-depth computer knowledge. Users are therefore +encouraged to load and test the software's suitability as regards their +requirements in conditions enabling the security of their systems and/or +data to be ensured and, more generally, to use and operate it in the +same conditions as regards security. + +The fact that you are presently reading this means that you have had +knowledge of the CeCILL license and that you accept its terms. ++*/ + #ifndef __COSMOTOOL_CONFIG_HPP #define __COSMOTOOL_CONFIG_HPP diff --git a/src/cosmopower.cpp b/src/cosmopower.cpp new file mode 100644 index 0000000..a96a6bb --- /dev/null +++ b/src/cosmopower.cpp @@ -0,0 +1,341 @@ +/*+ +This is CosmoTool (./src/cosmopower.cpp) -- Copyright (C) Guilhem Lavaux (2007-2013) + +guilhem.lavaux@gmail.com + +This software is a computer program whose purpose is to provide a toolbox for cosmological +data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...) + +This software is governed by the CeCILL license under French law and +abiding by the rules of distribution of free software. You can use, +modify and/ or redistribute the software under the terms of the CeCILL +license as circulated by CEA, CNRS and INRIA at the following URL +"http://www.cecill.info". + +As a counterpart to the access to the source code and rights to copy, +modify and redistribute granted by the license, users are provided only +with a limited warranty and the software's author, the holder of the +economic rights, and the successive licensors have only limited +liability. + +In this respect, the user's attention is drawn to the risks associated +with loading, using, modifying and/or developing or reproducing the +software by the user in light of its specific status of free software, +that may mean that it is complicated to manipulate, and that also +therefore means that it is reserved for developers and experienced +professionals having in-depth computer knowledge. Users are therefore +encouraged to load and test the software's suitability as regards their +requirements in conditions enabling the security of their systems and/or +data to be ensured and, more generally, to use and operate it in the +same conditions as regards security. + +The fact that you are presently reading this means that you have had +knowledge of the CeCILL license and that you accept its terms. ++*/ + +#include +#include +#include +#include +#include +#include +#include +#include "cosmopower.hpp" + +using namespace std; +using namespace CosmoTool; + +#define USE_GSL +#define TOLERANCE 1e-6 +#define NUM_ITERATION 8000 + + +CosmoPower::CosmoPower() +{ + eval = &CosmoPower::powerEfstathiou; + + n = 1.0; + K0 = 1; + V_LG_CMB = 627; + + CMB_VECTOR[0] = 56.759; + CMB_VECTOR[1] = -540.02; + CMB_VECTOR[2] = 313.50; + + h = 0.719; + SIGMA8 = 0.77; + OMEGA_B = 0.043969; + OMEGA_C = 0.21259; + + Theta_27 = 2.728/2.7; + + updateCosmology(); +} + +/* + * This is \hat{tophat} + */ +static double tophatFilter(double u) +{ + if (u != 0) + return 3 / (u*u*u) * (sin(u) - u * cos(u)); + else + return 1; +} + +static double powC(double q, double alpha_c) +{ + return 14.2 / alpha_c + 386 / (1 + 69.9 * pow(q, 1.08)); +} + +static double T_tilde_0(double q, double alpha_c, double beta_c) +{ + double a = log(M_E + 1.8 * beta_c * q); + return a / ( a + powC(q, alpha_c) * q * q); +} + +static double j_0(double x) +{ + if (x == 0) + return 1.0; + return sin(x)/x; +} + +static double powG(double y) +{ + double a = sqrt(1 + y); + + return y * (-6 * a + (2 + 3 * y) *log((a + 1)/(a - 1))); +} + +double CosmoPower::powerEfstathiou(double k) +{ + double a = 6.4/Gamma0; + double b = 3/Gamma0; + double c = 1.7/Gamma0; + double nu = 1.13; + + double f = (a*k) + pow(b*k,1.5) + pow(c*k,2); + + // EFSTATHIOU ET AL. (1992) + return normPower * pow(k,n) * pow(1+pow(f,nu),(-2/nu)); +} + +double CosmoPower::powerHuWiggles(double k) +{ + // EISENSTEIN ET HU (1998) + // FULL POWER SPECTRUM WITH BARYONS AND WIGGLES + + double k_silk = 1.6 * pow(OMEGA_B * h * h, 0.52) * pow(OmegaEff, 0.73) * (1 + pow(10.4 * OmegaEff, -0.95)); + double z_eq = 2.50e4 * OmegaEff * pow(Theta_27, -4); + double s = 44.5 * log(9.83 / OmegaEff) / (sqrt(1 + 10 * pow(OMEGA_B * h * h, 0.75))); + double f = 1 / (1 + pow(k * s / 5.4, 4)); + double k_eq = 7.46e-2 * OmegaEff * pow(Theta_27, -2); + double a1 = pow(46.9 * OmegaEff, 0.670) * (1 + pow(32.1 * OmegaEff, -0.532)); + double a2 = pow(12.0 * OmegaEff, 0.424) * (1 + pow(45.0 * OmegaEff, -0.582)); + double alpha_c = pow(a1, -OMEGA_B/ OMEGA_0) * pow(a2, -pow(OMEGA_B / OMEGA_0, 3)); + + double q = k / (13.41 * k_eq); + double b1_betac = 0.944 * 1/(1 + pow(458 * OmegaEff, -0.708)); + double b2_betac = pow(0.395 * OmegaEff, -0.0266); + double beta_c = 1/ ( 1 + b1_betac * (pow(OMEGA_C / OMEGA_0, b2_betac) - 1) ); + double T_c = f * T_tilde_0(q, 1, beta_c) + (1 - f) * T_tilde_0(q, alpha_c, beta_c); + + double b1_zd = 0.313 * pow(OmegaEff, -0.419) * (1 + 0.607 * pow(OmegaEff, 0.674)); + double b2_zd = 0.238 * pow(OmegaEff, 0.223); + double z_d = 1291 * pow(OmegaEff, 0.251) / (1 + 0.659 * pow(OmegaEff, 0.828)) * (1 + b1_zd * pow(OmegaEff, b2_zd)); + double R_d = 31.5 * OMEGA_B * h * h * pow(Theta_27, -4) * 1e3 / z_d; + + double alpha_b = 2.07 * k_eq * s * pow(1 + R_d, -0.75) * powG((1 + z_eq)/(1 + z_d)); + double beta_b = 0.5 + OMEGA_B / OMEGA_0 + (3 - 2 * OMEGA_B / OMEGA_0) * sqrt(pow(17.2 * OmegaEff, 2) + 1); + double beta_node = 8.41 * pow(OmegaEff, 0.435); + double s_tilde = s * pow(1 + pow(beta_node / (k * s), 3), -1./3); + + double T_b = (T_tilde_0(q, 1, 1) / (1 + pow(k * s / 5.2, 2)) + alpha_b / (1 + pow(beta_b / (k * s), 3)) * exp(-pow(k/k_silk, 1.4))) * j_0(k * s_tilde); + + double T_k = OMEGA_B/OMEGA_0 * T_b + OMEGA_C/OMEGA_0 * T_c; + + return normPower * pow(k,n) * T_k * T_k; +} + +double CosmoPower::powerHuBaryons(double k) +{ + double s = 44.5 * log(9.83 / OmegaEff) / (sqrt(1 + 10 * pow(OMEGA_B * h * h, 0.75))); + double alpha_Gamma = 1 - 0.328 * log(431 * OmegaEff) * OMEGA_B / OMEGA_0 + 0.38 * log(22.3 * OmegaEff) * pow(OMEGA_B / OMEGA_0, 2); + double GammaEff = OMEGA_0 * h * (alpha_Gamma + (1 - alpha_Gamma)/(1 + pow(0.43 * k * s, 4))); + double q = k/(h*GammaEff) * pow(Theta_27, 2); + double L_0 = log(2 * M_E + 1.8 * q); + double C_0 = 14.2 + 731 / (1 + 62.5 * q); + double T0 = L_0 / (L_0 + C_0 * q * q); + + return normPower * pow(k,n) * T0 * T0; +} + +double CosmoPower::powerOld(double k) +{ + static const double l = 1/(Omega * h*h); + static const double alpha = 1.7 * l, beta = 9.0 * pow(l, 1.5), gamma = l*l; + return normPower * pow(k,n) * pow(1 + alpha * k + beta * pow(k,1.5) + gamma *k*k,-2); +} + +double CosmoPower::powerSugiyama(double k) +{ + double q = k * Theta_27*Theta_27 / (OmegaEff * exp(-OMEGA_B - sqrt(h/0.5)*OMEGA_B/OMEGA_0)); + double L0 = log(2*M_E + 1.8 * q); + double C0 = 14.2 + 731 / (1 + 62.5 * q); + double T_k = L0 / (L0 + C0 * q*q); + + return normPower * pow(k,n) * T_k * T_k; +} + +double CosmoPower::powerBardeen(double k) +{ + double q = k / (OmegaEff); + + double poly = 1 + 3.89 * q + pow(16.1*q,2) + pow(5.46*q,3) + pow(6.71*q,4); + double T_k = log(1+2.34*q)/(2.34*q) * pow(poly,-0.25); + + return normPower * pow(k,n) * T_k * T_k; +} + +double CosmoPower::powerBDM(double k) +{ + k /= h*h; + double alpha1 = 190; + double Gmu = 4.6; + double alpha2 = 11.5; + double alpha3 = 11; + double alpha4 = 12.55; + double alpha5 = 0.0004; + return normPower*k*alpha1*alpha1*Gmu*Gmu/(1+(alpha2*k)+pow(alpha3*k,2)+pow(alpha4*k,3))*pow(1+pow(alpha5/k,2), -2); +} + +double CosmoPower::powerTest(double k) +{ + return 1/(1+k*k); +} + +/* + * This function computes the normalization of the power spectrum. It requests + * a sigma8 (density fluctuations within 8 Mpc/h) + */ +static double gslPowSpecNorm(double k, void *params) +{ + CosmoPower *c = (CosmoPower *)params; + + return c->integrandNormalize(k); +} + +double CosmoPower::integrandNormalize(double x) +{ + double k = (1-x)/x; + double f = tophatFilter(k*8.0/h); + return power(k)*k*k*f*f/(x*x); +} + +void CosmoPower::normalize() +{ + double normVal = 0; + double abserr; + gsl_integration_workspace *w = gsl_integration_workspace_alloc(NUM_ITERATION); + gsl_function f; + + f.function = gslPowSpecNorm; + f.params = this; + + normPower = 1; + + ofstream ff("PP_k.txt"); + for (int i = 0; i < 100; i++) + { + double k = pow(10.0, 4.0*i/100.-2); + ff << k << " " << power(k) << endl; + } + + // gsl_integration_qagiu(&f, 0, 0, TOLERANCE, NUM_ITERATION, w, &normVal, &abserr); + gsl_integration_qag(&f, 0, 1, 0, TOLERANCE, NUM_ITERATION, GSL_INTEG_GAUSS61, w, &normVal, &abserr); + gsl_integration_workspace_free(w); + + normVal /= (2*M_PI*M_PI); + + normPower = SIGMA8*SIGMA8/normVal; +} + +void CosmoPower::updateCosmology() +{ + OMEGA_0 = OMEGA_B+OMEGA_C; + Omega = OMEGA_0; + beta = pow(OMEGA_0, 5./9); + OmegaEff = OMEGA_0 * h * h; + Gamma0 = OMEGA_0 * h * h; + omega_B = OMEGA_B * h * h; + omega_C = OMEGA_C * h * h; +} + +void CosmoPower::updatePhysicalCosmology() +{ + OMEGA_B = omega_B / (h*h); + OMEGA_C = omega_C / (h*h); + OMEGA_0 = Gamma0 / (h*h); + beta = pow(OMEGA_0, 5./9); +} + +double CosmoPower::eval_theta_theta(double k) +{ + // Jennings (2012) fit + double P_deltadelta = power(k); + + static const double alpha0 = -12480.5, alpha1 = 1.824, alpha2 = 2165.87, alpha3=1.796; + if (k > 0.3) + return 0; + double r =(alpha0*sqrt(P_deltadelta) + alpha1*P_deltadelta*P_deltadelta)/(alpha2 + alpha3*P_deltadelta); + assert(P_deltadelta > 0); + + if (r < 0) + return 0; + return r; +} + +double CosmoPower::power(double k) +{ + return (this->*eval)(k); +} + + +void CosmoPower::setFunction(CosmoFunction f) +{ + switch (f) + { + case POWER_EFSTATHIOU: + eval = &CosmoPower::powerEfstathiou; + break; + case HU_WIGGLES: + eval = &CosmoPower::powerHuWiggles; + break; + case HU_BARYON: + eval = &CosmoPower::powerHuBaryons; + break; + case OLD_POWERSPECTRUM: + eval = &CosmoPower::powerOld; + break; + case POWER_BARDEEN: + eval = &CosmoPower::powerBardeen; + break; + case POWER_SUGIYAMA: + eval = &CosmoPower::powerSugiyama; + break; + case POWER_BDM: + eval = &CosmoPower::powerBDM; + break; + case POWER_TEST: + eval = &CosmoPower::powerTest; + break; + default: + abort(); + } +} + +void CosmoPower::setNormalization(double A_K) +{ + normPower = A_K/power(0.002); +} diff --git a/src/cosmopower.hpp b/src/cosmopower.hpp new file mode 100644 index 0000000..dd86cf1 --- /dev/null +++ b/src/cosmopower.hpp @@ -0,0 +1,109 @@ +/*+ +This is CosmoTool (./src/cosmopower.hpp) -- Copyright (C) Guilhem Lavaux (2007-2013) + +guilhem.lavaux@gmail.com + +This software is a computer program whose purpose is to provide a toolbox for cosmological +data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...) + +This software is governed by the CeCILL license under French law and +abiding by the rules of distribution of free software. You can use, +modify and/ or redistribute the software under the terms of the CeCILL +license as circulated by CEA, CNRS and INRIA at the following URL +"http://www.cecill.info". + +As a counterpart to the access to the source code and rights to copy, +modify and redistribute granted by the license, users are provided only +with a limited warranty and the software's author, the holder of the +economic rights, and the successive licensors have only limited +liability. + +In this respect, the user's attention is drawn to the risks associated +with loading, using, modifying and/or developing or reproducing the +software by the user in light of its specific status of free software, +that may mean that it is complicated to manipulate, and that also +therefore means that it is reserved for developers and experienced +professionals having in-depth computer knowledge. Users are therefore +encouraged to load and test the software's suitability as regards their +requirements in conditions enabling the security of their systems and/or +data to be ensured and, more generally, to use and operate it in the +same conditions as regards security. + +The fact that you are presently reading this means that you have had +knowledge of the CeCILL license and that you accept its terms. ++*/ + +#ifndef _COSMOPOWER_HPP +#define _COSMOPOWER_HPP + +namespace CosmoTool { + + class CosmoPower + { + public: + // PRIMARY VARIABLES + double n; + double K0; + double V_LG_CMB; + + double CMB_VECTOR[3]; + + double h; + double SIGMA8; + double OMEGA_B; + double OMEGA_C; + double omega_B; + double omega_C; + double Theta_27; + + // DERIVED VARIABLES + double OMEGA_0; + double Omega; + double beta; + double OmegaEff; + double Gamma0; + double normPower; + + enum CosmoFunction + { + POWER_EFSTATHIOU, + HU_WIGGLES, + HU_BARYON, + OLD_POWERSPECTRUM, + POWER_BARDEEN, + POWER_SUGIYAMA, + POWER_BDM, + POWER_TEST + }; + + CosmoPower(); + + void setFunction(CosmoFunction f); + + void updateCosmology(); + void updatePhysicalCosmology(); + void normalize(); + void setNormalization(double A_K); + + + double eval_theta_theta(double k); + double power(double k); + + double integrandNormalize(double k); + private: + double (CosmoPower::*eval)(double); + + double powerEfstathiou(double k); + double powerHuWiggles(double k); + double powerHuBaryons(double k); + double powerOld(double k); + double powerBardeen(double k); + double powerSugiyama(double k); + double powerBDM(double k); + double powerTest(double k); + + }; + +}; + +#endif diff --git a/src/dinterpolate.hpp b/src/dinterpolate.hpp index 9ce6281..1bf3322 100644 --- a/src/dinterpolate.hpp +++ b/src/dinterpolate.hpp @@ -1,3 +1,38 @@ +/*+ +This is CosmoTool (./src/dinterpolate.hpp) -- Copyright (C) Guilhem Lavaux (2007-2013) + +guilhem.lavaux@gmail.com + +This software is a computer program whose purpose is to provide a toolbox for cosmological +data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...) + +This software is governed by the CeCILL license under French law and +abiding by the rules of distribution of free software. You can use, +modify and/ or redistribute the software under the terms of the CeCILL +license as circulated by CEA, CNRS and INRIA at the following URL +"http://www.cecill.info". + +As a counterpart to the access to the source code and rights to copy, +modify and redistribute granted by the license, users are provided only +with a limited warranty and the software's author, the holder of the +economic rights, and the successive licensors have only limited +liability. + +In this respect, the user's attention is drawn to the risks associated +with loading, using, modifying and/or developing or reproducing the +software by the user in light of its specific status of free software, +that may mean that it is complicated to manipulate, and that also +therefore means that it is reserved for developers and experienced +professionals having in-depth computer knowledge. Users are therefore +encouraged to load and test the software's suitability as regards their +requirements in conditions enabling the security of their systems and/or +data to be ensured and, more generally, to use and operate it in the +same conditions as regards security. + +The fact that you are presently reading this means that you have had +knowledge of the CeCILL license and that you accept its terms. ++*/ + #ifndef __COSMO_DINTERPOLATE_HPP #define __COSMO_DINTERPOLATE_HPP @@ -29,6 +64,7 @@ namespace CosmoTool { uint32_t numSimplex; uint32_t *simplex_list; gsl_eigen_symmv_workspace *eigen_work; + bool *disable_simplex; /** * This construct the interpolator. The construction is time consuming so @@ -53,6 +89,7 @@ namespace CosmoTool { this->simplex_list = simplex_list; this->numPoints = numPoints; this->numSimplex = numSimplex; + this->disable_simplex = new bool[numSimplex]; buildPreweight(); buildQuickAccess(); @@ -67,6 +104,7 @@ namespace CosmoTool { delete quickAccess; delete[] point_to_simplex_list_base; delete[] all_preweight; + delete[] disable_simplex; gsl_eigen_symmv_free(eigen_work); } diff --git a/src/dinterpolate.tcc b/src/dinterpolate.tcc index e18a7e6..e0cf763 100644 --- a/src/dinterpolate.tcc +++ b/src/dinterpolate.tcc @@ -1,3 +1,6 @@ +#include +#include +#include #include #include #include @@ -21,7 +24,8 @@ namespace CosmoTool { for (uint32_t i = 0; i < (N+1)*numSimplex; i++) { assert(simplex_list[i] < numPoints); - numSimplex_by_point[simplex_list[i]]++; + if (!disable_simplex[i/(N+1)]) + numSimplex_by_point[simplex_list[i]]++; } // Compute the total number and the index for accessing lists. @@ -37,7 +41,11 @@ namespace CosmoTool { { for (int j = 0; j <= N; j++) { - uint32_t p = simplex_list[(N+1)*i+j]; + uint32_t s = (N+1)*i+j; + if (disable_simplex[i]) + continue; + + uint32_t p = simplex_list[s]; assert(index_by_point[p] < point_to_simplex_size); point_to_simplex_list_base[index_by_point[p]] = i; ++index_by_point[p]; @@ -48,7 +56,9 @@ namespace CosmoTool { for (uint32_t i = 0; i < numPoints; i++) { // check assertion - assert((i==0 && index_by_point[0]==numSimplex_by_point[0]) || ((index_by_point[i]-index_by_point[i-1]) == (numSimplex_by_point[i]+1))); + assert((i==0 && index_by_point[0]==numSimplex_by_point[0]) + || + ((index_by_point[i]-index_by_point[i-1]) == (numSimplex_by_point[i]+1))); assert(index_by_point[i] < point_to_simplex_size); point_to_simplex_list_base[index_by_point[i]] = -1; } @@ -80,6 +90,7 @@ namespace CosmoTool { double preweight[N*N]; double preweight_inverse[N*N]; gsl_permutation *p = gsl_permutation_alloc(N); + uint32_t numDisabled = 0; all_preweight = new PType[N*N*numSimplex]; @@ -105,13 +116,38 @@ namespace CosmoTool { gsl_linalg_LU_decomp(&M.matrix, p, &signum); double a = fabs(gsl_linalg_LU_det(&M.matrix, signum)); if (a < 1e-10) - throw InvalidArgumentException("Invalid tesselation. One simplex is coplanar."); - gsl_linalg_LU_invert(&M.matrix, p, &iM.matrix); + { +#ifdef DEBUG + for (int j = 0; j < N; j++) + { + PType xref = positions[pref][j]; + for (int k = 0; k < N; k++) + { + preweight[j*N + k] = positions[simplex_list[k+base+1]][j] - xref; + } + } + std::ofstream f("matrix.txt"); + for (int j = 0; j < N*N; j++) + f << std::setprecision(12) << preweight[j] << std::endl; + throw InvalidArgumentException("Invalid tesselation. One simplex is coplanar."); +#else + gsl_matrix_set_zero(&iM.matrix); + disable_simplex[i] = true; + numDisabled++; +#endif + } + else { + gsl_linalg_LU_invert(&M.matrix, p, &iM.matrix); + disable_simplex[i] = false; + } + for (int j = 0; j < N*N; j++) all_preweight[N*N*i + j] = preweight_inverse[j]; } + std::cout << "Number of disabled simplices: " << numDisabled << std::endl; + gsl_permutation_free(p); } @@ -166,6 +202,9 @@ namespace CosmoTool { template bool DelaunayInterpolate::checkPointInSimplex(const CoordType& pos, uint32_t simplex) { + if (disable_simplex[simplex]) + return false; + uint32_t *desc_simplex = &simplex_list[simplex*(N+1)]; CoordType *p[N+1], v[N], hyper; diff --git a/src/eskow.hpp b/src/eskow.hpp index ac9b292..987cb88 100644 --- a/src/eskow.hpp +++ b/src/eskow.hpp @@ -1,3 +1,38 @@ +/*+ +This is CosmoTool (./src/eskow.hpp) -- Copyright (C) Guilhem Lavaux (2007-2013) + +guilhem.lavaux@gmail.com + +This software is a computer program whose purpose is to provide a toolbox for cosmological +data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...) + +This software is governed by the CeCILL license under French law and +abiding by the rules of distribution of free software. You can use, +modify and/ or redistribute the software under the terms of the CeCILL +license as circulated by CEA, CNRS and INRIA at the following URL +"http://www.cecill.info". + +As a counterpart to the access to the source code and rights to copy, +modify and redistribute granted by the license, users are provided only +with a limited warranty and the software's author, the holder of the +economic rights, and the successive licensors have only limited +liability. + +In this respect, the user's attention is drawn to the risks associated +with loading, using, modifying and/or developing or reproducing the +software by the user in light of its specific status of free software, +that may mean that it is complicated to manipulate, and that also +therefore means that it is reserved for developers and experienced +professionals having in-depth computer knowledge. Users are therefore +encouraged to load and test the software's suitability as regards their +requirements in conditions enabling the security of their systems and/or +data to be ensured and, more generally, to use and operate it in the +same conditions as regards security. + +The fact that you are presently reading this means that you have had +knowledge of the CeCILL license and that you accept its terms. ++*/ + #ifndef __ESKOW_CHOLESKY_HPP #define __ESKOW_CHOLESKY_HPP diff --git a/src/field.hpp b/src/field.hpp index 25ec447..45f6eac 100644 --- a/src/field.hpp +++ b/src/field.hpp @@ -1,3 +1,38 @@ +/*+ +This is CosmoTool (./src/field.hpp) -- Copyright (C) Guilhem Lavaux (2007-2013) + +guilhem.lavaux@gmail.com + +This software is a computer program whose purpose is to provide a toolbox for cosmological +data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...) + +This software is governed by the CeCILL license under French law and +abiding by the rules of distribution of free software. You can use, +modify and/ or redistribute the software under the terms of the CeCILL +license as circulated by CEA, CNRS and INRIA at the following URL +"http://www.cecill.info". + +As a counterpart to the access to the source code and rights to copy, +modify and redistribute granted by the license, users are provided only +with a limited warranty and the software's author, the holder of the +economic rights, and the successive licensors have only limited +liability. + +In this respect, the user's attention is drawn to the risks associated +with loading, using, modifying and/or developing or reproducing the +software by the user in light of its specific status of free software, +that may mean that it is complicated to manipulate, and that also +therefore means that it is reserved for developers and experienced +professionals having in-depth computer knowledge. Users are therefore +encouraged to load and test the software's suitability as regards their +requirements in conditions enabling the security of their systems and/or +data to be ensured and, more generally, to use and operate it in the +same conditions as regards security. + +The fact that you are presently reading this means that you have had +knowledge of the CeCILL license and that you accept its terms. ++*/ + #ifndef __COSMOTOOL_FIELD #define __COSMOTOOL_FIELD diff --git a/src/fixArray.hpp b/src/fixArray.hpp index d1516e4..0d2e26d 100644 --- a/src/fixArray.hpp +++ b/src/fixArray.hpp @@ -1,3 +1,38 @@ +/*+ +This is CosmoTool (./src/fixArray.hpp) -- Copyright (C) Guilhem Lavaux (2007-2013) + +guilhem.lavaux@gmail.com + +This software is a computer program whose purpose is to provide a toolbox for cosmological +data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...) + +This software is governed by the CeCILL license under French law and +abiding by the rules of distribution of free software. You can use, +modify and/ or redistribute the software under the terms of the CeCILL +license as circulated by CEA, CNRS and INRIA at the following URL +"http://www.cecill.info". + +As a counterpart to the access to the source code and rights to copy, +modify and redistribute granted by the license, users are provided only +with a limited warranty and the software's author, the holder of the +economic rights, and the successive licensors have only limited +liability. + +In this respect, the user's attention is drawn to the risks associated +with loading, using, modifying and/or developing or reproducing the +software by the user in light of its specific status of free software, +that may mean that it is complicated to manipulate, and that also +therefore means that it is reserved for developers and experienced +professionals having in-depth computer knowledge. Users are therefore +encouraged to load and test the software's suitability as regards their +requirements in conditions enabling the security of their systems and/or +data to be ensured and, more generally, to use and operate it in the +same conditions as regards security. + +The fact that you are presently reading this means that you have had +knowledge of the CeCILL license and that you accept its terms. ++*/ + #ifndef __FIX_ARRAY_HPP #define __FIX_ARRAY_HPP diff --git a/src/fortran.cpp b/src/fortran.cpp index 727a6d0..79641ac 100644 --- a/src/fortran.cpp +++ b/src/fortran.cpp @@ -1,3 +1,38 @@ +/*+ +This is CosmoTool (./src/fortran.cpp) -- Copyright (C) Guilhem Lavaux (2007-2013) + +guilhem.lavaux@gmail.com + +This software is a computer program whose purpose is to provide a toolbox for cosmological +data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...) + +This software is governed by the CeCILL license under French law and +abiding by the rules of distribution of free software. You can use, +modify and/ or redistribute the software under the terms of the CeCILL +license as circulated by CEA, CNRS and INRIA at the following URL +"http://www.cecill.info". + +As a counterpart to the access to the source code and rights to copy, +modify and redistribute granted by the license, users are provided only +with a limited warranty and the software's author, the holder of the +economic rights, and the successive licensors have only limited +liability. + +In this respect, the user's attention is drawn to the risks associated +with loading, using, modifying and/or developing or reproducing the +software by the user in light of its specific status of free software, +that may mean that it is complicated to manipulate, and that also +therefore means that it is reserved for developers and experienced +professionals having in-depth computer knowledge. Users are therefore +encouraged to load and test the software's suitability as regards their +requirements in conditions enabling the security of their systems and/or +data to be ensured and, more generally, to use and operate it in the +same conditions as regards security. + +The fact that you are presently reading this means that you have had +knowledge of the CeCILL license and that you accept its terms. ++*/ + #include "config.hpp" #include #include diff --git a/src/fortran.hpp b/src/fortran.hpp index 651d0fd..d63e8d3 100644 --- a/src/fortran.hpp +++ b/src/fortran.hpp @@ -1,3 +1,38 @@ +/*+ +This is CosmoTool (./src/fortran.hpp) -- Copyright (C) Guilhem Lavaux (2007-2013) + +guilhem.lavaux@gmail.com + +This software is a computer program whose purpose is to provide a toolbox for cosmological +data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...) + +This software is governed by the CeCILL license under French law and +abiding by the rules of distribution of free software. You can use, +modify and/ or redistribute the software under the terms of the CeCILL +license as circulated by CEA, CNRS and INRIA at the following URL +"http://www.cecill.info". + +As a counterpart to the access to the source code and rights to copy, +modify and redistribute granted by the license, users are provided only +with a limited warranty and the software's author, the holder of the +economic rights, and the successive licensors have only limited +liability. + +In this respect, the user's attention is drawn to the risks associated +with loading, using, modifying and/or developing or reproducing the +software by the user in light of its specific status of free software, +that may mean that it is complicated to manipulate, and that also +therefore means that it is reserved for developers and experienced +professionals having in-depth computer knowledge. Users are therefore +encouraged to load and test the software's suitability as regards their +requirements in conditions enabling the security of their systems and/or +data to be ensured and, more generally, to use and operate it in the +same conditions as regards security. + +The fact that you are presently reading this means that you have had +knowledge of the CeCILL license and that you accept its terms. ++*/ + #ifndef __COSMO_FORTRAN_HPP #define __COSMO_FORTRAN_HPP diff --git a/src/fourier/base_types.hpp b/src/fourier/base_types.hpp index 9f54320..f39e3a5 100644 --- a/src/fourier/base_types.hpp +++ b/src/fourier/base_types.hpp @@ -1,3 +1,38 @@ +/*+ +This is CosmoTool (./src/fourier/base_types.hpp) -- Copyright (C) Guilhem Lavaux (2007-2013) + +guilhem.lavaux@gmail.com + +This software is a computer program whose purpose is to provide a toolbox for cosmological +data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...) + +This software is governed by the CeCILL license under French law and +abiding by the rules of distribution of free software. You can use, +modify and/ or redistribute the software under the terms of the CeCILL +license as circulated by CEA, CNRS and INRIA at the following URL +"http://www.cecill.info". + +As a counterpart to the access to the source code and rights to copy, +modify and redistribute granted by the license, users are provided only +with a limited warranty and the software's author, the holder of the +economic rights, and the successive licensors have only limited +liability. + +In this respect, the user's attention is drawn to the risks associated +with loading, using, modifying and/or developing or reproducing the +software by the user in light of its specific status of free software, +that may mean that it is complicated to manipulate, and that also +therefore means that it is reserved for developers and experienced +professionals having in-depth computer knowledge. Users are therefore +encouraged to load and test the software's suitability as regards their +requirements in conditions enabling the security of their systems and/or +data to be ensured and, more generally, to use and operate it in the +same conditions as regards security. + +The fact that you are presently reading this means that you have had +knowledge of the CeCILL license and that you accept its terms. ++*/ + #ifndef __BASE_FOURIER_TYPES_HPP #define __BASE_FOURIER_TYPES_HPP @@ -24,14 +59,28 @@ namespace CosmoTool typedef Eigen::Map MapType; typedef Eigen::Map ConstMapType; typedef FourierMap > FourierMapType; + typedef boost::shared_ptr FourierMapPtr; + typedef boost::shared_ptr > SpectrumFunctionPtr; - virtual boost::shared_ptr - newRandomFourier(gsl_rng *rng, const FourierMapType& like_map) const = 0; + virtual ~SpectrumFunction() {} - virtual void mul(FourierMap >& m) const = 0; + virtual void + newRandomFourier(gsl_rng *rng, FourierMapType& out_map) const = 0; + + virtual SpectrumFunctionPtr copy() const = 0; + + virtual const T *data() const { return 0; } + virtual const int *dof() const { return 0; } + virtual long size() const { return -1; } + + virtual void sqrt() = 0; + + virtual void mul(FourierMapType& m) const = 0; + virtual void mul_sqrt(FourierMapType& m) const = 0; + virtual void mul_inv(FourierMapType& m) const = 0; + virtual void mul_inv_sqrt(FourierMapType& m) const = 0; }; - template class FourierMap { @@ -60,6 +109,12 @@ namespace CosmoTool return ConstMapType(data(), size()); } + FourierMap& operator=(const FourierMap& m) + { + eigen() = m.eigen(); + return *this; + } + void sqrt() { MapType m = eigen(); @@ -76,7 +131,7 @@ namespace CosmoTool { assert(size() == map2->size()); MapType m(data(), size()); - MapType m2(map2->data(), map2->size()); + ConstMapType m2(map2->data(), map2->size()); m *= m2; } @@ -130,6 +185,18 @@ namespace CosmoTool virtual void synthesis_conjugate() = 0; }; + template + class MapUtilityFunction + { + public: + typedef SpectrumFunction Spectrum; + typedef boost::shared_ptr Spectrum_ptr; + typedef FourierMap > FMap; + + virtual Spectrum_ptr estimateSpectrumFromMap(const FMap& m) const = 0; + virtual Spectrum_ptr newSpectrumFromRaw(T *data, long size, + const Spectrum& like_spec) const = 0; + }; }; diff --git a/src/fourier/details/euclidian_maps.hpp b/src/fourier/details/euclidian_maps.hpp new file mode 100644 index 0000000..53b0278 --- /dev/null +++ b/src/fourier/details/euclidian_maps.hpp @@ -0,0 +1,267 @@ +/*+ +This is CosmoTool (./src/fourier/details/euclidian_maps.hpp) -- Copyright (C) Guilhem Lavaux (2007-2013) + +guilhem.lavaux@gmail.com + +This software is a computer program whose purpose is to provide a toolbox for cosmological +data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...) + +This software is governed by the CeCILL license under French law and +abiding by the rules of distribution of free software. You can use, +modify and/ or redistribute the software under the terms of the CeCILL +license as circulated by CEA, CNRS and INRIA at the following URL +"http://www.cecill.info". + +As a counterpart to the access to the source code and rights to copy, +modify and redistribute granted by the license, users are provided only +with a limited warranty and the software's author, the holder of the +economic rights, and the successive licensors have only limited +liability. + +In this respect, the user's attention is drawn to the risks associated +with loading, using, modifying and/or developing or reproducing the +software by the user in light of its specific status of free software, +that may mean that it is complicated to manipulate, and that also +therefore means that it is reserved for developers and experienced +professionals having in-depth computer knowledge. Users are therefore +encouraged to load and test the software's suitability as regards their +requirements in conditions enabling the security of their systems and/or +data to be ensured and, more generally, to use and operate it in the +same conditions as regards security. + +The fact that you are presently reading this means that you have had +knowledge of the CeCILL license and that you accept its terms. ++*/ + +#ifndef __DETAILS_EUCLIDIAN_MAPS +#define __DETAILS_EUCLIDIAN_MAPS + + +namespace CosmoTool +{ + + template + class EuclidianFourierMapBase: public FourierMap + { + public: + typedef std::vector DimArray; + private: + boost::shared_ptr m_data; + DimArray m_dims; + long m_size; + public: + + EuclidianFourierMapBase(boost::shared_ptr indata, const DimArray& indims) + { + m_data = indata; + m_dims = indims; + m_size = 1; + for (int i = 0; i < m_dims.size(); i++) + m_size *= m_dims[i]; + } + + virtual ~EuclidianFourierMapBase() + { + } + + const DimArray& getDims() const { return m_dims; } + + virtual const T *data() const { return m_data.get(); } + virtual T *data() { return m_data.get(); } + virtual long size() const { return m_size; } + + virtual FourierMap *copy() const + { + FourierMap *m = this->mimick(); + m->eigen() = this->eigen(); + return m; + } + + }; + + template + class EuclidianFourierMapReal: public EuclidianFourierMapBase + { + public: + typedef typename EuclidianFourierMapBase::DimArray DimArray; + + EuclidianFourierMapReal(boost::shared_ptr indata, const DimArray& indims) + : EuclidianFourierMapBase(indata, indims) + {} + + virtual FourierMap *mimick() const + { + return new EuclidianFourierMapReal( + boost::shared_ptr((T *)fftw_malloc(sizeof(T)*this->size()), + std::ptr_fun(fftw_free)), + this->getDims()); + } + + virtual T dot_product(const FourierMap& other) const + throw(std::bad_cast) + { + const EuclidianFourierMapReal& m2 = dynamic_cast&>(other); + if (this->size() != m2.size()) + throw std::bad_cast(); + + return (this->eigen()*m2.eigen()).sum(); + } + }; + + template + class EuclidianFourierMapComplex: public EuclidianFourierMapBase > + { + protected: + typedef boost::shared_ptr > ptr_t; + std::vector delta_k; + int m_dim0; + bool even0, alleven; + long plane_size; + public: + typedef typename EuclidianFourierMapBase >::DimArray DimArray; + + EuclidianFourierMapComplex(ptr_t indata, + int dim0, + const DimArray& indims, + const std::vector& dk) + : EuclidianFourierMapBase >(indata, indims), delta_k(dk), m_dim0(dim0), even0((dim0 % 2)==0) + { + assert(dk.size() == indims.size()); + plane_size = 1; + alleven = true; + for (int q = 1; q < indims.size(); q++) + { + plane_size *= indims[q]; + alleven = alleven && ((indims[q]%2)==0); + } + } + + virtual FourierMap > *mimick() const + { + return + new EuclidianFourierMapComplex( + ptr_t((std::complex *) + fftw_malloc(sizeof(std::complex)*this->size()), + std::ptr_fun(fftw_free)), + m_dim0, + this->getDims(), + this->delta_k); + } + + const std::vector& get_delta_k() const + { + return this->delta_k; + } + + template + void get_Kvec(const Array& ik, Array2& kvec) const + { + const DimArray& dims = this->getDims(); + assert(ik.size() == dims.size()); + assert(kvec.size() == dims.size()); + + kvec[0] = ik[0] * delta_k[0]; + for (int q = 1; q < ik.size(); q++) + { + int dk = ik[q]; + if (dk > dims[q]/2) + dk = dk - dims[q]; + + kvec[q] = dk*delta_k[q]; + } + } + + template + void get_Kvec(long p, Array2& kvec) const + { + const DimArray& dims = this->getDims(); + DimArray d(delta_k.size()); + get_IKvec(p, d); + get_Kvec(d, kvec); + } + + void get_IKvec(long p, DimArray& ikvec) const + { + const DimArray& dims = this->getDims(); + assert(dims.size()==ikvec.size()); + for (int q = 0; q < ikvec.size(); q++) + { + ikvec[q] = p%dims[q]; + p = (p-ikvec[q])/dims[q]; + } + } + + + template + double get_K(const Array& ik) const + { + const DimArray& dims = this->getDims(); + assert(ik.size() == dims.size()); + double k2 = 0; + k2 += CosmoTool::square(ik[0]*delta_k[0]); + + for (int q = 1; q < ik.size(); q++) + { + int dk = ik[q]; + + if (dk > dims[q]/2) + dk = dk - dims[q]; + + k2 += CosmoTool::square(delta_k[q]*dk); + } + return std::sqrt(k2); + } + + double get_K(long p) const + { + const DimArray& dims = this->getDims(); + DimArray d(delta_k.size()); + for (int q = 0; q < d.size(); q++) + { + d[q] = p%dims[q]; + p = (p-d[q])/dims[q]; + } + return get_K(d); + } + + bool allDimensionsEven() const { return alleven; } + bool firstDimensionEven() const { return even0; } + + virtual std::complex dot_product(const FourierMap >& other) const + throw(std::bad_cast) + { + const EuclidianFourierMapComplex& m2 = dynamic_cast&>(other); + if (this->size() != m2.size()) + throw std::bad_cast(); + + const std::complex *d1 = this->data(); + const std::complex *d2 = m2.data(); + const DimArray& dims = this->getDims(); + int N0 = dims[0] + (even0 ? 0 : 1); + std::complex result = 0; + + for (long q0 = 1; q0 < N0-1; q0++) + { + for (long p = 0; p < plane_size; p++) + { + long idx = q0+dims[0]*p; + assert(idx < this->size()); + result += 2*(conj(d1[idx]) * d2[idx]).real(); + } + } + if (even0) + { + for (long p = 0; p < plane_size; p++) + { + long q0 = N0*p, q1 = (p+1)*N0-1; + result += conj(d1[q0]) * d2[q0]; + result += conj(d1[q1]) * d2[q1]; + } + } + return result; + } + + }; +}; + +#endif diff --git a/src/fourier/details/euclidian_spectrum_1d.hpp b/src/fourier/details/euclidian_spectrum_1d.hpp new file mode 100644 index 0000000..36b524a --- /dev/null +++ b/src/fourier/details/euclidian_spectrum_1d.hpp @@ -0,0 +1,215 @@ +/*+ +This is CosmoTool (./src/fourier/details/euclidian_spectrum_1d.hpp) -- Copyright (C) Guilhem Lavaux (2007-2013) + +guilhem.lavaux@gmail.com + +This software is a computer program whose purpose is to provide a toolbox for cosmological +data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...) + +This software is governed by the CeCILL license under French law and +abiding by the rules of distribution of free software. You can use, +modify and/ or redistribute the software under the terms of the CeCILL +license as circulated by CEA, CNRS and INRIA at the following URL +"http://www.cecill.info". + +As a counterpart to the access to the source code and rights to copy, +modify and redistribute granted by the license, users are provided only +with a limited warranty and the software's author, the holder of the +economic rights, and the successive licensors have only limited +liability. + +In this respect, the user's attention is drawn to the risks associated +with loading, using, modifying and/or developing or reproducing the +software by the user in light of its specific status of free software, +that may mean that it is complicated to manipulate, and that also +therefore means that it is reserved for developers and experienced +professionals having in-depth computer knowledge. Users are therefore +encouraged to load and test the software's suitability as regards their +requirements in conditions enabling the security of their systems and/or +data to be ensured and, more generally, to use and operate it in the +same conditions as regards security. + +The fact that you are presently reading this means that you have had +knowledge of the CeCILL license and that you accept its terms. ++*/ + +#ifndef __DETAILS_EUCLIDIAN_SPECTRUM_1D +#define __DETAILS_EUCLIDIAN_SPECTRUM_1D + +#include +#include + + +namespace CosmoTool +{ + template + class EuclidianOperator + { + public: + typedef boost::function1 Function; + + Function base, op; + T operator()(T k) { + return op(base(k)); + } + }; + + template + class EuclidianSpectrum_1D: public SpectrumFunction + { + public: + typedef boost::function1 Function; + protected: + Function f; + + static T msqrt(T a) { return std::sqrt(a); } + public: + typedef typename SpectrumFunction::FourierMapType FourierMapType; + typedef typename SpectrumFunction::SpectrumFunctionPtr SpectrumFunctionPtr; + typedef boost::shared_ptr ptr_map; + + EuclidianSpectrum_1D(Function P) + : f(P) + { + } + + void newRandomFourier(gsl_rng *rng, FourierMapType& out_map) const; + + SpectrumFunctionPtr copy() const { + return SpectrumFunctionPtr(new EuclidianSpectrum_1D(f)); + } + + void sqrt() { + EuclidianOperator o; + o.base = f; + o.op = &EuclidianSpectrum_1D::msqrt; + f = (Function(o)); + } + + void mul(FourierMapType& m) const; + void mul_sqrt(FourierMapType& m) const; + void mul_inv(FourierMapType& m) const; + void mul_inv_sqrt(FourierMapType& m) const; + }; + + + template + void EuclidianSpectrum_1D::newRandomFourier(gsl_rng *rng, FourierMapType& out_map) const + { + typedef EuclidianFourierMapComplex MapT; + typedef typename EuclidianSpectrum_1D::ptr_map ptr_map; + typedef typename MapT::DimArray DimArray; + + MapT& rand_map = dynamic_cast(out_map); + + std::complex *d = rand_map.data(); + long idx; + const DimArray& dims = rand_map.getDims(); + const std::vector& delta_k = rand_map.get_delta_k(); + long plane_size; + bool alleven = rand_map.allDimensionsEven(); + double V = 1; + + for (int p = 0; p < delta_k.size(); p++) + V *= (2*M_PI/delta_k[p]); + + for (long p = 1; p < rand_map.size(); p++) + { + double A_k = std::sqrt(0.5*V*f(rand_map.get_K(p))); + d[p] = std::complex(gsl_ran_gaussian(rng, A_k), + gsl_ran_gaussian(rng, A_k)); + } + // Generate the mean value + d[0] = std::complex(gsl_ran_gaussian(rng, std::sqrt(V*f(0))), 0); + + if (!rand_map.firstDimensionEven()) + return; + + // Correct the Nyquist plane + idx = dims[0]-1; // Stick to the last element of the first dimension + d[idx] = std::complex(d[idx].real() + d[idx].imag(), 0); + // 1D is special case + if (dims.size() == 1) + return; + + plane_size = 1; + for (int q = 1; q < dims.size(); q++) + { + plane_size *= dims[q]; + } + + for (long p = 1; p < plane_size/2; p++) + { + long q = (p+1)*dims[0]-1; + long q2 = (plane_size-p+1)*dims[0]-1; + assert(q < plane_size*dims[0]); + assert(q2 < plane_size*dims[0]); + d[q] = conj(d[q2]); + } + + if (alleven) + { + long q = 0; + for (int i = dims.size()-1; i >= 1; i--) + q = dims[i]*q + dims[i]/2; + q += dims[0]-1; + d[q] = std::complex(d[q].real()+d[q].imag(),0); + } + } + + template + void EuclidianSpectrum_1D::mul(FourierMapType& m) const + { + EuclidianFourierMapComplex& m_c = dynamic_cast&>(m); + std::complex *d = m.data(); + + for (long p = 0; p < m_c.size(); p++) + d[p] *= f(m_c.get_K(p)); + } + + template + void EuclidianSpectrum_1D::mul_sqrt(FourierMapType& m) const + { + EuclidianFourierMapComplex& m_c = dynamic_cast&>(m); + std::complex *d = m.data(); + + for (long p = 0; p < m_c.size(); p++) + d[p] *= std::sqrt(f(m_c.get_K(p))); + } + + template + void EuclidianSpectrum_1D::mul_inv(FourierMapType& m) const + { + EuclidianFourierMapComplex& m_c = dynamic_cast&>(m); + std::complex *d = m.data(); + + for (long p = 0; p < m_c.size(); p++) + { + T A = f(m_c.get_K(p)); + if (A==0) + d[p] = 0; + else + d[p] /= A; + } + } + + template + void EuclidianSpectrum_1D::mul_inv_sqrt(FourierMapType& m) const + { + EuclidianFourierMapComplex& m_c = dynamic_cast&>(m); + std::complex *d = m.data(); + + for (long p = 0; p < m_c.size(); p++) + { + T A = std::sqrt(f(m_c.get_K(p))); + if (A == 0) + d[p] = 0; + else + d[p] /= A; + } + } + +}; + + +#endif diff --git a/src/fourier/details/euclidian_spectrum_1d_bin.hpp b/src/fourier/details/euclidian_spectrum_1d_bin.hpp new file mode 100644 index 0000000..ee6cd0b --- /dev/null +++ b/src/fourier/details/euclidian_spectrum_1d_bin.hpp @@ -0,0 +1,100 @@ +/*+ +This is CosmoTool (./src/fourier/details/euclidian_spectrum_1d_bin.hpp) -- Copyright (C) Guilhem Lavaux (2007-2013) + +guilhem.lavaux@gmail.com + +This software is a computer program whose purpose is to provide a toolbox for cosmological +data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...) + +This software is governed by the CeCILL license under French law and +abiding by the rules of distribution of free software. You can use, +modify and/ or redistribute the software under the terms of the CeCILL +license as circulated by CEA, CNRS and INRIA at the following URL +"http://www.cecill.info". + +As a counterpart to the access to the source code and rights to copy, +modify and redistribute granted by the license, users are provided only +with a limited warranty and the software's author, the holder of the +economic rights, and the successive licensors have only limited +liability. + +In this respect, the user's attention is drawn to the risks associated +with loading, using, modifying and/or developing or reproducing the +software by the user in light of its specific status of free software, +that may mean that it is complicated to manipulate, and that also +therefore means that it is reserved for developers and experienced +professionals having in-depth computer knowledge. Users are therefore +encouraged to load and test the software's suitability as regards their +requirements in conditions enabling the security of their systems and/or +data to be ensured and, more generally, to use and operate it in the +same conditions as regards security. + +The fact that you are presently reading this means that you have had +knowledge of the CeCILL license and that you accept its terms. ++*/ + +#ifndef __DETAILS_EUCLIDIAN_SPECTRUM_1D_BIN +#define __DETAILS_EUCLIDIAN_SPECTRUM_1D_BIN + +#include +#include + +namespace CosmoTool +{ + template + class EuclidianSpectrum_1D_Binned: public EuclidianSpectrum_1D + { + protected: + T *m_data; + long m_size; + T m_kmin, m_kmax, m_logdeltak; + std::vector m_dof; + public: + typedef typename SpectrumFunction::FourierMapType FourierMapType; + typedef typename SpectrumFunction::SpectrumFunctionPtr SpectrumFunctionPtr; + typedef boost::shared_ptr ptr_map; + + T interpolate_spectrum(T k) + { + T q = std::log(k/m_kmin)/m_logdeltak; + long ik = std::floor(q); + + if (ik >= m_size-1) + return m_data[m_size-1]; + else if (ik < 0) + return k/m_kmin*m_data[0]; + + return std::exp((q-ik)*m_data[ik+1] + (1-(q-ik))*m_data[ik]); + } + + EuclidianSpectrum_1D_Binned(int numBin, T kmin, T kmax) + : EuclidianSpectrum_1D(boost::bind(&EuclidianSpectrum_1D_Binned::interpolate_spectrum, this, _1)) + { + m_data = new T[numBin]; + m_kmin = kmin; + m_kmax = kmax; + m_size = numBin; + m_logdeltak = std::log(m_kmax/m_kmin); + } + + SpectrumFunctionPtr copy() const + { + EuclidianSpectrum_1D_Binned *s = new EuclidianSpectrum_1D_Binned(m_size, m_kmin, m_kmax); + std::copy(m_data, m_data+m_size, s->m_data); + return SpectrumFunctionPtr(s); + } + + void set_dof(std::vector& dof_array) + { + assert(m_size == dof_array.size()); + m_dof = dof_array; + } + + const T *data() const { return m_data; } + long size() const { return m_size; } + const T *dof() const { return &m_dof[0]; } + }; + +}; + +#endif diff --git a/src/fourier/details/euclidian_transform.hpp b/src/fourier/details/euclidian_transform.hpp new file mode 100644 index 0000000..527d38b --- /dev/null +++ b/src/fourier/details/euclidian_transform.hpp @@ -0,0 +1,216 @@ +/*+ +This is CosmoTool (./src/fourier/details/euclidian_transform.hpp) -- Copyright (C) Guilhem Lavaux (2007-2013) + +guilhem.lavaux@gmail.com + +This software is a computer program whose purpose is to provide a toolbox for cosmological +data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...) + +This software is governed by the CeCILL license under French law and +abiding by the rules of distribution of free software. You can use, +modify and/ or redistribute the software under the terms of the CeCILL +license as circulated by CEA, CNRS and INRIA at the following URL +"http://www.cecill.info". + +As a counterpart to the access to the source code and rights to copy, +modify and redistribute granted by the license, users are provided only +with a limited warranty and the software's author, the holder of the +economic rights, and the successive licensors have only limited +liability. + +In this respect, the user's attention is drawn to the risks associated +with loading, using, modifying and/or developing or reproducing the +software by the user in light of its specific status of free software, +that may mean that it is complicated to manipulate, and that also +therefore means that it is reserved for developers and experienced +professionals having in-depth computer knowledge. Users are therefore +encouraged to load and test the software's suitability as regards their +requirements in conditions enabling the security of their systems and/or +data to be ensured and, more generally, to use and operate it in the +same conditions as regards security. + +The fact that you are presently reading this means that you have had +knowledge of the CeCILL license and that you accept its terms. ++*/ + +#ifndef __DETAILS_EUCLIDIAN_TRANSFORM +#define __DETAILS_EUCLIDIAN_TRANSFORM + +namespace CosmoTool +{ + + template + class EuclidianFourierTransform: public FourierTransform + { + public: + typedef typename EuclidianFourierMapBase::DimArray DimArray; + private: + typedef FFTW_Calls calls; + EuclidianFourierMapReal *realMap; + EuclidianFourierMapComplex *fourierMap; + typename calls::plan_type m_analysis, m_synthesis; + double volume; + long N, Nc; + DimArray m_dims, m_dims_hc; + std::vector m_L; + public: + EuclidianFourierTransform(const DimArray& dims, const std::vector& L) + { + assert(L.size() == dims.size()); + std::vector dk(L.size()); + std::vector swapped_dims(dims.size()); + + m_dims = dims; + m_dims_hc = dims; + m_dims_hc[0] = dims[0]/2+1; + m_L = L; + + N = 1; + Nc = 1; + volume = 1; + for (int i = 0; i < dims.size(); i++) + { + N *= dims[i]; + Nc *= m_dims_hc[i]; + volume *= L[i]; + dk[i] = 2*M_PI/L[i]; + swapped_dims[dims.size()-1-i] = dims[i]; + } + + realMap = new EuclidianFourierMapReal( + boost::shared_ptr(calls::alloc_real(N), + std::ptr_fun(calls::free)), + m_dims); + fourierMap = new EuclidianFourierMapComplex( + boost::shared_ptr >((std::complex*)calls::alloc_complex(Nc), + std::ptr_fun(calls::free)), + dims[0], m_dims_hc, dk); + m_analysis = calls::plan_dft_r2c(dims.size(), &swapped_dims[0], + realMap->data(), (typename calls::complex_type *)fourierMap->data(), + FFTW_DESTROY_INPUT|FFTW_MEASURE); + m_synthesis = calls::plan_dft_c2r(dims.size(), &swapped_dims[0], + (typename calls::complex_type *)fourierMap->data(), realMap->data(), + FFTW_DESTROY_INPUT|FFTW_MEASURE); + } + + virtual ~EuclidianFourierTransform() + { + delete realMap; + delete fourierMap; + calls::destroy_plan(m_synthesis); + calls::destroy_plan(m_analysis); + } + + void synthesis() + { + calls::execute(m_synthesis); + realMap->scale(1/volume); + } + + void analysis() + { + calls::execute(m_analysis); + fourierMap->scale(volume/N); + } + + void synthesis_conjugate() + { + calls::execute(m_analysis); + fourierMap->scale(1/volume); + } + + void analysis_conjugate() + { + calls::execute(m_synthesis); + realMap->scale(volume/N); + } + + const FourierMap >& fourierSpace() const + { + return *fourierMap; + } + + FourierMap >& fourierSpace() + { + return *fourierMap; + } + + const FourierMap& realSpace() const + { + return *realMap; + } + + FourierMap& realSpace() + { + return *realMap; + } + + FourierTransform *mimick() const + { + return new EuclidianFourierTransform(m_dims, m_L); + } + }; + + template + class EuclidianFourierTransform_1d: public EuclidianFourierTransform + { + private: + template + static std::vector make_1d_vector(T2 a) + { + T2 arr[2] = { a}; + return std::vector(&arr[0],&arr[1]); + } + public: + EuclidianFourierTransform_1d(int Nx, double Lx) + : EuclidianFourierTransform(make_1d_vector(Nx), make_1d_vector(Lx)) + { + } + + virtual ~EuclidianFourierTransform_1d() {} + + }; + + template + class EuclidianFourierTransform_2d: public EuclidianFourierTransform + { + private: + template + static std::vector make_2d_vector(T2 a, T2 b) + { + T2 arr[2] = { a, b}; + return std::vector(&arr[0], &arr[2]); + } + public: + EuclidianFourierTransform_2d(int Nx, int Ny, double Lx, double Ly) + : EuclidianFourierTransform(make_2d_vector(Nx, Ny), make_2d_vector(Lx, Ly)) + { + } + + virtual ~EuclidianFourierTransform_2d() {} + + }; + + template + class EuclidianFourierTransform_3d: public EuclidianFourierTransform + { + private: + template + static std::vector make_3d_vector(T2 a, T2 b, T2 c) + { + T2 arr[3] = { a, b, c}; + return std::vector(&arr[0], &arr[3]); + } + + public: + EuclidianFourierTransform_3d(int Nx, int Ny, int Nz, double Lx, double Ly, double Lz) + : EuclidianFourierTransform(make_3d_vector(Nx, Ny, Nz), make_3d_vector(Lx, Ly, Lz)) + { + } + + virtual ~EuclidianFourierTransform_3d() {} + }; + +}; + +#endif diff --git a/src/fourier/details/healpix_alms.hpp b/src/fourier/details/healpix_alms.hpp new file mode 100644 index 0000000..fb438a8 --- /dev/null +++ b/src/fourier/details/healpix_alms.hpp @@ -0,0 +1,111 @@ +/*+ +This is CosmoTool (./src/fourier/details/healpix_alms.hpp) -- Copyright (C) Guilhem Lavaux (2007-2013) + +guilhem.lavaux@gmail.com + +This software is a computer program whose purpose is to provide a toolbox for cosmological +data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...) + +This software is governed by the CeCILL license under French law and +abiding by the rules of distribution of free software. You can use, +modify and/ or redistribute the software under the terms of the CeCILL +license as circulated by CEA, CNRS and INRIA at the following URL +"http://www.cecill.info". + +As a counterpart to the access to the source code and rights to copy, +modify and redistribute granted by the license, users are provided only +with a limited warranty and the software's author, the holder of the +economic rights, and the successive licensors have only limited +liability. + +In this respect, the user's attention is drawn to the risks associated +with loading, using, modifying and/or developing or reproducing the +software by the user in light of its specific status of free software, +that may mean that it is complicated to manipulate, and that also +therefore means that it is reserved for developers and experienced +professionals having in-depth computer knowledge. Users are therefore +encouraged to load and test the software's suitability as regards their +requirements in conditions enabling the security of their systems and/or +data to be ensured and, more generally, to use and operate it in the +same conditions as regards security. + +The fact that you are presently reading this means that you have had +knowledge of the CeCILL license and that you accept its terms. ++*/ + +#ifndef __COSMOTOOL_FOURIER_HEALPIX_DETAILS_ALM_HPP +#define __COSMOTOOL_FOURIER_HEALPIX_DETAILS_ALM_HPP + +namespace CosmoTool +{ + template + class HealpixFourierALM: public FourierMap > + { + private: + std::complex *alms; + long m_size; + long Lmax_, Mmax_, TVal_; + Eigen::aligned_allocator > alloc; + public: + typedef unsigned long LType; + + LType Lmax() const { return Lmax_; } + LType Mmax() const { return Mmax_; } + + LType Num_Alms() const + { + return ((Mmax_+1)*(Mmax_+2))/2 + (Mmax_+1)*(Lmax_-Mmax_); + } + + LType index_l0(LType m) const + { + return ((m*(TVal_-m))/2); + } + + LType index(LType l, LType m) const + { + return index_l0(m) + l; + } + + HealpixFourierALM(LType lmax, LType mmax) + : Lmax_(lmax), Mmax_(mmax), TVal_(2*lmax+1) + { + m_size = Num_Alms(); + alms = alloc.allocate(m_size); + } + + virtual ~HealpixFourierALM() + { + alloc.deallocate(alms, m_size); + } + + virtual const std::complex* data() const { return alms; } + virtual std::complex * data() { return alms;} + virtual long size() const { return m_size; } + + virtual FourierMap > *mimick() const + { + return new HealpixFourierALM(Lmax_, Mmax_); + } + + virtual std::complex dot_product(const FourierMap >& other) const + throw(std::bad_cast) + { + const HealpixFourierALM& mfm = dynamic_cast&>(other); + typedef typename FourierMap >::MapType MapType; + std::complex S; + + if (m_size != mfm.m_size) + throw std::bad_cast(); + + MapType m1(alms, m_size); + MapType m2(mfm.alms, mfm.m_size); + + S = (m1.block(0,0,1,Lmax_+1).conjugate() * m2.block(0,0,1,Lmax_+1)).sum(); + S += std::complex(2,0)*(m1.block(0,1+Lmax_,1,m_size-1-Lmax_).conjugate() * m2.block(0,1+Lmax_,1,m_size-1-Lmax_)).sum(); + return S; + } + }; +}; + +#endif diff --git a/src/fourier/details/healpix_map.hpp b/src/fourier/details/healpix_map.hpp new file mode 100644 index 0000000..a2a3f18 --- /dev/null +++ b/src/fourier/details/healpix_map.hpp @@ -0,0 +1,89 @@ +/*+ +This is CosmoTool (./src/fourier/details/healpix_map.hpp) -- Copyright (C) Guilhem Lavaux (2007-2013) + +guilhem.lavaux@gmail.com + +This software is a computer program whose purpose is to provide a toolbox for cosmological +data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...) + +This software is governed by the CeCILL license under French law and +abiding by the rules of distribution of free software. You can use, +modify and/ or redistribute the software under the terms of the CeCILL +license as circulated by CEA, CNRS and INRIA at the following URL +"http://www.cecill.info". + +As a counterpart to the access to the source code and rights to copy, +modify and redistribute granted by the license, users are provided only +with a limited warranty and the software's author, the holder of the +economic rights, and the successive licensors have only limited +liability. + +In this respect, the user's attention is drawn to the risks associated +with loading, using, modifying and/or developing or reproducing the +software by the user in light of its specific status of free software, +that may mean that it is complicated to manipulate, and that also +therefore means that it is reserved for developers and experienced +professionals having in-depth computer knowledge. Users are therefore +encouraged to load and test the software's suitability as regards their +requirements in conditions enabling the security of their systems and/or +data to be ensured and, more generally, to use and operate it in the +same conditions as regards security. + +The fact that you are presently reading this means that you have had +knowledge of the CeCILL license and that you accept its terms. ++*/ + +#ifndef __COSMOTOOL_FOURIER_HEALPIX_DETAILS_MAP_HPP +#define __COSMOTOOL_FOURIER_HEALPIX_DETAILS_MAP_HPP + +namespace CosmoTool +{ + + template + class HealpixFourierMap: public FourierMap + { + private: + T *m_data; + long Npix, m_Nside; + Eigen::aligned_allocator alloc; + public: + HealpixFourierMap(long nSide) + : Npix(12*nSide*nSide), m_Nside(nSide) + { + m_data = alloc.allocate(Npix); + } + + virtual ~HealpixFourierMap() + { + alloc.deallocate(m_data, Npix); + } + + long Nside() const { return m_Nside; } + virtual const T* data() const { return m_data; } + virtual T *data() { return m_data; } + virtual long size() const { return Npix; } + + virtual T dot_product(const FourierMap& other) const + throw(std::bad_cast) + { + typedef typename FourierMap::MapType MapType; + + const HealpixFourierMap& mfm = dynamic_cast&>(other); + if (Npix != mfm.size()) + throw std::bad_cast(); + + MapType m1(m_data, Npix); + MapType m2(mfm.m_data, mfm.Npix); + + return (m1*m2).sum(); + } + + virtual FourierMap *mimick() const + { + return new HealpixFourierMap(m_Nside); + } + }; + +}; + +#endif diff --git a/src/fourier/details/healpix_spectrum.hpp b/src/fourier/details/healpix_spectrum.hpp new file mode 100644 index 0000000..92f7085 --- /dev/null +++ b/src/fourier/details/healpix_spectrum.hpp @@ -0,0 +1,179 @@ +/*+ +This is CosmoTool (./src/fourier/details/healpix_spectrum.hpp) -- Copyright (C) Guilhem Lavaux (2007-2013) + +guilhem.lavaux@gmail.com + +This software is a computer program whose purpose is to provide a toolbox for cosmological +data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...) + +This software is governed by the CeCILL license under French law and +abiding by the rules of distribution of free software. You can use, +modify and/ or redistribute the software under the terms of the CeCILL +license as circulated by CEA, CNRS and INRIA at the following URL +"http://www.cecill.info". + +As a counterpart to the access to the source code and rights to copy, +modify and redistribute granted by the license, users are provided only +with a limited warranty and the software's author, the holder of the +economic rights, and the successive licensors have only limited +liability. + +In this respect, the user's attention is drawn to the risks associated +with loading, using, modifying and/or developing or reproducing the +software by the user in light of its specific status of free software, +that may mean that it is complicated to manipulate, and that also +therefore means that it is reserved for developers and experienced +professionals having in-depth computer knowledge. Users are therefore +encouraged to load and test the software's suitability as regards their +requirements in conditions enabling the security of their systems and/or +data to be ensured and, more generally, to use and operate it in the +same conditions as regards security. + +The fact that you are presently reading this means that you have had +knowledge of the CeCILL license and that you accept its terms. ++*/ + +#ifndef __COSMOTOOL_FOURIER_HEALPIX_DETAILS_SPECTRUM_HPP +#define __COSMOTOOL_FOURIER_HEALPIX_DETAILS_SPECTRUM_HPP + +namespace CosmoTool +{ + template + class HealpixSpectrum: public SpectrumFunction + { + protected: + std::vector cls; + int *m_dof; + public: + typedef typename SpectrumFunction::FourierMapType FourierMapType; + typedef boost::shared_ptr ptr_map; + typedef typename SpectrumFunction::SpectrumFunctionPtr SpectrumFunctionPtr; + typedef unsigned long LType; + + HealpixSpectrum(LType Lmax) + : cls(Lmax+1), m_dof(new int[Lmax+1]) + { + for (int l = 0; l <= Lmax; l++) + m_dof[l] = l + 1; + } + + T *data() { return &cls[0]; } + const T *data() const { return &cls[0]; } + const int *dof() const { return m_dof; } + + void set_dof(LType l, int dof) { m_dof[l] = dof; } + + LType Lmax() const { return cls.size()-1; } + long size() const { return cls.size(); } + + void newRandomFourier(gsl_rng *rng, FourierMapType& like_map) const; + + SpectrumFunctionPtr copy() const { + HealpixSpectrum *s = new HealpixSpectrum(Lmax()); + s->cls = cls; + return SpectrumFunctionPtr(s); + } + + void sqrt() { + std::transform(cls.begin(), cls.end(), cls.begin(), std::ptr_fun(std::sqrt)); + } + + void mul(FourierMapType& m) const; + void mul_sqrt(FourierMapType& m) const; + void mul_inv(FourierMapType& m) const; + void mul_inv_sqrt(FourierMapType& m) const; + }; + + template + void HealpixSpectrum::newRandomFourier(gsl_rng *rng, FourierMapType& out_map) const + { + HealpixFourierALM& alms = dynamic_cast&>(out_map); + long lmaxGen = std::min(cls.size()-1, alms.Lmax()); + std::complex *new_data = alms.data(); + + for (LType l = 0; l <= lmaxGen; l++) + { + double Al = std::sqrt(cls[l]); + + new_data[alms.index(l,0)] = gsl_ran_gaussian(rng, Al); + Al *= M_SQRT1_2; + for (LType m = 1; m <= std::min(l,alms.Mmax()); m++) + { + std::complex& c = new_data[alms.index(l,m)]; + c.real() = gsl_ran_gaussian(rng, Al); + c.imag() = gsl_ran_gaussian(rng, Al); + } + } + } + + template + void HealpixSpectrum::mul(FourierMapType& like_map) const + { + HealpixFourierALM& alms = dynamic_cast&>(like_map); + std::complex *data = alms.data(); + + for (LType l = 0; l <= alms.Lmax(); l++) + { + double Al = cls[l]; + + for (LType m = 0; m <= std::min(l,alms.Mmax()); m++) + { + data[alms.index(l,m)] *= Al; + } + } + } + + template + void HealpixSpectrum::mul_sqrt(FourierMapType& like_map) const + { + HealpixFourierALM& alms = dynamic_cast&>(like_map); + std::complex *data = alms.data(); + + for (LType l = 0; l <= alms.Lmax(); l++) + { + double Al = std::sqrt(cls[l]); + + for (LType m = 0; m <= std::min(l,alms.Mmax()); m++) + { + data[alms.index(l,m)] *= Al; + } + } + } + + template + void HealpixSpectrum::mul_inv(FourierMapType& like_map) const + { + HealpixFourierALM& alms = dynamic_cast&>(like_map); + std::complex *data = alms.data(); + + for (LType l = 0; l <= alms.Lmax(); l++) + { + double Al = (cls[l] <= 0) ? 0 : (1/cls[l]); + + for (LType m = 0; m <= std::min(l,alms.Mmax()); m++) + { + data[alms.index(l,m)] *= Al; + } + } + } + + template + void HealpixSpectrum::mul_inv_sqrt(FourierMapType& like_map) const + { + HealpixFourierALM& alms = dynamic_cast&>(like_map); + std::complex *data = alms.data(); + + for (LType l = 0; l <= alms.Lmax(); l++) + { + double Al = (cls[l] <= 0) ? 0 : std::sqrt(1/cls[l]); + + for (LType m = 0; m <= std::min(l,alms.Mmax()); m++) + { + data[alms.index(l,m)] *= Al; + } + } + } + +}; + +#endif diff --git a/src/fourier/details/healpix_transform.hpp b/src/fourier/details/healpix_transform.hpp new file mode 100644 index 0000000..6e89213 --- /dev/null +++ b/src/fourier/details/healpix_transform.hpp @@ -0,0 +1,132 @@ +/*+ +This is CosmoTool (./src/fourier/details/healpix_transform.hpp) -- Copyright (C) Guilhem Lavaux (2007-2013) + +guilhem.lavaux@gmail.com + +This software is a computer program whose purpose is to provide a toolbox for cosmological +data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...) + +This software is governed by the CeCILL license under French law and +abiding by the rules of distribution of free software. You can use, +modify and/ or redistribute the software under the terms of the CeCILL +license as circulated by CEA, CNRS and INRIA at the following URL +"http://www.cecill.info". + +As a counterpart to the access to the source code and rights to copy, +modify and redistribute granted by the license, users are provided only +with a limited warranty and the software's author, the holder of the +economic rights, and the successive licensors have only limited +liability. + +In this respect, the user's attention is drawn to the risks associated +with loading, using, modifying and/or developing or reproducing the +software by the user in light of its specific status of free software, +that may mean that it is complicated to manipulate, and that also +therefore means that it is reserved for developers and experienced +professionals having in-depth computer knowledge. Users are therefore +encouraged to load and test the software's suitability as regards their +requirements in conditions enabling the security of their systems and/or +data to be ensured and, more generally, to use and operate it in the +same conditions as regards security. + +The fact that you are presently reading this means that you have had +knowledge of the CeCILL license and that you accept its terms. ++*/ + +#ifndef __COSMOTOOL_FOURIER_HEALPIX_DETAILS_TRANSFORM_HPP +#define __COSMOTOOL_FOURIER_HEALPIX_DETAILS_TRANSFORM_HPP + +namespace CosmoTool +{ + + template struct HealpixJobHelper__ {}; + + template<> struct HealpixJobHelper__ + { enum {val=1}; }; + + template<> struct HealpixJobHelper__ + { enum {val=0}; }; + + + template + class HealpixFourierTransform: public FourierTransform + { + private: + sharp_alm_info *ainfo; + sharp_geom_info *ginfo; + HealpixFourierMap realMap; + HealpixFourierALM fourierMap; + int m_iterate; + public: + HealpixFourierTransform(long nSide, long Lmax, long Mmax, int iterate = 0) + : realMap(nSide), fourierMap(Lmax, Mmax), ainfo(0), ginfo(0), m_iterate(iterate) + { + sharp_make_healpix_geom_info (nSide, 1, &ginfo); + sharp_make_triangular_alm_info (Lmax, Mmax, 1, &ainfo); + } + + virtual ~HealpixFourierTransform() + { + sharp_destroy_geom_info(ginfo); + sharp_destroy_alm_info(ainfo); + } + + virtual const FourierMap >& fourierSpace() const { return fourierMap; } + + virtual FourierMap >& fourierSpace() { return fourierMap; } + + virtual const FourierMap& realSpace() const { return realMap; } + + virtual FourierMap& realSpace() { return realMap; } + + virtual FourierTransform *mimick() const + { + return new HealpixFourierTransform(realMap.Nside(), fourierMap.Lmax(), fourierMap.Mmax()); + } + + virtual void analysis() + { + void *aptr=reinterpret_cast(fourierMap.data()), *mptr=reinterpret_cast(realMap.data()); + + sharp_execute (SHARP_MAP2ALM, 0, 0, &aptr, &mptr, ginfo, ainfo, 1, + HealpixJobHelper__::val,0,0,0); + for (int i = 0; i < m_iterate; i++) + { + HealpixFourierMap tmp_map(realMap.Nside()); + void *tmp_ptr=reinterpret_cast(tmp_map.data()); + typename HealpixFourierMap::MapType m0 = tmp_map.eigen(); + typename HealpixFourierMap::MapType m1 = realMap.eigen(); + + sharp_execute (SHARP_ALM2MAP, 0, 0, &aptr, &tmp_ptr, ginfo, ainfo, 1, + HealpixJobHelper__::val,0,0,0); + m0 = m1 - m0; + sharp_execute (SHARP_MAP2ALM, 0, 1, &aptr, &tmp_ptr, ginfo, ainfo, 1, + HealpixJobHelper__::val,0,0,0); + } + } + + virtual void synthesis() + { + void *aptr=reinterpret_cast(fourierMap.data()), *mptr=reinterpret_cast(realMap.data()); + + sharp_execute (SHARP_ALM2MAP, 0, 0, &aptr, &mptr, ginfo, ainfo, 1, + HealpixJobHelper__::val,0,0,0); + } + + virtual void analysis_conjugate() + { + synthesis(); + realMap.scale(4*M_PI/realMap.size()); + } + + virtual void synthesis_conjugate() + { + analysis(); + fourierMap.scale(realMap.size()/(4*M_PI)); + } + + }; + +}; + +#endif diff --git a/src/fourier/details/healpix_utility.hpp b/src/fourier/details/healpix_utility.hpp new file mode 100644 index 0000000..8803ae9 --- /dev/null +++ b/src/fourier/details/healpix_utility.hpp @@ -0,0 +1,98 @@ +/*+ +This is CosmoTool (./src/fourier/details/healpix_utility.hpp) -- Copyright (C) Guilhem Lavaux (2007-2013) + +guilhem.lavaux@gmail.com + +This software is a computer program whose purpose is to provide a toolbox for cosmological +data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...) + +This software is governed by the CeCILL license under French law and +abiding by the rules of distribution of free software. You can use, +modify and/ or redistribute the software under the terms of the CeCILL +license as circulated by CEA, CNRS and INRIA at the following URL +"http://www.cecill.info". + +As a counterpart to the access to the source code and rights to copy, +modify and redistribute granted by the license, users are provided only +with a limited warranty and the software's author, the holder of the +economic rights, and the successive licensors have only limited +liability. + +In this respect, the user's attention is drawn to the risks associated +with loading, using, modifying and/or developing or reproducing the +software by the user in light of its specific status of free software, +that may mean that it is complicated to manipulate, and that also +therefore means that it is reserved for developers and experienced +professionals having in-depth computer knowledge. Users are therefore +encouraged to load and test the software's suitability as regards their +requirements in conditions enabling the security of their systems and/or +data to be ensured and, more generally, to use and operate it in the +same conditions as regards security. + +The fact that you are presently reading this means that you have had +knowledge of the CeCILL license and that you accept its terms. ++*/ + +#ifndef __COSMOTOOL_FOURIER_HEALPIX_DETAILS_SPECTRUM_HP +#define __COSMOTOOL_FOURIER_HEALPIX_DETAILS_SPECTRUM_HP + +namespace CosmoTool +{ + + template + class HealpixUtilityFunction: public MapUtilityFunction + { + public: + typedef typename MapUtilityFunction::Spectrum Spectrum; + typedef typename MapUtilityFunction::Spectrum_ptr Spectrum_ptr; + typedef typename MapUtilityFunction::FMap FMap; + typedef typename HealpixSpectrum::LType LType; + + Spectrum_ptr estimateSpectrumFromMap(const FMap& m) const + { + const HealpixFourierALM& alms = dynamic_cast&>(m); + LType Lmax = alms.Lmax(), Mmax = alms.Mmax(); + HealpixSpectrum *spectrum = new HealpixSpectrum(Lmax); + T *data_cls = spectrum->data(); + const std::complex *data_alms = alms.data(); + + // make an estimate of the cls + std::fill(data_cls, data_cls+Lmax+1, 0); + LType q = 0; + for (LType m = 0; m <= Mmax; ++m ) + { + for (LType l = m; l <= Lmax; ++l) + { + // Triangular storage + data_cls[l] += std::norm(data_alms[l+q]); + } + q += Lmax+1-m; + } + assert(q==alms.size()); + + for (LType l = 0; l <= Lmax; ++l) + { + int dof = 1 + std::min(l, Mmax); + spectrum->set_dof(l, dof); + data_cls[l] /= dof; + } + + return Spectrum_ptr(spectrum); + } + + Spectrum_ptr newSpectrumFromRaw(T *data, long size, + const Spectrum& like_spec) const + { + const HealpixSpectrum& in_spec = dynamic_cast&>(like_spec); + HealpixSpectrum *new_spectrum = new HealpixSpectrum(in_spec.Lmax()); + T *out_d = new_spectrum->data(); + + std::copy(data, data + min(size,new_spectrum->size()), out_d); + + return Spectrum_ptr(new_spectrum); + } + }; + +}; + +#endif diff --git a/src/fourier/euclidian.hpp b/src/fourier/euclidian.hpp index fac27f6..efc76e7 100644 --- a/src/fourier/euclidian.hpp +++ b/src/fourier/euclidian.hpp @@ -1,6 +1,42 @@ +/*+ +This is CosmoTool (./src/fourier/euclidian.hpp) -- Copyright (C) Guilhem Lavaux (2007-2013) + +guilhem.lavaux@gmail.com + +This software is a computer program whose purpose is to provide a toolbox for cosmological +data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...) + +This software is governed by the CeCILL license under French law and +abiding by the rules of distribution of free software. You can use, +modify and/ or redistribute the software under the terms of the CeCILL +license as circulated by CEA, CNRS and INRIA at the following URL +"http://www.cecill.info". + +As a counterpart to the access to the source code and rights to copy, +modify and redistribute granted by the license, users are provided only +with a limited warranty and the software's author, the holder of the +economic rights, and the successive licensors have only limited +liability. + +In this respect, the user's attention is drawn to the risks associated +with loading, using, modifying and/or developing or reproducing the +software by the user in light of its specific status of free software, +that may mean that it is complicated to manipulate, and that also +therefore means that it is reserved for developers and experienced +professionals having in-depth computer knowledge. Users are therefore +encouraged to load and test the software's suitability as regards their +requirements in conditions enabling the security of their systems and/or +data to be ensured and, more generally, to use and operate it in the +same conditions as regards security. + +The fact that you are presently reading this means that you have had +knowledge of the CeCILL license and that you accept its terms. ++*/ + #ifndef __COSMOTOOL_FOURIER_EUCLIDIAN_HPP #define __COSMOTOOL_FOURIER_EUCLIDIAN_HPP +#include #include #include #include @@ -8,400 +44,9 @@ #include "base_types.hpp" #include "fft/fftw_calls.hpp" #include "../algo.hpp" - -namespace CosmoTool -{ - template - class EuclidianSpectrum_1D: public SpectrumFunction - { - public: - typedef boost::function1 Function; - protected: - Function f; - public: - typedef typename SpectrumFunction::FourierMapType FourierMapType; - typedef boost::shared_ptr ptr_map; - - EuclidianSpectrum_1D(Function P) - : f(P) - { - } - - ptr_map newRandomFourier(gsl_rng *rng, const FourierMapType& like_map) const; - - void mul(FourierMap >& m) const; - }; - - template - class EuclidianFourierMapBase: public FourierMap - { - public: - typedef std::vector DimArray; - private: - boost::shared_ptr m_data; - DimArray m_dims; - long m_size; - public: - - EuclidianFourierMapBase(boost::shared_ptr indata, const DimArray& indims) - { - m_data = indata; - m_dims = indims; - m_size = 1; - for (int i = 0; i < m_dims.size(); i++) - m_size *= m_dims[i]; - } - - virtual ~EuclidianFourierMapBase() - { - } - - const DimArray& getDims() const { return m_dims; } - - virtual const T *data() const { return m_data.get(); } - virtual T *data() { return m_data.get(); } - virtual long size() const { return m_size; } - - virtual FourierMap *copy() const - { - FourierMap *m = this->mimick(); - m->eigen() = this->eigen(); - return m; - } - - }; - - template - class EuclidianFourierMapReal: public EuclidianFourierMapBase - { - public: - typedef typename EuclidianFourierMapBase::DimArray DimArray; - - EuclidianFourierMapReal(boost::shared_ptr indata, const DimArray& indims) - : EuclidianFourierMapBase(indata, indims) - {} - - virtual FourierMap *mimick() const - { - return new EuclidianFourierMapReal( - boost::shared_ptr((T *)fftw_malloc(sizeof(T)*this->size()), - std::ptr_fun(fftw_free)), - this->getDims()); - } - - virtual T dot_product(const FourierMap& other) const - throw(std::bad_cast) - { - const EuclidianFourierMapReal& m2 = dynamic_cast&>(other); - if (this->size() != m2.size()) - throw std::bad_cast(); - - return (this->eigen()*m2.eigen()).sum(); - } - }; - - template - class EuclidianFourierMapComplex: public EuclidianFourierMapBase > - { - protected: - typedef boost::shared_ptr > ptr_t; - std::vector delta_k; - long plane_size; - public: - typedef typename EuclidianFourierMapBase >::DimArray DimArray; - - EuclidianFourierMapComplex(ptr_t indata, - const DimArray& indims, - const std::vector& dk) - : EuclidianFourierMapBase >(indata, indims), delta_k(dk) - { - assert(dk.size() == indims.size()); - plane_size = 1; - for (int q = 1; q < indims.size(); q++) - plane_size *= indims[q]; - } - - virtual FourierMap > *mimick() const - { - return - new EuclidianFourierMapComplex( - ptr_t((std::complex *) - fftw_malloc(sizeof(std::complex)*this->size()), - std::ptr_fun(fftw_free)), - this->getDims(), - this->delta_k); - } - - template - double get_K(const Array& ik) - { - const DimArray& dims = this->getDims(); - assert(ik.size() == dims.size()); - double k2 = 0; - for (int q = 0; q < ik.size(); q++) - { - int dk = ik; - - if (dk > dims[q]/2) - dk = dk - dims[q]; - - k2 += CosmoTool::square(delta_k[q]*dk); - } - return std::sqrt(k2); - } - - double get_K(long p) - { - const DimArray& dims = this->getDims(); - DimArray d(delta_k.size()); - for (int q = 0; q < d.size(); q++) - { - d[q] = p%dims[q]; - p = (p-d[q])/dims[q]; - } - return get_K(d); - } - - virtual std::complex dot_product(const FourierMap >& other) const - throw(std::bad_cast) - { - const EuclidianFourierMapComplex& m2 = dynamic_cast&>(other); - if (this->size() != m2.size()) - throw std::bad_cast(); - - const std::complex *d1 = this->data(); - const std::complex *d2 = m2.data(); - const DimArray& dims = this->getDims(); - int N0 = dims[0]; - std::complex result = 0; - - for (long q0 = 1; q0 < N0-1; q0++) - { - for (long p = 0; p < plane_size; p++) - { - long idx = q0+N0*p; - assert(idx < this->size()); - result += 2*(conj(d1[idx]) * d2[idx]).real(); - } - } - for (long p = 0; p < plane_size; p++) - { - long q0 = N0*p, q1 = (p+1)*N0-1; - result += conj(d1[q0]) * d2[q0]; - result += conj(d1[q1]) * d2[q1]; - } - return result; - } - - }; - - template - class EuclidianFourierTransform: public FourierTransform - { - public: - typedef typename EuclidianFourierMapBase::DimArray DimArray; - private: - typedef FFTW_Calls calls; - EuclidianFourierMapReal *realMap; - EuclidianFourierMapComplex *fourierMap; - typename calls::plan_type m_analysis, m_synthesis; - double volume; - long N; - DimArray m_dims, m_dims_hc; - std::vector m_L; - public: - EuclidianFourierTransform(const DimArray& dims, const std::vector& L) - { - assert(L.size() == dims.size()); - std::vector dk(L.size()); - - m_dims = dims; - m_dims_hc = dims; - m_dims_hc[0] = dims[0]/2+1; - m_L = L; - - N = 1; - volume = 1; - for (int i = 0; i < dims.size(); i++) - { - N *= dims[i]; - volume *= L[i]; - dk[i] = 2*M_PI/L[i]; - } - - realMap = new EuclidianFourierMapReal( - boost::shared_ptr(calls::alloc_real(N), - std::ptr_fun(calls::free)), - m_dims); - fourierMap = new EuclidianFourierMapComplex( - boost::shared_ptr >((std::complex*)calls::alloc_complex(N), - std::ptr_fun(calls::free)), - m_dims_hc, dk); - m_analysis = calls::plan_dft_r2c(dims.size(), &dims[0], - realMap->data(), (typename calls::complex_type *)fourierMap->data(), - FFTW_MEASURE); - m_synthesis = calls::plan_dft_c2r(dims.size(), &dims[0], - (typename calls::complex_type *)fourierMap->data(), realMap->data(), - FFTW_MEASURE); - } - - virtual ~EuclidianFourierTransform() - { - delete realMap; - delete fourierMap; - calls::destroy_plan(m_synthesis); - calls::destroy_plan(m_analysis); - } - - void synthesis() - { - calls::execute(m_synthesis); - realMap->scale(1/volume); - } - - void analysis() - { - calls::execute(m_analysis); - fourierMap->scale(volume/N); - } - - void synthesis_conjugate() - { - calls::execute(m_analysis); - fourierMap->scale(1/volume); - } - - void analysis_conjugate() - { - calls::execute(m_synthesis); - realMap->scale(volume/N); - } - - const FourierMap >& fourierSpace() const - { - return *fourierMap; - } - - FourierMap >& fourierSpace() - { - return *fourierMap; - } - - const FourierMap& realSpace() const - { - return *realMap; - } - - FourierMap& realSpace() - { - return *realMap; - } - - FourierTransform *mimick() const - { - return new EuclidianFourierTransform(m_dims, m_L); - } - }; - - template - class EuclidianFourierTransform_2d: public EuclidianFourierTransform - { - private: - template - static std::vector make_2d_vector(T2 a, T2 b) - { - T2 arr[2] = { a, b}; - return std::vector(&arr[0], &arr[2]); - } - public: - EuclidianFourierTransform_2d(int Nx, int Ny, double Lx, double Ly) - : EuclidianFourierTransform(make_2d_vector(Nx, Ny), make_2d_vector(Lx, Ly)) - { - } - - virtual ~EuclidianFourierTransform_2d() {} - - }; - - template - class EuclidianFourierTransform_3d: public EuclidianFourierTransform - { - private: - template - static std::vector make_3d_vector(T2 a, T2 b, T2 c) - { - T2 arr[2] = { a, b, c}; - return std::vector(&arr[0], &arr[3]); - } - - public: - EuclidianFourierTransform_3d(int Nx, int Ny, int Nz, double Lx, double Ly, double Lz) - : EuclidianFourierTransform(make_3d_vector(Nx, Ny, Nz), make_3d_vector(Lx, Ly, Lz)) - { - } - - virtual ~EuclidianFourierTransform_3d() {} - }; - - - template - typename EuclidianSpectrum_1D::ptr_map - EuclidianSpectrum_1D::newRandomFourier(gsl_rng *rng, const FourierMapType& like_map) const - { - typedef EuclidianFourierMapComplex MapT; - typedef typename MapT::DimArray DimArray; - - MapT& m_c = dynamic_cast(like_map); - MapT *rand_map = m_c.mimick(); - std::complex *d = rand_map->data(); - long idx; - const DimArray& dims = rand_map->getDims(); - long plane_size; - - for (long p = 1; p < m_c.size(); p++) - { - double A_k = std::sqrt(0.5*f(rand_map->get_K(p))); - d[p] = std::complex(gsl_ran_gaussian(rng, A_k), - gsl_ran_gaussian(rng, A_k)); - } - // Generate the mean value - d[0] = std::complex(std::sqrt(f(0)), 0); - // Correct the Nyquist plane - idx = dims[0]-1; // Stick to the last element of the first dimension - // 1D is special case - if (dims.size() == 1) - { - d[idx] = std::complex(d[idx].real() + d[idx].imag(), 0); - return boost::shared_ptr::FourierMapType>(rand_map); - } - - plane_size = 1; - for (int q = 1; q < dims.size(); q++) - { - plane_size *= dims[q]; - } - - for (long p = 1; p < plane_size/2; p++) - { - long q = (p+1)*dims[0]-1; - long q2 = (plane_size-p+1)*dims[0]-1; - assert(q < plane_size*dims[0]); - assert(q2 < plane_size*dims[0]); - d[q] = conj(d[q2]); - } - long q = dims[0]; - d[q] = std::complex(d[q].real() + d[q].imag()); - } - - template - void EuclidianSpectrum_1D::mul(FourierMap >& m) const - { - EuclidianFourierMapComplex& m_c = dynamic_cast&>(m); - std::complex *d = m.data(); - - for (long p = 0; p < m_c.size(); p++) - d[p] *= f(m.get_K(p)); - } -}; +#include "details/euclidian_maps.hpp" +#include "details/euclidian_spectrum_1d.hpp" +#include "details/euclidian_spectrum_1d_bin.hpp" +#include "details/euclidian_transform.hpp" #endif diff --git a/src/fourier/fft/fftw_calls.hpp b/src/fourier/fft/fftw_calls.hpp index 0106055..c811728 100644 --- a/src/fourier/fft/fftw_calls.hpp +++ b/src/fourier/fft/fftw_calls.hpp @@ -1,3 +1,38 @@ +/*+ +This is CosmoTool (./src/fourier/fft/fftw_calls.hpp) -- Copyright (C) Guilhem Lavaux (2007-2013) + +guilhem.lavaux@gmail.com + +This software is a computer program whose purpose is to provide a toolbox for cosmological +data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...) + +This software is governed by the CeCILL license under French law and +abiding by the rules of distribution of free software. You can use, +modify and/ or redistribute the software under the terms of the CeCILL +license as circulated by CEA, CNRS and INRIA at the following URL +"http://www.cecill.info". + +As a counterpart to the access to the source code and rights to copy, +modify and redistribute granted by the license, users are provided only +with a limited warranty and the software's author, the holder of the +economic rights, and the successive licensors have only limited +liability. + +In this respect, the user's attention is drawn to the risks associated +with loading, using, modifying and/or developing or reproducing the +software by the user in light of its specific status of free software, +that may mean that it is complicated to manipulate, and that also +therefore means that it is reserved for developers and experienced +professionals having in-depth computer knowledge. Users are therefore +encouraged to load and test the software's suitability as regards their +requirements in conditions enabling the security of their systems and/or +data to be ensured and, more generally, to use and operate it in the +same conditions as regards security. + +The fact that you are presently reading this means that you have had +knowledge of the CeCILL license and that you accept its terms. ++*/ + #ifndef __FFTW_UNIFIED_CALLS_HPP #define __FFTW_UNIFIED_CALLS_HPP @@ -28,8 +63,8 @@ public: \ typedef prefix ## _complex complex_type; \ typedef prefix ## _plan plan_type; \ \ - static complex_type *alloc_complex(int N) { return prefix ## _alloc_complex(N); } \ - static real_type *alloc_real(int N) { return prefix ## _alloc_real(N); } \ + static complex_type *alloc_complex(size_t N) { return prefix ## _alloc_complex(N); } \ + static real_type *alloc_real(size_t N) { return prefix ## _alloc_real(N); } \ static void free(void *p) { fftw_free(p); } \ \ static void execute(plan_type p) { prefix ## _execute(p); } \ diff --git a/src/fourier/healpix.hpp b/src/fourier/healpix.hpp index 8aac379..e8b8b2d 100644 --- a/src/fourier/healpix.hpp +++ b/src/fourier/healpix.hpp @@ -1,3 +1,38 @@ +/*+ +This is CosmoTool (./src/fourier/healpix.hpp) -- Copyright (C) Guilhem Lavaux (2007-2013) + +guilhem.lavaux@gmail.com + +This software is a computer program whose purpose is to provide a toolbox for cosmological +data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...) + +This software is governed by the CeCILL license under French law and +abiding by the rules of distribution of free software. You can use, +modify and/ or redistribute the software under the terms of the CeCILL +license as circulated by CEA, CNRS and INRIA at the following URL +"http://www.cecill.info". + +As a counterpart to the access to the source code and rights to copy, +modify and redistribute granted by the license, users are provided only +with a limited warranty and the software's author, the holder of the +economic rights, and the successive licensors have only limited +liability. + +In this respect, the user's attention is drawn to the risks associated +with loading, using, modifying and/or developing or reproducing the +software by the user in light of its specific status of free software, +that may mean that it is complicated to manipulate, and that also +therefore means that it is reserved for developers and experienced +professionals having in-depth computer knowledge. Users are therefore +encouraged to load and test the software's suitability as regards their +requirements in conditions enabling the security of their systems and/or +data to be ensured and, more generally, to use and operate it in the +same conditions as regards security. + +The fact that you are presently reading this means that you have had +knowledge of the CeCILL license and that you accept its terms. ++*/ + #ifndef __COSMOTOOL_FOURIER_HEALPIX_HPP #define __COSMOTOOL_FOURIER_HEALPIX_HPP @@ -12,257 +47,13 @@ #include #include #include +#include +#include -namespace CosmoTool -{ - template - class HealpixSpectrum: public SpectrumFunction - { - protected: - std::vector cls; - public: - typedef typename SpectrumFunction::FourierMapType FourierMapType; - typedef boost::shared_ptr ptr_map; - - - HealpixSpectrum(long Lmax) - : cls(Lmax+1) {} - - T *data() { return &cls[0]; } - const T *data() const { return &cls[0]; } - long size() const { return cls.size(); } - - ptr_map newRandomFourier(gsl_rng *rng, const FourierMapType& like_map) const = 0; - }; - - template - class HealpixFourierMap: public FourierMap - { - private: - T *m_data; - long Npix, m_Nside; - Eigen::aligned_allocator alloc; - public: - HealpixFourierMap(long nSide) - : Npix(12*nSide*nSide), m_Nside(nSide) - { - m_data = alloc.allocate(Npix); - } - - virtual ~HealpixFourierMap() - { - alloc.deallocate(m_data, Npix); - } - - long Nside() const { return m_Nside; } - virtual const T* data() const { return m_data; } - virtual T *data() { return m_data; } - virtual long size() const { return Npix; } - - virtual T dot_product(const FourierMap& other) const - throw(std::bad_cast) - { - typedef typename FourierMap::MapType MapType; - - const HealpixFourierMap& mfm = dynamic_cast&>(other); - if (Npix != mfm.size()) - throw std::bad_cast(); - - MapType m1(m_data, Npix); - MapType m2(mfm.m_data, mfm.Npix); - - return (m1*m2).sum(); - } - - virtual FourierMap *mimick() const - { - return new HealpixFourierMap(m_Nside); - } - }; - - - template - class HealpixFourierALM: public FourierMap > - { - private: - std::complex *alms; - long m_size; - long Lmax_, Mmax_, TVal_; - Eigen::aligned_allocator > alloc; - public: - typedef unsigned long LType; - - LType Lmax() const { return Lmax_; } - LType Mmax() const { return Mmax_; } - - LType Num_Alms() const - { - return ((Mmax_+1)*(Mmax_+2))/2 + (Mmax_+1)*(Lmax_-Mmax_); - } - - LType index_l0(LType m) const - { - return ((m*(TVal_-m))/2); - } - - LType index(LType l, LType m) const - { - return index_l0(m) + l; - } - - HealpixFourierALM(LType lmax, LType mmax) - : Lmax_(lmax), Mmax_(mmax), TVal_(2*lmax+1) - { - m_size = Num_Alms(); - alms = alloc.allocate(m_size); - } - - virtual ~HealpixFourierALM() - { - alloc.deallocate(alms, m_size); - } - - virtual const std::complex* data() const { return alms; } - virtual std::complex * data() { return alms;} - virtual long size() const { return m_size; } - - virtual FourierMap > *mimick() const - { - return new HealpixFourierALM(Lmax_, Mmax_); - } - - virtual std::complex dot_product(const FourierMap >& other) const - throw(std::bad_cast) - { - const HealpixFourierALM& mfm = dynamic_cast&>(other); - typedef typename FourierMap >::MapType MapType; - std::complex S; - - if (m_size != mfm.m_size) - throw std::bad_cast(); - - MapType m1(alms, m_size); - MapType m2(mfm.alms, mfm.m_size); - - S = (m1.block(0,0,1,Lmax_+1).conjugate() * m2.block(0,0,1,Lmax_+1)).sum(); - S += std::complex(2,0)*(m1.block(0,1+Lmax_,1,m_size-1-Lmax_).conjugate() * m2.block(0,1+Lmax_,1,m_size-1-Lmax_)).sum(); - return S; - } - }; - - template struct HealpixJobHelper__ {}; - - template<> struct HealpixJobHelper__ - { enum {val=1}; }; - - template<> struct HealpixJobHelper__ - { enum {val=0}; }; - - - template - class HealpixFourierTransform: public FourierTransform - { - private: - sharp_alm_info *ainfo; - sharp_geom_info *ginfo; - HealpixFourierMap realMap; - HealpixFourierALM fourierMap; - int m_iterate; - public: - HealpixFourierTransform(long nSide, long Lmax, long Mmax, int iterate = 0) - : realMap(nSide), fourierMap(Lmax, Mmax), ainfo(0), ginfo(0), m_iterate(iterate) - { - sharp_make_healpix_geom_info (nSide, 1, &ginfo); - sharp_make_triangular_alm_info (Lmax, Mmax, 1, &ainfo); - } - - virtual ~HealpixFourierTransform() - { - sharp_destroy_geom_info(ginfo); - sharp_destroy_alm_info(ainfo); - } - - virtual const FourierMap >& fourierSpace() const { return fourierMap; } - - virtual FourierMap >& fourierSpace() { return fourierMap; } - - virtual const FourierMap& realSpace() const { return realMap; } - - virtual FourierMap& realSpace() { return realMap; } - - virtual FourierTransform *mimick() const - { - return new HealpixFourierTransform(realMap.Nside(), fourierMap.Lmax(), fourierMap.Mmax()); - } - - virtual void analysis() - { - void *aptr=reinterpret_cast(fourierMap.data()), *mptr=reinterpret_cast(realMap.data()); - - sharp_execute (SHARP_MAP2ALM, 0, 0, &aptr, &mptr, ginfo, ainfo, 1, - HealpixJobHelper__::val,0,0,0); - for (int i = 0; i < m_iterate; i++) - { - HealpixFourierMap tmp_map(realMap.Nside()); - void *tmp_ptr=reinterpret_cast(tmp_map.data()); - typename HealpixFourierMap::MapType m0 = tmp_map.eigen(); - typename HealpixFourierMap::MapType m1 = realMap.eigen(); - - sharp_execute (SHARP_ALM2MAP, 0, 0, &aptr, &tmp_ptr, ginfo, ainfo, 1, - HealpixJobHelper__::val,0,0,0); - m0 = m1 - m0; - sharp_execute (SHARP_MAP2ALM, 0, 1, &aptr, &tmp_ptr, ginfo, ainfo, 1, - HealpixJobHelper__::val,0,0,0); - } - } - - virtual void synthesis() - { - void *aptr=reinterpret_cast(fourierMap.data()), *mptr=reinterpret_cast(realMap.data()); - - sharp_execute (SHARP_ALM2MAP, 0, 0, &aptr, &mptr, ginfo, ainfo, 1, - HealpixJobHelper__::val,0,0,0); - } - - virtual void analysis_conjugate() - { - synthesis(); - realMap.scale(4*M_PI/realMap.size()); - } - - virtual void synthesis_conjugate() - { - analysis(); - fourierMap.scale(realMap.size()/(4*M_PI)); - } - - }; - - template - typename HealpixSpectrum::ptr_map - HealpixSpectrum::newRandomFourier(gsl_rng *rng, const FourierMapType& like_map) const - { - const HealpixFourierALM& alms = dynamic_cast&>(like_map); - HealpixFourierALM *new_alms; - ptr_map r(new_alms = new HealpixFourierALM(alms.Lmax(), alms.Mmax())); - long lmaxGen = std::min(cls.size()-1, alms.Lmax()); - std::complex *new_data = new_alms->data(); - - for (long l = 0; l < lmaxGen; l++) - { - double Al = std::sqrt(cls[l]); - - new_data[alms.index(l,0)] = gsl_ran_gaussian(rng, Al); - Al *= M_SQRT1_2; - for (long m = 1; m < alms.Mmax(); m++) - { - std::complex& c = new_data[alms.index(l,m)]; - c.real() = gsl_ran_gaussian(rng, Al); - c.imag() = gsl_ran_gaussian(rng, Al); - } - } - return r; - } -}; +#include "details/healpix_map.hpp" +#include "details/healpix_alms.hpp" +#include "details/healpix_transform.hpp" +#include "details/healpix_spectrum.hpp" +#include "details/healpix_utility.hpp" #endif diff --git a/src/growthFactor.cpp b/src/growthFactor.cpp index 485d876..b718d4e 100644 --- a/src/growthFactor.cpp +++ b/src/growthFactor.cpp @@ -1,3 +1,38 @@ +/*+ +This is CosmoTool (./src/growthFactor.cpp) -- Copyright (C) Guilhem Lavaux (2007-2013) + +guilhem.lavaux@gmail.com + +This software is a computer program whose purpose is to provide a toolbox for cosmological +data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...) + +This software is governed by the CeCILL license under French law and +abiding by the rules of distribution of free software. You can use, +modify and/ or redistribute the software under the terms of the CeCILL +license as circulated by CEA, CNRS and INRIA at the following URL +"http://www.cecill.info". + +As a counterpart to the access to the source code and rights to copy, +modify and redistribute granted by the license, users are provided only +with a limited warranty and the software's author, the holder of the +economic rights, and the successive licensors have only limited +liability. + +In this respect, the user's attention is drawn to the risks associated +with loading, using, modifying and/or developing or reproducing the +software by the user in light of its specific status of free software, +that may mean that it is complicated to manipulate, and that also +therefore means that it is reserved for developers and experienced +professionals having in-depth computer knowledge. Users are therefore +encouraged to load and test the software's suitability as regards their +requirements in conditions enabling the security of their systems and/or +data to be ensured and, more generally, to use and operate it in the +same conditions as regards security. + +The fact that you are presently reading this means that you have had +knowledge of the CeCILL license and that you accept its terms. ++*/ + #include #include #include "interpolate.hpp" diff --git a/src/growthFactor.hpp b/src/growthFactor.hpp index 2643fd0..1a455df 100644 --- a/src/growthFactor.hpp +++ b/src/growthFactor.hpp @@ -1,3 +1,38 @@ +/*+ +This is CosmoTool (./src/growthFactor.hpp) -- Copyright (C) Guilhem Lavaux (2007-2013) + +guilhem.lavaux@gmail.com + +This software is a computer program whose purpose is to provide a toolbox for cosmological +data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...) + +This software is governed by the CeCILL license under French law and +abiding by the rules of distribution of free software. You can use, +modify and/ or redistribute the software under the terms of the CeCILL +license as circulated by CEA, CNRS and INRIA at the following URL +"http://www.cecill.info". + +As a counterpart to the access to the source code and rights to copy, +modify and redistribute granted by the license, users are provided only +with a limited warranty and the software's author, the holder of the +economic rights, and the successive licensors have only limited +liability. + +In this respect, the user's attention is drawn to the risks associated +with loading, using, modifying and/or developing or reproducing the +software by the user in light of its specific status of free software, +that may mean that it is complicated to manipulate, and that also +therefore means that it is reserved for developers and experienced +professionals having in-depth computer knowledge. Users are therefore +encouraged to load and test the software's suitability as regards their +requirements in conditions enabling the security of their systems and/or +data to be ensured and, more generally, to use and operate it in the +same conditions as regards security. + +The fact that you are presently reading this means that you have had +knowledge of the CeCILL license and that you accept its terms. ++*/ + #ifndef COSMO_GROWTH_FACTOR_HPP #define COSMO_GROWTH_FACTOR_HPP diff --git a/src/h5_readFlash.cpp b/src/h5_readFlash.cpp index 0c78f58..e460bdb 100644 --- a/src/h5_readFlash.cpp +++ b/src/h5_readFlash.cpp @@ -1,3 +1,11 @@ +/*+ +!! + +This file has been developped by P. M. Sutter + +!! + ++*/ /* This file contains the functions that read the data from the HDF5 file * The functions accept the PARAMESH data through arguments. */ @@ -185,7 +193,7 @@ void h5_read_flash3_particles (H5File* file, float *vel1, float *vel2, float *vel3, - int *id) + long *id) { herr_t status; @@ -341,7 +349,7 @@ void h5_read_flash3_particles (H5File* file, if (id) { for(p=0; p < (pcount); p++) { - id[p+poffset] = (int) *(partBuffer+iptag-1+p*numProps); + id[p+poffset] = (long) *(partBuffer+iptag-1+p*numProps); } } if (pos1 && pos2 && pos3) { diff --git a/src/h5_readFlash.hpp b/src/h5_readFlash.hpp index 5187232..abcf9aa 100644 --- a/src/h5_readFlash.hpp +++ b/src/h5_readFlash.hpp @@ -1,3 +1,12 @@ +/*+ + +!!! NOTE !!! + +This file has been developped by P. M. Sutter. + +!!!! + ++*/ /* This file contains the functions that read the data from the HDF5 file * The functions accept the PARAMESH data through arguments. */ @@ -29,7 +38,7 @@ void h5_read_flash3_particles (H5File* file, float *vel1, float *vel2, float *vel3, - int *id); + long *id); void h5_read_flash3_header_info(H5File* file, double* time, /* simulation time */ diff --git a/src/hdf5_flash.h b/src/hdf5_flash.h index aab54e6..bf53f36 100644 --- a/src/hdf5_flash.h +++ b/src/hdf5_flash.h @@ -1,3 +1,10 @@ +/*+ +!! + +This particular file has been developped by P. M. Sutter + +!! ++*/ /* general header file for the HDF 5 IO in FLASH */ diff --git a/src/interpolate.cpp b/src/interpolate.cpp index 9a9fc70..ec6575c 100644 --- a/src/interpolate.cpp +++ b/src/interpolate.cpp @@ -1,3 +1,38 @@ +/*+ +This is CosmoTool (./src/interpolate.cpp) -- Copyright (C) Guilhem Lavaux (2007-2013) + +guilhem.lavaux@gmail.com + +This software is a computer program whose purpose is to provide a toolbox for cosmological +data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...) + +This software is governed by the CeCILL license under French law and +abiding by the rules of distribution of free software. You can use, +modify and/ or redistribute the software under the terms of the CeCILL +license as circulated by CEA, CNRS and INRIA at the following URL +"http://www.cecill.info". + +As a counterpart to the access to the source code and rights to copy, +modify and redistribute granted by the license, users are provided only +with a limited warranty and the software's author, the holder of the +economic rights, and the successive licensors have only limited +liability. + +In this respect, the user's attention is drawn to the risks associated +with loading, using, modifying and/or developing or reproducing the +software by the user in light of its specific status of free software, +that may mean that it is complicated to manipulate, and that also +therefore means that it is reserved for developers and experienced +professionals having in-depth computer knowledge. Users are therefore +encouraged to load and test the software's suitability as regards their +requirements in conditions enabling the security of their systems and/or +data to be ensured and, more generally, to use and operate it in the +same conditions as regards security. + +The fact that you are presently reading this means that you have had +knowledge of the CeCILL license and that you accept its terms. ++*/ + #include #include #include diff --git a/src/interpolate.hpp b/src/interpolate.hpp index a0ffacc..2e30e29 100644 --- a/src/interpolate.hpp +++ b/src/interpolate.hpp @@ -1,3 +1,38 @@ +/*+ +This is CosmoTool (./src/interpolate.hpp) -- Copyright (C) Guilhem Lavaux (2007-2013) + +guilhem.lavaux@gmail.com + +This software is a computer program whose purpose is to provide a toolbox for cosmological +data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...) + +This software is governed by the CeCILL license under French law and +abiding by the rules of distribution of free software. You can use, +modify and/ or redistribute the software under the terms of the CeCILL +license as circulated by CEA, CNRS and INRIA at the following URL +"http://www.cecill.info". + +As a counterpart to the access to the source code and rights to copy, +modify and redistribute granted by the license, users are provided only +with a limited warranty and the software's author, the holder of the +economic rights, and the successive licensors have only limited +liability. + +In this respect, the user's attention is drawn to the risks associated +with loading, using, modifying and/or developing or reproducing the +software by the user in light of its specific status of free software, +that may mean that it is complicated to manipulate, and that also +therefore means that it is reserved for developers and experienced +professionals having in-depth computer knowledge. Users are therefore +encouraged to load and test the software's suitability as regards their +requirements in conditions enabling the security of their systems and/or +data to be ensured and, more generally, to use and operate it in the +same conditions as regards security. + +The fact that you are presently reading this means that you have had +knowledge of the CeCILL license and that you accept its terms. ++*/ + #ifndef __CTOOL_INTERPOLATE_HPP #define __CTOOL_INTERPOLATE_HPP diff --git a/src/interpolate3d.hpp b/src/interpolate3d.hpp index 8c165ff..0f66d02 100644 --- a/src/interpolate3d.hpp +++ b/src/interpolate3d.hpp @@ -1,3 +1,38 @@ +/*+ +This is CosmoTool (./src/interpolate3d.hpp) -- Copyright (C) Guilhem Lavaux (2007-2013) + +guilhem.lavaux@gmail.com + +This software is a computer program whose purpose is to provide a toolbox for cosmological +data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...) + +This software is governed by the CeCILL license under French law and +abiding by the rules of distribution of free software. You can use, +modify and/ or redistribute the software under the terms of the CeCILL +license as circulated by CEA, CNRS and INRIA at the following URL +"http://www.cecill.info". + +As a counterpart to the access to the source code and rights to copy, +modify and redistribute granted by the license, users are provided only +with a limited warranty and the software's author, the holder of the +economic rights, and the successive licensors have only limited +liability. + +In this respect, the user's attention is drawn to the risks associated +with loading, using, modifying and/or developing or reproducing the +software by the user in light of its specific status of free software, +that may mean that it is complicated to manipulate, and that also +therefore means that it is reserved for developers and experienced +professionals having in-depth computer knowledge. Users are therefore +encouraged to load and test the software's suitability as regards their +requirements in conditions enabling the security of their systems and/or +data to be ensured and, more generally, to use and operate it in the +same conditions as regards security. + +The fact that you are presently reading this means that you have had +knowledge of the CeCILL license and that you accept its terms. ++*/ + #ifndef __COSMO_INTERPOLATE3D_HPP #define __COSMO_INTERPOLATE3D_HPP diff --git a/src/kdtree_leaf.hpp b/src/kdtree_leaf.hpp index 145fa57..0481889 100644 --- a/src/kdtree_leaf.hpp +++ b/src/kdtree_leaf.hpp @@ -1,3 +1,38 @@ +/*+ +This is CosmoTool (./src/kdtree_leaf.hpp) -- Copyright (C) Guilhem Lavaux (2007-2013) + +guilhem.lavaux@gmail.com + +This software is a computer program whose purpose is to provide a toolbox for cosmological +data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...) + +This software is governed by the CeCILL license under French law and +abiding by the rules of distribution of free software. You can use, +modify and/ or redistribute the software under the terms of the CeCILL +license as circulated by CEA, CNRS and INRIA at the following URL +"http://www.cecill.info". + +As a counterpart to the access to the source code and rights to copy, +modify and redistribute granted by the license, users are provided only +with a limited warranty and the software's author, the holder of the +economic rights, and the successive licensors have only limited +liability. + +In this respect, the user's attention is drawn to the risks associated +with loading, using, modifying and/or developing or reproducing the +software by the user in light of its specific status of free software, +that may mean that it is complicated to manipulate, and that also +therefore means that it is reserved for developers and experienced +professionals having in-depth computer knowledge. Users are therefore +encouraged to load and test the software's suitability as regards their +requirements in conditions enabling the security of their systems and/or +data to be ensured and, more generally, to use and operate it in the +same conditions as regards security. + +The fact that you are presently reading this means that you have had +knowledge of the CeCILL license and that you accept its terms. ++*/ + #ifndef __LEAF_KDTREE_HPP #define __LEAF_KDTREE_HPP diff --git a/src/kdtree_splitters.hpp b/src/kdtree_splitters.hpp index 4092cde..c0364a6 100644 --- a/src/kdtree_splitters.hpp +++ b/src/kdtree_splitters.hpp @@ -1,3 +1,38 @@ +/*+ +This is CosmoTool (./src/kdtree_splitters.hpp) -- Copyright (C) Guilhem Lavaux (2007-2013) + +guilhem.lavaux@gmail.com + +This software is a computer program whose purpose is to provide a toolbox for cosmological +data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...) + +This software is governed by the CeCILL license under French law and +abiding by the rules of distribution of free software. You can use, +modify and/ or redistribute the software under the terms of the CeCILL +license as circulated by CEA, CNRS and INRIA at the following URL +"http://www.cecill.info". + +As a counterpart to the access to the source code and rights to copy, +modify and redistribute granted by the license, users are provided only +with a limited warranty and the software's author, the holder of the +economic rights, and the successive licensors have only limited +liability. + +In this respect, the user's attention is drawn to the risks associated +with loading, using, modifying and/or developing or reproducing the +software by the user in light of its specific status of free software, +that may mean that it is complicated to manipulate, and that also +therefore means that it is reserved for developers and experienced +professionals having in-depth computer knowledge. Users are therefore +encouraged to load and test the software's suitability as regards their +requirements in conditions enabling the security of their systems and/or +data to be ensured and, more generally, to use and operate it in the +same conditions as regards security. + +The fact that you are presently reading this means that you have had +knowledge of the CeCILL license and that you accept its terms. ++*/ + #ifndef __KDTREE_SPLITTERS_HPP #define __KDTREE_SPLITTERS_HPP diff --git a/src/loadFlash.cpp b/src/loadFlash.cpp index ada498d..d437f25 100644 --- a/src/loadFlash.cpp +++ b/src/loadFlash.cpp @@ -1,3 +1,38 @@ +/*+ +This is CosmoTool (./src/loadFlash.cpp) -- Copyright (C) Guilhem Lavaux (2007-2013) + +guilhem.lavaux@gmail.com + +This software is a computer program whose purpose is to provide a toolbox for cosmological +data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...) + +This software is governed by the CeCILL license under French law and +abiding by the rules of distribution of free software. You can use, +modify and/ or redistribute the software under the terms of the CeCILL +license as circulated by CEA, CNRS and INRIA at the following URL +"http://www.cecill.info". + +As a counterpart to the access to the source code and rights to copy, +modify and redistribute granted by the license, users are provided only +with a limited warranty and the software's author, the holder of the +economic rights, and the successive licensors have only limited +liability. + +In this respect, the user's attention is drawn to the risks associated +with loading, using, modifying and/or developing or reproducing the +software by the user in light of its specific status of free software, +that may mean that it is complicated to manipulate, and that also +therefore means that it is reserved for developers and experienced +professionals having in-depth computer knowledge. Users are therefore +encouraged to load and test the software's suitability as regards their +requirements in conditions enabling the security of their systems and/or +data to be ensured and, more generally, to use and operate it in the +same conditions as regards security. + +The fact that you are presently reading this means that you have had +knowledge of the CeCILL license and that you accept its terms. ++*/ + /* Reads in FLASH v3 files in HDF5 format */ #include @@ -71,7 +106,7 @@ SimuData *CosmoTool::loadFlashMulti(const char *fname, int id, int loadflags) } } if (loadflags & NEED_GADGET_ID) { - data->Id = new int[data->NumPart]; + data->Id = new long[data->NumPart]; if (data->Id == 0) { delete data; return 0; diff --git a/src/loadFlash.hpp b/src/loadFlash.hpp index 3d7d8af..62f511b 100644 --- a/src/loadFlash.hpp +++ b/src/loadFlash.hpp @@ -1,3 +1,38 @@ +/*+ +This is CosmoTool (./src/loadFlash.hpp) -- Copyright (C) Guilhem Lavaux (2007-2013) + +guilhem.lavaux@gmail.com + +This software is a computer program whose purpose is to provide a toolbox for cosmological +data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...) + +This software is governed by the CeCILL license under French law and +abiding by the rules of distribution of free software. You can use, +modify and/ or redistribute the software under the terms of the CeCILL +license as circulated by CEA, CNRS and INRIA at the following URL +"http://www.cecill.info". + +As a counterpart to the access to the source code and rights to copy, +modify and redistribute granted by the license, users are provided only +with a limited warranty and the software's author, the holder of the +economic rights, and the successive licensors have only limited +liability. + +In this respect, the user's attention is drawn to the risks associated +with loading, using, modifying and/or developing or reproducing the +software by the user in light of its specific status of free software, +that may mean that it is complicated to manipulate, and that also +therefore means that it is reserved for developers and experienced +professionals having in-depth computer knowledge. Users are therefore +encouraged to load and test the software's suitability as regards their +requirements in conditions enabling the security of their systems and/or +data to be ensured and, more generally, to use and operate it in the +same conditions as regards security. + +The fact that you are presently reading this means that you have had +knowledge of the CeCILL license and that you accept its terms. ++*/ + #ifndef __COSMO_LOAD_FLASH_HPP #define __COSMO_LOAD_FLASH_HPP diff --git a/src/loadFlash_dummy.cpp b/src/loadFlash_dummy.cpp index 2df6604..2c6a86c 100644 --- a/src/loadFlash_dummy.cpp +++ b/src/loadFlash_dummy.cpp @@ -1,3 +1,38 @@ +/*+ +This is CosmoTool (./src/loadFlash_dummy.cpp) -- Copyright (C) Guilhem Lavaux (2007-2013) + +guilhem.lavaux@gmail.com + +This software is a computer program whose purpose is to provide a toolbox for cosmological +data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...) + +This software is governed by the CeCILL license under French law and +abiding by the rules of distribution of free software. You can use, +modify and/ or redistribute the software under the terms of the CeCILL +license as circulated by CEA, CNRS and INRIA at the following URL +"http://www.cecill.info". + +As a counterpart to the access to the source code and rights to copy, +modify and redistribute granted by the license, users are provided only +with a limited warranty and the software's author, the holder of the +economic rights, and the successive licensors have only limited +liability. + +In this respect, the user's attention is drawn to the risks associated +with loading, using, modifying and/or developing or reproducing the +software by the user in light of its specific status of free software, +that may mean that it is complicated to manipulate, and that also +therefore means that it is reserved for developers and experienced +professionals having in-depth computer knowledge. Users are therefore +encouraged to load and test the software's suitability as regards their +requirements in conditions enabling the security of their systems and/or +data to be ensured and, more generally, to use and operate it in the +same conditions as regards security. + +The fact that you are presently reading this means that you have had +knowledge of the CeCILL license and that you accept its terms. ++*/ + #include "load_data.hpp" #include "loadFlash.hpp" diff --git a/src/loadGadget.cpp b/src/loadGadget.cpp index 6e95bad..f9f1072 100644 --- a/src/loadGadget.cpp +++ b/src/loadGadget.cpp @@ -1,3 +1,38 @@ +/*+ +This is CosmoTool (./src/loadGadget.cpp) -- Copyright (C) Guilhem Lavaux (2007-2013) + +guilhem.lavaux@gmail.com + +This software is a computer program whose purpose is to provide a toolbox for cosmological +data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...) + +This software is governed by the CeCILL license under French law and +abiding by the rules of distribution of free software. You can use, +modify and/ or redistribute the software under the terms of the CeCILL +license as circulated by CEA, CNRS and INRIA at the following URL +"http://www.cecill.info". + +As a counterpart to the access to the source code and rights to copy, +modify and redistribute granted by the license, users are provided only +with a limited warranty and the software's author, the holder of the +economic rights, and the successive licensors have only limited +liability. + +In this respect, the user's attention is drawn to the risks associated +with loading, using, modifying and/or developing or reproducing the +software by the user in light of its specific status of free software, +that may mean that it is complicated to manipulate, and that also +therefore means that it is reserved for developers and experienced +professionals having in-depth computer knowledge. Users are therefore +encouraged to load and test the software's suitability as regards their +requirements in conditions enabling the security of their systems and/or +data to be ensured and, more generally, to use and operate it in the +same conditions as regards security. + +The fact that you are presently reading this means that you have had +knowledge of the CeCILL license and that you accept its terms. ++*/ + #include #include #include @@ -10,61 +45,41 @@ using namespace CosmoTool; using namespace std; -PurePositionData *CosmoTool::loadGadgetPosition(const char *fname) -{ - PurePositionData *data; - int p, n; - UnformattedRead f(fname); - GadgetHeader h; - - data = new PurePositionData; - f.beginCheckpoint(); - for (int i = 0; i < 6; i++) - h.npart[i] = f.readInt32(); - for (int i = 0; i < 6; i++) - h.mass[i] = f.readReal64(); - h.time = f.readReal64(); - h.redshift = f.readReal64(); - h.flag_sfr = f.readInt32(); - h.flag_feedback = f.readInt32(); - for (int i = 0; i < 6; i++) - h.npartTotal[i] = f.readInt32(); - h.flag_cooling = f.readInt32(); - h.num_files = f.readInt32(); - data->BoxSize = h.BoxSize = f.readReal64(); - h.Omega0 = f.readReal64(); - h.OmegaLambda = f.readReal64(); - h.HubbleParam = f.readReal64(); - f.endCheckpoint(true); - - data->NumPart = 0; - for(int k=0; k<5; k++) - data->NumPart += h.npart[k]; - - data->pos = new FCoordinates[data->NumPart]; - - f.beginCheckpoint(); - for(int k = 0, p = 0; k < 5; k++) { - for(int n = 0; n < h.npart[k]; n++) { - data->pos[p][0] = f.readReal32(); - data->pos[p][1] = f.readReal32(); - data->pos[p][2] = f.readReal32(); - p++; - } - } - f.endCheckpoint(); - - // Skip velocities - f.skip((long)data->NumPart*3+2*4); - // Skip ids - return data; +void loadGadgetHeader(UnformattedRead *f, GadgetHeader& h, SimuData *data, int id) +{ + f->beginCheckpoint(); + for (int i = 0; i < 6; i++) + h.npart[i] = f->readInt32(); + for (int i = 0; i < 6; i++) + h.mass[i] = f->readReal64(); + data->time = h.time = f->readReal64(); + h.redshift = f->readReal64(); + h.flag_sfr = f->readInt32(); + h.flag_feedback = f->readInt32(); + for (int i = 0; i < 6; i++) + h.npartTotal[i] = f->readInt32(); + h.flag_cooling = f->readInt32(); + h.num_files = f->readInt32(); + data->BoxSize = h.BoxSize = f->readReal64(); + data->Omega_M = h.Omega0 = f->readReal64(); + data->Omega_Lambda = h.OmegaLambda = f->readReal64(); + data->Hubble = h.HubbleParam = f->readReal64(); + f->endCheckpoint(true); + + long NumPart = 0, NumPartTotal = 0; + for(int k=0; k<6; k++) + { + NumPart += h.npart[k]; + NumPartTotal += (id < 0) ? h.npart[k] : h.npartTotal[k]; + } + data->NumPart = NumPart; + data->TotalNumPart = NumPartTotal; } - - - -SimuData *CosmoTool::loadGadgetMulti(const char *fname, int id, int loadflags, int GadgetFormat) +SimuData *CosmoTool::loadGadgetMulti(const char *fname, int id, + int loadflags, int GadgetFormat, + SimuFilter filter) { SimuData *data; int p, n; @@ -98,35 +113,11 @@ SimuData *CosmoTool::loadGadgetMulti(const char *fname, int id, int loadflags, i } long NumPart = 0, NumPartTotal = 0; - + try { - f->beginCheckpoint(); - for (int i = 0; i < 6; i++) - h.npart[i] = f->readInt32(); - for (int i = 0; i < 6; i++) - h.mass[i] = f->readReal64(); - data->time = h.time = f->readReal64(); - h.redshift = f->readReal64(); - h.flag_sfr = f->readInt32(); - h.flag_feedback = f->readInt32(); - for (int i = 0; i < 6; i++) - h.npartTotal[i] = f->readInt32(); - h.flag_cooling = f->readInt32(); - h.num_files = f->readInt32(); - data->BoxSize = h.BoxSize = f->readReal64(); - data->Omega_M = h.Omega0 = f->readReal64(); - data->Omega_Lambda = h.OmegaLambda = f->readReal64(); - data->Hubble = h.HubbleParam = f->readReal64(); - f->endCheckpoint(true); - - for(int k=0; k<6; k++) - { - NumPart += h.npart[k]; - NumPartTotal += (id < 0) ? h.npart[k] : h.npartTotal[k]; - } - data->NumPart = NumPart; - data->TotalNumPart = NumPartTotal; + loadGadgetHeader(f, h, data, id); + if (GadgetFormat == 1) velmul = sqrt(h.time); else if (GadgetFormat == 2) @@ -135,6 +126,9 @@ SimuData *CosmoTool::loadGadgetMulti(const char *fname, int id, int loadflags, i cerr << "unknown gadget format" << endl; abort(); } + + NumPart = data->NumPart; + NumPartTotal = data->TotalNumPart; } catch (const InvalidUnformattedAccess& e) { @@ -236,7 +230,7 @@ SimuData *CosmoTool::loadGadgetMulti(const char *fname, int id, int loadflags, i try { f->beginCheckpoint(); - data->Id = new int[data->NumPart]; + data->Id = new long[data->NumPart]; if (data->Id == 0) { delete f; diff --git a/src/loadGadget.hpp b/src/loadGadget.hpp index eec0fc9..e05a187 100644 --- a/src/loadGadget.hpp +++ b/src/loadGadget.hpp @@ -1,3 +1,38 @@ +/*+ +This is CosmoTool (./src/loadGadget.hpp) -- Copyright (C) Guilhem Lavaux (2007-2013) + +guilhem.lavaux@gmail.com + +This software is a computer program whose purpose is to provide a toolbox for cosmological +data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...) + +This software is governed by the CeCILL license under French law and +abiding by the rules of distribution of free software. You can use, +modify and/ or redistribute the software under the terms of the CeCILL +license as circulated by CEA, CNRS and INRIA at the following URL +"http://www.cecill.info". + +As a counterpart to the access to the source code and rights to copy, +modify and redistribute granted by the license, users are provided only +with a limited warranty and the software's author, the holder of the +economic rights, and the successive licensors have only limited +liability. + +In this respect, the user's attention is drawn to the risks associated +with loading, using, modifying and/or developing or reproducing the +software by the user in light of its specific status of free software, +that may mean that it is complicated to manipulate, and that also +therefore means that it is reserved for developers and experienced +professionals having in-depth computer knowledge. Users are therefore +encouraged to load and test the software's suitability as regards their +requirements in conditions enabling the security of their systems and/or +data to be ensured and, more generally, to use and operate it in the +same conditions as regards security. + +The fact that you are presently reading this means that you have had +knowledge of the CeCILL license and that you accept its terms. ++*/ + #ifndef __COSMO_LOAD_GADGET_HPP #define __COSMO_LOAD_GADGET_HPP @@ -7,10 +42,9 @@ #include "loadSimu.hpp" namespace CosmoTool { - - PurePositionData *loadGadgetPosition(const char *fname); - SimuData *loadGadgetMulti(const char *fname, int id, int flags, int GadgetFormat = 1); + SimuData *loadGadgetMulti(const char *fname, int id, int flags, + int GadgetFormat = 1, SimuFilter filter = 0); // Only single snapshot supported void writeGadget(const std::string& fname, SimuData *data, int GadgetFormat = 1); diff --git a/src/loadRamses.cpp b/src/loadRamses.cpp index 4dcbad2..3eb909c 100644 --- a/src/loadRamses.cpp +++ b/src/loadRamses.cpp @@ -1,3 +1,38 @@ +/*+ +This is CosmoTool (./src/loadRamses.cpp) -- Copyright (C) Guilhem Lavaux (2007-2013) + +guilhem.lavaux@gmail.com + +This software is a computer program whose purpose is to provide a toolbox for cosmological +data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...) + +This software is governed by the CeCILL license under French law and +abiding by the rules of distribution of free software. You can use, +modify and/ or redistribute the software under the terms of the CeCILL +license as circulated by CEA, CNRS and INRIA at the following URL +"http://www.cecill.info". + +As a counterpart to the access to the source code and rights to copy, +modify and redistribute granted by the license, users are provided only +with a limited warranty and the software's author, the holder of the +economic rights, and the successive licensors have only limited +liability. + +In this respect, the user's attention is drawn to the risks associated +with loading, using, modifying and/or developing or reproducing the +software by the user in light of its specific status of free software, +that may mean that it is complicated to manipulate, and that also +therefore means that it is reserved for developers and experienced +professionals having in-depth computer knowledge. Users are therefore +encouraged to load and test the software's suitability as regards their +requirements in conditions enabling the security of their systems and/or +data to be ensured and, more generally, to use and operate it in the +same conditions as regards security. + +The fact that you are presently reading this means that you have had +knowledge of the CeCILL license and that you accept its terms. ++*/ + #include #include #include @@ -318,6 +353,9 @@ CosmoTool::SimuData *CosmoTool::loadRamsesSimu(const char *basename, int outputI data->Omega_M = info.omega_m; data->Omega_Lambda = info.omega_lambda; + if (flags == 0) + return data; + if (cpuid < 0) cpuid = 1; else cpuid++; @@ -381,7 +419,7 @@ CosmoTool::SimuData *CosmoTool::loadRamsesSimu(const char *basename, int outputI } if (flags & NEED_GADGET_ID) { - data->Id = new int[nPar]; + data->Id = new long[nPar]; } for (int k = 0; k < 3; k++) diff --git a/src/loadRamses.hpp b/src/loadRamses.hpp index 15e6275..5282a3c 100644 --- a/src/loadRamses.hpp +++ b/src/loadRamses.hpp @@ -1,3 +1,38 @@ +/*+ +This is CosmoTool (./src/loadRamses.hpp) -- Copyright (C) Guilhem Lavaux (2007-2013) + +guilhem.lavaux@gmail.com + +This software is a computer program whose purpose is to provide a toolbox for cosmological +data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...) + +This software is governed by the CeCILL license under French law and +abiding by the rules of distribution of free software. You can use, +modify and/ or redistribute the software under the terms of the CeCILL +license as circulated by CEA, CNRS and INRIA at the following URL +"http://www.cecill.info". + +As a counterpart to the access to the source code and rights to copy, +modify and redistribute granted by the license, users are provided only +with a limited warranty and the software's author, the holder of the +economic rights, and the successive licensors have only limited +liability. + +In this respect, the user's attention is drawn to the risks associated +with loading, using, modifying and/or developing or reproducing the +software by the user in light of its specific status of free software, +that may mean that it is complicated to manipulate, and that also +therefore means that it is reserved for developers and experienced +professionals having in-depth computer knowledge. Users are therefore +encouraged to load and test the software's suitability as regards their +requirements in conditions enabling the security of their systems and/or +data to be ensured and, more generally, to use and operate it in the +same conditions as regards security. + +The fact that you are presently reading this means that you have had +knowledge of the CeCILL license and that you accept its terms. ++*/ + #ifndef _LOAD_RAMSES_HPP #define _LOAD_RAMSES_HPP diff --git a/src/loadSimu.hpp b/src/loadSimu.hpp index a7db941..1c360a4 100644 --- a/src/loadSimu.hpp +++ b/src/loadSimu.hpp @@ -1,6 +1,43 @@ +/*+ +This is CosmoTool (./src/loadSimu.hpp) -- Copyright (C) Guilhem Lavaux (2007-2013) + +guilhem.lavaux@gmail.com + +This software is a computer program whose purpose is to provide a toolbox for cosmological +data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...) + +This software is governed by the CeCILL license under French law and +abiding by the rules of distribution of free software. You can use, +modify and/ or redistribute the software under the terms of the CeCILL +license as circulated by CEA, CNRS and INRIA at the following URL +"http://www.cecill.info". + +As a counterpart to the access to the source code and rights to copy, +modify and redistribute granted by the license, users are provided only +with a limited warranty and the software's author, the holder of the +economic rights, and the successive licensors have only limited +liability. + +In this respect, the user's attention is drawn to the risks associated +with loading, using, modifying and/or developing or reproducing the +software by the user in light of its specific status of free software, +that may mean that it is complicated to manipulate, and that also +therefore means that it is reserved for developers and experienced +professionals having in-depth computer knowledge. Users are therefore +encouraged to load and test the software's suitability as regards their +requirements in conditions enabling the security of their systems and/or +data to be ensured and, more generally, to use and operate it in the +same conditions as regards security. + +The fact that you are presently reading this means that you have had +knowledge of the CeCILL license and that you accept its terms. ++*/ + #ifndef __COSMOTOOLBOX_HPP #define __COSMOTOOLBOX_HPP +#include +#include namespace CosmoTool { @@ -9,9 +46,24 @@ namespace CosmoTool static const int NEED_VELOCITY = 4; static const int NEED_TYPE = 8; + struct SimuParticle + { + float Pos[3]; + float Vel[3]; + int type; + int id; + + bool flag_vel, flag_type, flag_id; + }; + + typedef bool (*SimuFilter)(const SimuParticle& p); + class SimuData { public: + typedef void (*FreeFunction)(void *); + typedef std::map > AttributeMap; + float BoxSize; float time; float Hubble; @@ -21,10 +73,13 @@ namespace CosmoTool long NumPart; long TotalNumPart; - int *Id; + long *Id; float *Pos[3]; float *Vel[3]; int *type; + + AttributeMap attributes; + public: SimuData() : Id(0),NumPart(0),type(0) { Pos[0]=Pos[1]=Pos[2]=0; Vel[0]=Vel[1]=Vel[2]=0; } ~SimuData() @@ -40,7 +95,37 @@ namespace CosmoTool delete[] type; if (Id) delete[] Id; + + for (AttributeMap::iterator i = attributes.begin(); + i != attributes.end(); + ++i) + { + if (i->second.second) + i->second.second(i->second.first); + } } + + template + T *as(const std::string& n) + { + AttributeMap::iterator i = attributes.find(n); + if (i == attributes.end()) + return 0; + + return reinterpret_cast(i->second.first); + } + + void new_attribute(const std::string& n, void *p, FreeFunction free_func) + { + AttributeMap::iterator i = attributes.find(n); + if (i != attributes.end()) + { + if (i->second.second) + i->second.second(i->second.first); + } + attributes[n] = std::make_pair(p, free_func); + } + }; }; diff --git a/src/load_data.cpp b/src/load_data.cpp index 28c7dc6..f23f909 100644 --- a/src/load_data.cpp +++ b/src/load_data.cpp @@ -1,3 +1,38 @@ +/*+ +This is CosmoTool (./src/load_data.cpp) -- Copyright (C) Guilhem Lavaux (2007-2013) + +guilhem.lavaux@gmail.com + +This software is a computer program whose purpose is to provide a toolbox for cosmological +data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...) + +This software is governed by the CeCILL license under French law and +abiding by the rules of distribution of free software. You can use, +modify and/ or redistribute the software under the terms of the CeCILL +license as circulated by CEA, CNRS and INRIA at the following URL +"http://www.cecill.info". + +As a counterpart to the access to the source code and rights to copy, +modify and redistribute granted by the license, users are provided only +with a limited warranty and the software's author, the holder of the +economic rights, and the successive licensors have only limited +liability. + +In this respect, the user's attention is drawn to the risks associated +with loading, using, modifying and/or developing or reproducing the +software by the user in light of its specific status of free software, +that may mean that it is complicated to manipulate, and that also +therefore means that it is reserved for developers and experienced +professionals having in-depth computer knowledge. Users are therefore +encouraged to load and test the software's suitability as regards their +requirements in conditions enabling the security of their systems and/or +data to be ensured and, more generally, to use and operate it in the +same conditions as regards security. + +The fact that you are presently reading this means that you have had +knowledge of the CeCILL license and that you accept its terms. ++*/ + #include #include #include diff --git a/src/load_data.hpp b/src/load_data.hpp index c4bc611..6be3bbb 100644 --- a/src/load_data.hpp +++ b/src/load_data.hpp @@ -1,3 +1,38 @@ +/*+ +This is CosmoTool (./src/load_data.hpp) -- Copyright (C) Guilhem Lavaux (2007-2013) + +guilhem.lavaux@gmail.com + +This software is a computer program whose purpose is to provide a toolbox for cosmological +data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...) + +This software is governed by the CeCILL license under French law and +abiding by the rules of distribution of free software. You can use, +modify and/ or redistribute the software under the terms of the CeCILL +license as circulated by CEA, CNRS and INRIA at the following URL +"http://www.cecill.info". + +As a counterpart to the access to the source code and rights to copy, +modify and redistribute granted by the license, users are provided only +with a limited warranty and the software's author, the holder of the +economic rights, and the successive licensors have only limited +liability. + +In this respect, the user's attention is drawn to the risks associated +with loading, using, modifying and/or developing or reproducing the +software by the user in light of its specific status of free software, +that may mean that it is complicated to manipulate, and that also +therefore means that it is reserved for developers and experienced +professionals having in-depth computer knowledge. Users are therefore +encouraged to load and test the software's suitability as regards their +requirements in conditions enabling the security of their systems and/or +data to be ensured and, more generally, to use and operate it in the +same conditions as regards security. + +The fact that you are presently reading this means that you have had +knowledge of the CeCILL license and that you accept its terms. ++*/ + #ifndef _LOAD_GADGET_DATA_HPP #define _LOAD_GADGET_DATA_HPP diff --git a/src/mach.hpp b/src/mach.hpp index abcb22a..6bea525 100644 --- a/src/mach.hpp +++ b/src/mach.hpp @@ -1,3 +1,38 @@ +/*+ +This is CosmoTool (./src/mach.hpp) -- Copyright (C) Guilhem Lavaux (2007-2013) + +guilhem.lavaux@gmail.com + +This software is a computer program whose purpose is to provide a toolbox for cosmological +data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...) + +This software is governed by the CeCILL license under French law and +abiding by the rules of distribution of free software. You can use, +modify and/ or redistribute the software under the terms of the CeCILL +license as circulated by CEA, CNRS and INRIA at the following URL +"http://www.cecill.info". + +As a counterpart to the access to the source code and rights to copy, +modify and redistribute granted by the license, users are provided only +with a limited warranty and the software's author, the holder of the +economic rights, and the successive licensors have only limited +liability. + +In this respect, the user's attention is drawn to the risks associated +with loading, using, modifying and/or developing or reproducing the +software by the user in light of its specific status of free software, +that may mean that it is complicated to manipulate, and that also +therefore means that it is reserved for developers and experienced +professionals having in-depth computer knowledge. Users are therefore +encouraged to load and test the software's suitability as regards their +requirements in conditions enabling the security of their systems and/or +data to be ensured and, more generally, to use and operate it in the +same conditions as regards security. + +The fact that you are presently reading this means that you have had +knowledge of the CeCILL license and that you accept its terms. ++*/ + #ifndef __COSMO_MACHINE_TEST_HPP #define __COSMO_MACHINE_TEST_HPP diff --git a/src/miniargs.cpp b/src/miniargs.cpp index 4ad71c6..7c23230 100644 --- a/src/miniargs.cpp +++ b/src/miniargs.cpp @@ -1,3 +1,38 @@ +/*+ +This is CosmoTool (./src/miniargs.cpp) -- Copyright (C) Guilhem Lavaux (2007-2013) + +guilhem.lavaux@gmail.com + +This software is a computer program whose purpose is to provide a toolbox for cosmological +data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...) + +This software is governed by the CeCILL license under French law and +abiding by the rules of distribution of free software. You can use, +modify and/ or redistribute the software under the terms of the CeCILL +license as circulated by CEA, CNRS and INRIA at the following URL +"http://www.cecill.info". + +As a counterpart to the access to the source code and rights to copy, +modify and redistribute granted by the license, users are provided only +with a limited warranty and the software's author, the holder of the +economic rights, and the successive licensors have only limited +liability. + +In this respect, the user's attention is drawn to the risks associated +with loading, using, modifying and/or developing or reproducing the +software by the user in light of its specific status of free software, +that may mean that it is complicated to manipulate, and that also +therefore means that it is reserved for developers and experienced +professionals having in-depth computer knowledge. Users are therefore +encouraged to load and test the software's suitability as regards their +requirements in conditions enabling the security of their systems and/or +data to be ensured and, more generally, to use and operate it in the +same conditions as regards security. + +The fact that you are presently reading this means that you have had +knowledge of the CeCILL license and that you accept its terms. ++*/ + #include #include #include diff --git a/src/miniargs.hpp b/src/miniargs.hpp index b746e55..9f04b58 100644 --- a/src/miniargs.hpp +++ b/src/miniargs.hpp @@ -1,3 +1,38 @@ +/*+ +This is CosmoTool (./src/miniargs.hpp) -- Copyright (C) Guilhem Lavaux (2007-2013) + +guilhem.lavaux@gmail.com + +This software is a computer program whose purpose is to provide a toolbox for cosmological +data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...) + +This software is governed by the CeCILL license under French law and +abiding by the rules of distribution of free software. You can use, +modify and/ or redistribute the software under the terms of the CeCILL +license as circulated by CEA, CNRS and INRIA at the following URL +"http://www.cecill.info". + +As a counterpart to the access to the source code and rights to copy, +modify and redistribute granted by the license, users are provided only +with a limited warranty and the software's author, the holder of the +economic rights, and the successive licensors have only limited +liability. + +In this respect, the user's attention is drawn to the risks associated +with loading, using, modifying and/or developing or reproducing the +software by the user in light of its specific status of free software, +that may mean that it is complicated to manipulate, and that also +therefore means that it is reserved for developers and experienced +professionals having in-depth computer knowledge. Users are therefore +encouraged to load and test the software's suitability as regards their +requirements in conditions enabling the security of their systems and/or +data to be ensured and, more generally, to use and operate it in the +same conditions as regards security. + +The fact that you are presently reading this means that you have had +knowledge of the CeCILL license and that you accept its terms. ++*/ + #ifndef _MAK_MINIARGS_HPP #define _MAK_MINIARGS_HPP diff --git a/src/mykdtree.hpp b/src/mykdtree.hpp index c2087b5..27d66ae 100644 --- a/src/mykdtree.hpp +++ b/src/mykdtree.hpp @@ -1,3 +1,38 @@ +/*+ +This is CosmoTool (./src/mykdtree.hpp) -- Copyright (C) Guilhem Lavaux (2007-2013) + +guilhem.lavaux@gmail.com + +This software is a computer program whose purpose is to provide a toolbox for cosmological +data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...) + +This software is governed by the CeCILL license under French law and +abiding by the rules of distribution of free software. You can use, +modify and/ or redistribute the software under the terms of the CeCILL +license as circulated by CEA, CNRS and INRIA at the following URL +"http://www.cecill.info". + +As a counterpart to the access to the source code and rights to copy, +modify and redistribute granted by the license, users are provided only +with a limited warranty and the software's author, the holder of the +economic rights, and the successive licensors have only limited +liability. + +In this respect, the user's attention is drawn to the risks associated +with loading, using, modifying and/or developing or reproducing the +software by the user in light of its specific status of free software, +that may mean that it is complicated to manipulate, and that also +therefore means that it is reserved for developers and experienced +professionals having in-depth computer knowledge. Users are therefore +encouraged to load and test the software's suitability as regards their +requirements in conditions enabling the security of their systems and/or +data to be ensured and, more generally, to use and operate it in the +same conditions as regards security. + +The fact that you are presently reading this means that you have had +knowledge of the CeCILL license and that you accept its terms. ++*/ + #ifndef __HV_KDTREE_HPP #define __HV_KDTREE_HPP @@ -65,15 +100,16 @@ namespace CosmoTool { class RecursionMultipleInfo { public: - const typename KDDef::KDCoordinates& x; + typename KDDef::KDCoordinates x; BoundedQueue< KDCell *, typename KDDef::CoordType> queue; int traversed; RecursionMultipleInfo(const typename KDDef::KDCoordinates& rx, KDCell **cells, uint32_t numCells) - : x(rx), queue(cells, numCells, INFINITY),traversed(0) - { + : queue(cells, numCells, INFINITY),traversed(0) + { + std::copy(rx, rx+N, x); } }; @@ -97,6 +133,12 @@ namespace CosmoTool { KDTree(Cell *cells, uint32_t Ncells); ~KDTree(); + void setPeriodic(bool on, CoordType replicate) + { + periodic = on; + this->replicate = replicate; + } + uint32_t getIntersection(const coords& x, CoordType r, Cell **cells, uint32_t numCells) @@ -150,6 +192,9 @@ namespace CosmoTool { Cell **sortingHelper; Cell *base_cell; + bool periodic; + CoordType replicate; + Node *buildTree(Cell **cell0, uint32_t NumCells, uint32_t depth, diff --git a/src/mykdtree.tcc b/src/mykdtree.tcc index a5e3abb..e276cad 100644 --- a/src/mykdtree.tcc +++ b/src/mykdtree.tcc @@ -1,3 +1,4 @@ +#include "replicateGenerator.hpp" #include #include #include @@ -31,7 +32,8 @@ namespace CosmoTool { template KDTree::KDTree(Cell *cells, uint32_t Ncells) { - + periodic = false; + base_cell = cells; numNodes = Ncells; nodes = new Node[numNodes]; @@ -93,6 +95,19 @@ namespace CosmoTool { info.distances = 0; recursiveIntersectionCells(info, root, 0); + if (periodic) + { + ReplicateGenerator r(x, replicate); + + do + { + coords x_new; + r.getPosition(info.x); + recursiveIntersectionCells(info, root, 0); + } + while (r.next()); + } + return info.currentRank; } @@ -114,6 +129,18 @@ namespace CosmoTool { info.distances = distances; recursiveIntersectionCells(info, root, 0); + if (periodic) + { + ReplicateGenerator r(x, replicate); + + do + { + coords x_new; + r.getPosition(info.x); + recursiveIntersectionCells(info, root, 0); + } + while (r.next()); + } return info.currentRank; } @@ -131,6 +158,19 @@ namespace CosmoTool { info.distances = 0; recursiveIntersectionCells(info, root, 0); + if (periodic) + { + ReplicateGenerator r(x, replicate); + + do + { + coords x_new; + r.getPosition(info.x); + recursiveIntersectionCells(info, root, 0); + } + while (r.next()); + } + return info.currentRank; } @@ -365,7 +405,19 @@ namespace CosmoTool { CoordType R2 = INFINITY; Cell *best = 0; - recursiveNearest(root, 0, x, R2, best); + recursiveNearest(root, 0, x, R2, best); + if (periodic) + { + ReplicateGenerator r(x, replicate); + + do + { + coords x_new; + r.getPosition(x_new); + recursiveNearest(root, 0, x_new, R2, best); + } + while (r.next()); + } return best; } @@ -437,6 +489,18 @@ namespace CosmoTool { cells[i] = 0; recursiveMultipleNearest(info, root, 0); + if (periodic) + { + ReplicateGenerator r(x, replicate); + + do + { + coords x_new; + r.getPosition(info.x); + recursiveMultipleNearest(info, root, 0); + } + while (r.next()); + } // std::cout << "Traversed = " << info.traversed << std::endl; } @@ -452,9 +516,19 @@ namespace CosmoTool { cells[i] = 0; recursiveMultipleNearest(info, root, 0); - memcpy(distances, info.queue.getPriorities(), sizeof(CoordType)*N2); + if (periodic) + { + ReplicateGenerator r(x, replicate); - // std::cout << "Traversed = " << info.traversed << std::endl; + do + { + coords x_new; + r.getPosition(info.x); + recursiveMultipleNearest(info, root, 0); + } + while (r.next()); + } + memcpy(distances, info.queue.getPriorities(), sizeof(CoordType)*N2); } #ifdef __KD_TREE_SAVE_ON_DISK diff --git a/src/newton.hpp b/src/newton.hpp index 3e861f0..d3975fe 100644 --- a/src/newton.hpp +++ b/src/newton.hpp @@ -1,3 +1,38 @@ +/*+ +This is CosmoTool (./src/newton.hpp) -- Copyright (C) Guilhem Lavaux (2007-2013) + +guilhem.lavaux@gmail.com + +This software is a computer program whose purpose is to provide a toolbox for cosmological +data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...) + +This software is governed by the CeCILL license under French law and +abiding by the rules of distribution of free software. You can use, +modify and/ or redistribute the software under the terms of the CeCILL +license as circulated by CEA, CNRS and INRIA at the following URL +"http://www.cecill.info". + +As a counterpart to the access to the source code and rights to copy, +modify and redistribute granted by the license, users are provided only +with a limited warranty and the software's author, the holder of the +economic rights, and the successive licensors have only limited +liability. + +In this respect, the user's attention is drawn to the risks associated +with loading, using, modifying and/or developing or reproducing the +software by the user in light of its specific status of free software, +that may mean that it is complicated to manipulate, and that also +therefore means that it is reserved for developers and experienced +professionals having in-depth computer knowledge. Users are therefore +encouraged to load and test the software's suitability as regards their +requirements in conditions enabling the security of their systems and/or +data to be ensured and, more generally, to use and operate it in the +same conditions as regards security. + +The fact that you are presently reading this means that you have had +knowledge of the CeCILL license and that you accept its terms. ++*/ + #ifndef _COSMOTOOL_NEWTON_HPP #define _COSMOTOOL_NEWTON_HPP diff --git a/src/octTree.cpp b/src/octTree.cpp index 4ef929d..2c76609 100644 --- a/src/octTree.cpp +++ b/src/octTree.cpp @@ -1,3 +1,38 @@ +/*+ +This is CosmoTool (./src/octTree.cpp) -- Copyright (C) Guilhem Lavaux (2007-2013) + +guilhem.lavaux@gmail.com + +This software is a computer program whose purpose is to provide a toolbox for cosmological +data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...) + +This software is governed by the CeCILL license under French law and +abiding by the rules of distribution of free software. You can use, +modify and/ or redistribute the software under the terms of the CeCILL +license as circulated by CEA, CNRS and INRIA at the following URL +"http://www.cecill.info". + +As a counterpart to the access to the source code and rights to copy, +modify and redistribute granted by the license, users are provided only +with a limited warranty and the software's author, the holder of the +economic rights, and the successive licensors have only limited +liability. + +In this respect, the user's attention is drawn to the risks associated +with loading, using, modifying and/or developing or reproducing the +software by the user in light of its specific status of free software, +that may mean that it is complicated to manipulate, and that also +therefore means that it is reserved for developers and experienced +professionals having in-depth computer knowledge. Users are therefore +encouraged to load and test the software's suitability as regards their +requirements in conditions enabling the security of their systems and/or +data to be ensured and, more generally, to use and operate it in the +same conditions as regards security. + +The fact that you are presently reading this means that you have had +knowledge of the CeCILL license and that you accept its terms. ++*/ + #include #include #include diff --git a/src/octTree.hpp b/src/octTree.hpp index 447a6b0..425d8fd 100644 --- a/src/octTree.hpp +++ b/src/octTree.hpp @@ -1,3 +1,38 @@ +/*+ +This is CosmoTool (./src/octTree.hpp) -- Copyright (C) Guilhem Lavaux (2007-2013) + +guilhem.lavaux@gmail.com + +This software is a computer program whose purpose is to provide a toolbox for cosmological +data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...) + +This software is governed by the CeCILL license under French law and +abiding by the rules of distribution of free software. You can use, +modify and/ or redistribute the software under the terms of the CeCILL +license as circulated by CEA, CNRS and INRIA at the following URL +"http://www.cecill.info". + +As a counterpart to the access to the source code and rights to copy, +modify and redistribute granted by the license, users are provided only +with a limited warranty and the software's author, the holder of the +economic rights, and the successive licensors have only limited +liability. + +In this respect, the user's attention is drawn to the risks associated +with loading, using, modifying and/or developing or reproducing the +software by the user in light of its specific status of free software, +that may mean that it is complicated to manipulate, and that also +therefore means that it is reserved for developers and experienced +professionals having in-depth computer knowledge. Users are therefore +encouraged to load and test the software's suitability as regards their +requirements in conditions enabling the security of their systems and/or +data to be ensured and, more generally, to use and operate it in the +same conditions as regards security. + +The fact that you are presently reading this means that you have had +knowledge of the CeCILL license and that you accept its terms. ++*/ + #ifndef __COSMOTOOL_AMR_HPP #define __COSMOTOOL_AMR_HPP diff --git a/src/pool.hpp b/src/pool.hpp index 5d2203b..b3e2637 100644 --- a/src/pool.hpp +++ b/src/pool.hpp @@ -1,3 +1,38 @@ +/*+ +This is CosmoTool (./src/pool.hpp) -- Copyright (C) Guilhem Lavaux (2007-2013) + +guilhem.lavaux@gmail.com + +This software is a computer program whose purpose is to provide a toolbox for cosmological +data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...) + +This software is governed by the CeCILL license under French law and +abiding by the rules of distribution of free software. You can use, +modify and/ or redistribute the software under the terms of the CeCILL +license as circulated by CEA, CNRS and INRIA at the following URL +"http://www.cecill.info". + +As a counterpart to the access to the source code and rights to copy, +modify and redistribute granted by the license, users are provided only +with a limited warranty and the software's author, the holder of the +economic rights, and the successive licensors have only limited +liability. + +In this respect, the user's attention is drawn to the risks associated +with loading, using, modifying and/or developing or reproducing the +software by the user in light of its specific status of free software, +that may mean that it is complicated to manipulate, and that also +therefore means that it is reserved for developers and experienced +professionals having in-depth computer knowledge. Users are therefore +encouraged to load and test the software's suitability as regards their +requirements in conditions enabling the security of their systems and/or +data to be ensured and, more generally, to use and operate it in the +same conditions as regards security. + +The fact that you are presently reading this means that you have had +knowledge of the CeCILL license and that you accept its terms. ++*/ + #ifndef __COSMO_POOL_HPP #define __COSMO_POOL_HPP diff --git a/src/powerSpectrum.cpp b/src/powerSpectrum.cpp index 9f3d7c8..b2a34be 100644 --- a/src/powerSpectrum.cpp +++ b/src/powerSpectrum.cpp @@ -1,3 +1,39 @@ +/*+ +This is CosmoTool (./src/powerSpectrum.cpp) -- Copyright (C) Guilhem Lavaux (2007-2013) + +guilhem.lavaux@gmail.com + +This software is a computer program whose purpose is to provide a toolbox for cosmological +data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...) + +This software is governed by the CeCILL license under French law and +abiding by the rules of distribution of free software. You can use, +modify and/ or redistribute the software under the terms of the CeCILL +license as circulated by CEA, CNRS and INRIA at the following URL +"http://www.cecill.info". + +As a counterpart to the access to the source code and rights to copy, +modify and redistribute granted by the license, users are provided only +with a limited warranty and the software's author, the holder of the +economic rights, and the successive licensors have only limited +liability. + +In this respect, the user's attention is drawn to the risks associated +with loading, using, modifying and/or developing or reproducing the +software by the user in light of its specific status of free software, +that may mean that it is complicated to manipulate, and that also +therefore means that it is reserved for developers and experienced +professionals having in-depth computer knowledge. Users are therefore +encouraged to load and test the software's suitability as regards their +requirements in conditions enabling the security of their systems and/or +data to be ensured and, more generally, to use and operate it in the +same conditions as regards security. + +The fact that you are presently reading this means that you have had +knowledge of the CeCILL license and that you accept its terms. ++*/ + +#include #include #include #include @@ -21,7 +57,7 @@ using namespace std; #define POWER_BDM 7 #define POWER_TEST 8 -#define POWER_SPECTRUM HU_WIGGLES +#define POWER_SPECTRUM POWER_EFSTATHIOU namespace Cosmology { @@ -651,4 +687,17 @@ double computeCorrel2(double powNorm, double topHatRad1, double topHatRad2) #endif } + double vvCorrection(double P_deltadelta, double k) + { + static const double alpha0 = -12480.5, alpha1 = 1.824, alpha2 = 2165.87, alpha3=1.796; + if (k > 0.3) + return 0; + double r =(alpha0*sqrt(P_deltadelta) + alpha1*P_deltadelta*P_deltadelta)/(alpha2 + alpha3*P_deltadelta); + assert(P_deltadelta > 0); + + if (r < 0) + return 0; + return r; + } + }; diff --git a/src/powerSpectrum.hpp b/src/powerSpectrum.hpp index ad0b0a9..d4c6761 100644 --- a/src/powerSpectrum.hpp +++ b/src/powerSpectrum.hpp @@ -1,3 +1,38 @@ +/*+ +This is CosmoTool (./src/powerSpectrum.hpp) -- Copyright (C) Guilhem Lavaux (2007-2013) + +guilhem.lavaux@gmail.com + +This software is a computer program whose purpose is to provide a toolbox for cosmological +data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...) + +This software is governed by the CeCILL license under French law and +abiding by the rules of distribution of free software. You can use, +modify and/ or redistribute the software under the terms of the CeCILL +license as circulated by CEA, CNRS and INRIA at the following URL +"http://www.cecill.info". + +As a counterpart to the access to the source code and rights to copy, +modify and redistribute granted by the license, users are provided only +with a limited warranty and the software's author, the holder of the +economic rights, and the successive licensors have only limited +liability. + +In this respect, the user's attention is drawn to the risks associated +with loading, using, modifying and/or developing or reproducing the +software by the user in light of its specific status of free software, +that may mean that it is complicated to manipulate, and that also +therefore means that it is reserved for developers and experienced +professionals having in-depth computer knowledge. Users are therefore +encouraged to load and test the software's suitability as regards their +requirements in conditions enabling the security of their systems and/or +data to be ensured and, more generally, to use and operate it in the +same conditions as regards security. + +The fact that you are presently reading this means that you have had +knowledge of the CeCILL license and that you accept its terms. ++*/ + #ifndef _POWERSPECTRUM_HPP #define _POWERSPECTRUM_HPP @@ -30,6 +65,8 @@ namespace Cosmology { double computeVarianceZero(double powNorm); double computeCorrel(double powNorm, double topHatRad1); double computeCorrel2(double powNorm, double topHatRad1, double topHatRad2); + + double vvCorrection(double P_deltadelta, double k); }; #endif diff --git a/src/replicateGenerator.hpp b/src/replicateGenerator.hpp new file mode 100644 index 0000000..fd7035f --- /dev/null +++ b/src/replicateGenerator.hpp @@ -0,0 +1,100 @@ +/*+ +This is CosmoTool (./src/replicateGenerator.hpp) -- Copyright (C) Guilhem Lavaux (2007-2013) + +guilhem.lavaux@gmail.com + +This software is a computer program whose purpose is to provide a toolbox for cosmological +data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...) + +This software is governed by the CeCILL license under French law and +abiding by the rules of distribution of free software. You can use, +modify and/ or redistribute the software under the terms of the CeCILL +license as circulated by CEA, CNRS and INRIA at the following URL +"http://www.cecill.info". + +As a counterpart to the access to the source code and rights to copy, +modify and redistribute granted by the license, users are provided only +with a limited warranty and the software's author, the holder of the +economic rights, and the successive licensors have only limited +liability. + +In this respect, the user's attention is drawn to the risks associated +with loading, using, modifying and/or developing or reproducing the +software by the user in light of its specific status of free software, +that may mean that it is complicated to manipulate, and that also +therefore means that it is reserved for developers and experienced +professionals having in-depth computer knowledge. Users are therefore +encouraged to load and test the software's suitability as regards their +requirements in conditions enabling the security of their systems and/or +data to be ensured and, more generally, to use and operate it in the +same conditions as regards security. + +The fact that you are presently reading this means that you have had +knowledge of the CeCILL license and that you accept its terms. ++*/ +#ifndef __REPLICATE_GENERATOR_HPP +#define __REPLICATE_GENERATOR_HPP + +#include +#include "algo.hpp" + +namespace CosmoTool +{ + + template + class ReplicateGenerator + { + public: + typedef Coord Coords[N]; + Coord replicate; + + ReplicateGenerator(const Coords& x, Coord shift) + { + face = 0; + replicate = shift; + numFaces = spower(3); + std::copy(x, x+N, x_base); + if (!next()) + abort(); + } + + bool next() + { + if (face == numFaces) + return false; + + face++; + + bool no_move = true; + int q_face = face; + for (int i = 0; i < N; i++) + { + int c_face; + c_face = q_face % 3; + q_face /= 3; + x_shifted[i] = x_base[i] + (c_face-1)*replicate; + no_move = no_move && (c_face == 1); + } + if (no_move) + return next(); + return true; + } + + const Coord *getPosition() + { + return x_shifted; + } + + void getPosition(Coords& x_out) + { + std::copy(x_shifted, x_shifted+N, x_out); + } + + private: + Coord x_shifted[N], x_base[N]; + long face, numFaces; + }; + +}; + +#endif diff --git a/src/sparseGrid.hpp b/src/sparseGrid.hpp index bbe3abf..516be9d 100644 --- a/src/sparseGrid.hpp +++ b/src/sparseGrid.hpp @@ -1,3 +1,38 @@ +/*+ +This is CosmoTool (./src/sparseGrid.hpp) -- Copyright (C) Guilhem Lavaux (2007-2013) + +guilhem.lavaux@gmail.com + +This software is a computer program whose purpose is to provide a toolbox for cosmological +data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...) + +This software is governed by the CeCILL license under French law and +abiding by the rules of distribution of free software. You can use, +modify and/ or redistribute the software under the terms of the CeCILL +license as circulated by CEA, CNRS and INRIA at the following URL +"http://www.cecill.info". + +As a counterpart to the access to the source code and rights to copy, +modify and redistribute granted by the license, users are provided only +with a limited warranty and the software's author, the holder of the +economic rights, and the successive licensors have only limited +liability. + +In this respect, the user's attention is drawn to the risks associated +with loading, using, modifying and/or developing or reproducing the +software by the user in light of its specific status of free software, +that may mean that it is complicated to manipulate, and that also +therefore means that it is reserved for developers and experienced +professionals having in-depth computer knowledge. Users are therefore +encouraged to load and test the software's suitability as regards their +requirements in conditions enabling the security of their systems and/or +data to be ensured and, more generally, to use and operate it in the +same conditions as regards security. + +The fact that you are presently reading this means that you have had +knowledge of the CeCILL license and that you accept its terms. ++*/ + #ifndef __VOID_SPARSE_GRID_HPP #define __VOID_SPARSE_GRID_HPP diff --git a/src/sphSmooth.hpp b/src/sphSmooth.hpp index 7d14b8a..f6fe582 100644 --- a/src/sphSmooth.hpp +++ b/src/sphSmooth.hpp @@ -1,3 +1,38 @@ +/*+ +This is CosmoTool (./src/sphSmooth.hpp) -- Copyright (C) Guilhem Lavaux (2007-2013) + +guilhem.lavaux@gmail.com + +This software is a computer program whose purpose is to provide a toolbox for cosmological +data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...) + +This software is governed by the CeCILL license under French law and +abiding by the rules of distribution of free software. You can use, +modify and/ or redistribute the software under the terms of the CeCILL +license as circulated by CEA, CNRS and INRIA at the following URL +"http://www.cecill.info". + +As a counterpart to the access to the source code and rights to copy, +modify and redistribute granted by the license, users are provided only +with a limited warranty and the software's author, the holder of the +economic rights, and the successive licensors have only limited +liability. + +In this respect, the user's attention is drawn to the risks associated +with loading, using, modifying and/or developing or reproducing the +software by the user in light of its specific status of free software, +that may mean that it is complicated to manipulate, and that also +therefore means that it is reserved for developers and experienced +professionals having in-depth computer knowledge. Users are therefore +encouraged to load and test the software's suitability as regards their +requirements in conditions enabling the security of their systems and/or +data to be ensured and, more generally, to use and operate it in the +same conditions as regards security. + +The fact that you are presently reading this means that you have had +knowledge of the CeCILL license and that you accept its terms. ++*/ + #ifndef __COSMOTOOL_SPH_SMOOTH_HPP #define __COSMOTOOL_SPH_SMOOTH_HPP diff --git a/src/yorick.hpp b/src/yorick.hpp index 25734f7..39fd20d 100644 --- a/src/yorick.hpp +++ b/src/yorick.hpp @@ -1,3 +1,38 @@ +/*+ +This is CosmoTool (./src/yorick.hpp) -- Copyright (C) Guilhem Lavaux (2007-2013) + +guilhem.lavaux@gmail.com + +This software is a computer program whose purpose is to provide a toolbox for cosmological +data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...) + +This software is governed by the CeCILL license under French law and +abiding by the rules of distribution of free software. You can use, +modify and/ or redistribute the software under the terms of the CeCILL +license as circulated by CEA, CNRS and INRIA at the following URL +"http://www.cecill.info". + +As a counterpart to the access to the source code and rights to copy, +modify and redistribute granted by the license, users are provided only +with a limited warranty and the software's author, the holder of the +economic rights, and the successive licensors have only limited +liability. + +In this respect, the user's attention is drawn to the risks associated +with loading, using, modifying and/or developing or reproducing the +software by the user in light of its specific status of free software, +that may mean that it is complicated to manipulate, and that also +therefore means that it is reserved for developers and experienced +professionals having in-depth computer knowledge. Users are therefore +encouraged to load and test the software's suitability as regards their +requirements in conditions enabling the security of their systems and/or +data to be ensured and, more generally, to use and operate it in the +same conditions as regards security. + +The fact that you are presently reading this means that you have had +knowledge of the CeCILL license and that you accept its terms. ++*/ + #ifndef __YORICK_HELPERS_HPP #define __YORICK_HELPERS_HPP diff --git a/src/yorick_nc3.cpp b/src/yorick_nc3.cpp index d760009..d80781e 100644 --- a/src/yorick_nc3.cpp +++ b/src/yorick_nc3.cpp @@ -1,3 +1,38 @@ +/*+ +This is CosmoTool (./src/yorick_nc3.cpp) -- Copyright (C) Guilhem Lavaux (2007-2013) + +guilhem.lavaux@gmail.com + +This software is a computer program whose purpose is to provide a toolbox for cosmological +data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...) + +This software is governed by the CeCILL license under French law and +abiding by the rules of distribution of free software. You can use, +modify and/ or redistribute the software under the terms of the CeCILL +license as circulated by CEA, CNRS and INRIA at the following URL +"http://www.cecill.info". + +As a counterpart to the access to the source code and rights to copy, +modify and redistribute granted by the license, users are provided only +with a limited warranty and the software's author, the holder of the +economic rights, and the successive licensors have only limited +liability. + +In this respect, the user's attention is drawn to the risks associated +with loading, using, modifying and/or developing or reproducing the +software by the user in light of its specific status of free software, +that may mean that it is complicated to manipulate, and that also +therefore means that it is reserved for developers and experienced +professionals having in-depth computer knowledge. Users are therefore +encouraged to load and test the software's suitability as regards their +requirements in conditions enabling the security of their systems and/or +data to be ensured and, more generally, to use and operate it in the +same conditions as regards security. + +The fact that you are presently reading this means that you have had +knowledge of the CeCILL license and that you accept its terms. ++*/ + #include "ctool_netcdf_ver.hpp" #include "config.hpp" #ifdef NETCDFCPP4 @@ -151,7 +186,7 @@ namespace CosmoTool { ProgressiveOutput::saveArrayProgressive(const std::string& fname, uint32_t *dimList, uint32_t rank) { - NcFile *f = new NcFile(fname.c_str(), NcFile::Replace); + NcFile *f = new NcFile(fname.c_str(), NcFile::Replace, 0, 0, NcFile::Netcdf4); assert(f->is_valid()); @@ -202,7 +237,7 @@ namespace CosmoTool { void saveArray(const std::string& fname, T *array, uint32_t *dimList, uint32_t rank) { - NcFile f(fname.c_str(), NcFile::Replace); + NcFile f(fname.c_str(), NcFile::Replace, 0, 0, NcFile::Netcdf4); assert(f.is_valid()); diff --git a/src/yorick_nc4.cpp b/src/yorick_nc4.cpp index 8bdeda0..9a6f501 100644 --- a/src/yorick_nc4.cpp +++ b/src/yorick_nc4.cpp @@ -1,3 +1,38 @@ +/*+ +This is CosmoTool (./src/yorick_nc4.cpp) -- Copyright (C) Guilhem Lavaux (2007-2013) + +guilhem.lavaux@gmail.com + +This software is a computer program whose purpose is to provide a toolbox for cosmological +data analysis (e.g. filters, generalized Fourier transforms, power spectra, ...) + +This software is governed by the CeCILL license under French law and +abiding by the rules of distribution of free software. You can use, +modify and/ or redistribute the software under the terms of the CeCILL +license as circulated by CEA, CNRS and INRIA at the following URL +"http://www.cecill.info". + +As a counterpart to the access to the source code and rights to copy, +modify and redistribute granted by the license, users are provided only +with a limited warranty and the software's author, the holder of the +economic rights, and the successive licensors have only limited +liability. + +In this respect, the user's attention is drawn to the risks associated +with loading, using, modifying and/or developing or reproducing the +software by the user in light of its specific status of free software, +that may mean that it is complicated to manipulate, and that also +therefore means that it is reserved for developers and experienced +professionals having in-depth computer knowledge. Users are therefore +encouraged to load and test the software's suitability as regards their +requirements in conditions enabling the security of their systems and/or +data to be ensured and, more generally, to use and operate it in the +same conditions as regards security. + +The fact that you are presently reading this means that you have had +knowledge of the CeCILL license and that you accept its terms. ++*/ + #include "ctool_netcdf_ver.hpp" #include "config.hpp" #include