Merge pull request #283 from LSSTDESC/non_limber_speedup

Non limber speedup

.travis.yml

@@ -37,6 +37,15 @@ before_install:
   - pip install coveralls
   - ./install-gsl.sh
   - python class_install.py
+  - python angpow_install.py
+#addons:
+#    apt:
+#        packages:
+#            - texlive-latex-recommended
+#            - texlive-latex-extra
+#            - texlive-fonts-recommended
+#            - texlive-fonts-extra
+#            - dvipng
 
 install: python setup.py install
 script:

Makefile.am

@@ -18,18 +18,37 @@ libccl_la_CPPFLAGS = -I./include/ -I./class/include -D__CCL_DATA_DIR__=$(datadir
 
 else
 
+libccl_la_LIBADD =
 libccl_la_CPPFLAGS = -I./include/ -I$(CLASSDIR)/include -D__CCL_DATA_DIR__=$(datadir)/ccl
 
 endif
 
+# To force libtool to use C++ compilers, a .cc files must be seen in the source files
+# To make the link between Angpow (C++) and CCL (C), this dummy file forces the use of g++
+if LOCALANGPOW
+
+ EXTRA_DIST += ./examples/ccl_sample_angpow.c
+ libccl_la_SOURCES += ./angpow/src/angpow_cpp_placeholder.cc 
+
+endif
+
+
 noinst_PROGRAMS = tests/check_ccl
 TESTS = tests/check_ccl tests/check_swig_timestamp.sh
 
-tests_check_ccl_SOURCES = tests/ccl_test.c tests/ccl_test_utils.c tests/ccl_test_params.c tests/ccl_test_params_mnu.c tests/ccl_test_cosmology.c tests/ccl_test_distances.c tests/ccl_test_distances_hiz.c tests/ccl_test_growth.c tests/ccl_test_growth_hiz.c tests/ccl_test_bbks.c tests/ccl_test_eh.c tests/ccl_test_cls.c tests/ccl_test_cmblens.c tests/ccl_test_sigmaM.c tests/ccl_test_massfunc.c tests/ccl_test_correlation.c tests/ccl_test_correlation_3d.c tests/ccl_test_bcm.c tests/ccl_test_emu.c tests/ccl_test_emu_nu.c tests/ccl_test_distances_mnu.c tests/ccl_test_distances_hiz_mnu.c tests/ccl_test_power_nu.c 
+tests_check_ccl_SOURCES = tests/ccl_test.c tests/ccl_test_utils.c tests/ccl_test_params.c tests/ccl_test_params_mnu.c tests/ccl_test_cosmology.c tests/ccl_test_distances.c tests/ccl_test_distances_hiz.c tests/ccl_test_growth.c tests/ccl_test_growth_hiz.c tests/ccl_test_bbks.c tests/ccl_test_eh.c tests/ccl_test_cls.c tests/ccl_test_cmblens.c tests/ccl_test_sigmaM.c tests/ccl_test_massfunc.c tests/ccl_test_correlation.c tests/ccl_test_correlation_3d.c tests/ccl_test_bcm.c tests/ccl_test_emu.c tests/ccl_test_emu_nu.c tests/ccl_test_distances_mnu.c tests/ccl_test_distances_hiz_mnu.c tests/ccl_test_power_nu.c tests/ccl_test_nonlimber.c
 
 tests_check_ccl_CPPFLAGS = -I./tests/ -I./include/
 tests_check_ccl_LDFLAGS =  -L./ -lccl -lgsl -lgslcblas -lfftw3 $(CLASS_LIB)
 
+# Add a specific unit test for Angpow
+if LOCALANGPOW
+
+ tests_check_ccl_SOURCES += tests/ccl_test_angpow.c
+
+endif
+
+
 dist-hook:
 	rm -f $(distdir)/clean_code.sh
 

README.md

@@ -32,6 +32,11 @@ This procedure has one final caveat: if you already have a working installation
 
 Note that, if you want to use your own version of CLASS, you should follow the steps described in the section "Compiling against an external version of CLASS" below.
 
+## Installing Angpow
+CCL provides an optional link to the Angpow library that enables fast and accurate computations of the angular power spectra without using the Limber approximation (written in C++). We provide a python script to automatically install Angpow and make the link with CCL. You should first run this script (`python angpow_install.py`) before carrying out the next steps. The installation downloads the last release of Angpow from this [DESC repository](https://github.com/LSSTDESC/Angpow4CCL) and looks for a C++ compiler compatible with OpenMP. It provides a dedicated file `./angpow/src/angpow_ccl.cc` that makes the interface between the CCL structures and the Angpow classes.
+
+To remove Angpow, run `python angpow_install.py --clean` and install CCL again.
+
 ## C-only installation
 Once the CLASS library is installed, `CCL` can be easily installed using an *autotools*-generated configuration file. To install `CCL`, from the base directory (the one where this file is located) run:
 ```sh
@@ -270,35 +275,54 @@ double ccl_sigma8(ccl_cosmology *cosmo, int * status);
 These and other functions for different matter power spectra can be found in file ***include/ccl_power.h***.
 
 ### Angular power spectra
-`CCL` can compute angular power spectra for two tracer types: galaxy number counts and galaxy weak lensing. Tracer parameters are defined in structure **`CCL_ClTracer`**. In general, you can create this object with function **`ccl_cl_tracer_new`**
+`CCL` can compute angular power spectra for three tracer types: galaxy number counts, galaxy weak lensing and CMB lensing. Tracer parameters are defined in structure **`CCL_ClTracer`**. In general, you can create this object with function **`ccl_cl_tracer`**
 ````c
-CCL_ClTracer *ccl_cl_tracer_new(ccl_cosmology *cosmo,int tracer_type,
+CCL_ClTracer *ccl_cl_tracer(ccl_cosmology *cosmo,int tracer_type,
 				int has_rsd,int has_magnification,int has_intrinsic_alignment,
 				int nz_n,double *z_n,double *n,
 				int nz_b,double *z_b,double *b,
 				int nz_s,double *z_s,double *s,
 				int nz_ba,double *z_ba,double *ba,
-				int nz_rf,double *z_rf,double *rf, int * status);
+				int nz_rf,double *z_rf,double *rf, 
+				double z_source, int * status);
 ````
-Exact definition of these parameters are described in file ***include/ccl_cls.h***. Usually you can use simplified versions of this function, namely **`ccl_cl_tracer_number_counts_new`, `ccl_cl_tracer_number_counts_simple_new`, `ccl_cl_tracer_lensing_new`** or **`ccl_cl_tracer_lensing_simple_new`**. Two most simplified versions (one for number counts and one for shear) take parameters:
+Exact definition of these parameters are described in file ***include/ccl_cls.h***. Usually you can use simplified versions of this function, namely **`ccl_cl_tracer_number_counts`, `ccl_cl_tracer_number_counts_simple`, `ccl_cl_tracer_lensing`, `ccl_cl_tracer_lensing_simple`** or **`ccl_cl_tracer_cmblens`**. Two most simplified versions (one for number counts and one for shear) take parameters:
 ````c
-CCL_ClTracer *ccl_cl_tracer_number_counts_simple_new(ccl_cosmology *cosmo,
+CCL_ClTracer *ccl_cl_tracer_number_counts_simple(ccl_cosmology *cosmo,
 						     int nz_n,double *z_n,double *n,
-                             int nz_b,double *z_b,double *b, int * status);
-CCL_ClTracer *ccl_cl_tracer_lensing_simple_new(ccl_cosmology *cosmo,
+                                                     int nz_b,double *z_b,double *b, int * status);
+CCL_ClTracer *ccl_cl_tracer_lensing_simple(ccl_cosmology *cosmo,
 					       int nz_n,double *z_n,double *n, int * status);
 ````
 where `nz_n` is dimension (number of bins) of arrays `z_n` and `n`. `z_n` and `n` are arrays for the number count of objects per redshift interval (arbitrary normalization - renormalized inside). `nz_b`, `z_b` and `b` are the same for the clustering bias.
 
-With initialized tracers you can compute limber power spectrum with **`ccl_angular_cl`**
+Before computing the angular power spectrum, users must define a workspace structure that contains the relevant parameters for the computation:
 ````c
-double ccl_angular_cl(ccl_cosmology *cosmo,int l,CCL_ClTracer *clt1,CCL_ClTracer *clt2, int * status);
+CCL_ClWorkspace *ccl_cl_workspace_default(int lmax,int l_limber,int non_limber_method,
+					  double l_logstep,int l_linstep,
+					  double dchi,double dlk,double zmin,int *status)
 ````
-After you are done working with tracers, you should free its work space by **`ccl_cl_tracer_free`**
+where `lmax` sets the maximum multipole, `l_limber` the limit multipole from which the Limber approximation is used (`l_limber=-1` means that the Liber approximation is never used). The `non_limber_method` variable can be set to `CCL_NONLIMBER_METHOD_NATIVE` or `CCL_NONLIMBER_METHOD_ANGPOW` to choose the method to compute the non-Limber part of the angular power spectrum (either the native `CCL` code or the [`Angpow` library](https://github.com/LSSTDESC/CCL/blob/non_limber_speedup/README.md#installing-angpow)). Then `l_linstep` sets the maximum multipole until which the angular power spectrum is computed at each multipole, and `l_logstep` the logarithmic stepping to use above `l_linstep` (then the power spectrum is interpolated at each multipole). `dchi` sets the interval in comoving distance to use for the native non-Limber computation and `dlk`the logarithmic stepping for the Fourier k-integration (`Angpow` is not concerned by these two parameters). A simplified workspace is provided for computations that use only the Limber approximation at each multipole:
+````c
+CCL_ClWorkspace *ccl_cl_workspace_default_limber(int lmax,double l_logstep,int l_linstep,
+						 double dlk,int *status)
+````
+
+With initialized tracers and workspace you can compute limber power spectrum with **`ccl_angular_cls`**
+````c
+double ccl_angular_cls(ccl_cosmology *cosmo,int l,CCL_ClTracer *clt1,CCL_ClTracer *clt2, 
+						int nl_out,int *l_out,double *cl_out,int * status);
+````
+with `l_out` and `cl_out` arrays of size `nl_out` that contains the multipoles and the angular power spectrum. 
+
+After you are done working with tracers, you should free its work space by **`ccl_cl_tracer_free`** and **`ccl_cl_workspace_free`** 
 ````c
 void ccl_cl_tracer_free(CCL_ClTracer *clt);
+void ccl_cl_workspace_free(CCL_ClWorkspace *w);
 ````
 
+Note that for the moment `Angpow` can not handle the magnification lensing term for the galaxy number count tracers, and has not been tested for the weak lensing tracer. This limitations will be removed in the near future.
+
 ### Halo mass function
 The halo mass function *dN/dM* can be obtained by function **`ccl_massfunc`**
 ````c