Merge pull request #371 from LSSTDESC/gsl_params_struct

User-changable GSL parameters
LSSTDESC · May 29, 2018 · 46857fd · 46857fd
2 parents d82725b + d008341
commit 46857fd
Show file tree

Hide file tree

Showing 24 changed files with 4,141 additions and 3,949 deletions.
diff --git a/doc/0000-ccl_note/0000-ccl_note.pdf b/doc/0000-ccl_note/0000-ccl_note.pdf
diff --git a/doc/0000-ccl_note/authors.csv b/doc/0000-ccl_note/authors.csv
@@ -58,6 +58,7 @@ Plaszczynski,St\'ephane,St\'ephane Plaszczynski,Contributor,"Laboratoire de l'Ac
 Sanchez,Javier,Javier Sanchez,Contributor,"Department of Physics and Astronomy, University of California, Irvine, CA 92697, USA",Modified setup.py to allow pip installation and uninstall.,francs1@uci.edu
 Singh,Sukhdeep,Sukhdeep Singh,Contributor,"McWilliams Center for Cosmology, Department of Physics, Carnegie Mellon University, Pittsburgh, PA 15213, USA",Contributed to the correlation functions code.,sukhdeep@cmu.edu
 Slosar,Anze,An\v{z}e Slosar,Contributor,"Brookhaven National Laboratory, Physics Department, Upton, NY 11973, USA",Wrote and reviewed code.,anze@bnl.gov
+Tr\"oster,Tilman,Tilman Tr\"oster,Contributor,"Institute for Astronomy, Royal Observatory Edinburgh, Edinburgh EH9 3HJ, UK","Wrote code for user-changable precision parameters, found and fixed bugs.",ttr@roe.ac.uk
 Villarreal,Antonio,Antonio Villarreal,Contributor,"Department of Physics and Astronomy, University of Pittsburgh, Pittsburgh PA 15260","Contributed to initial benchmarking, halo mass function code, and general code and issues review.",asv13@pitt.edu
 Vrastil,Michal,Michal Vrastil,Contributor,"Institute of Physics CAS, Prague, 182 21, CZ","Wrote documentation and example code, reviewed code.",vrastil@fzu.cz
 Zuntz,Joe,Joe Zuntz,Contributor,"Institute for Astronomy, Royal Observatory Edinburgh, Edinburgh EH9 3HJ, UK","Wrote initial infrastructure, C testing setup, and reviewed code.",joezuntz@googlemail.com
diff --git a/doc/0000-ccl_note/authors.tex b/doc/0000-ccl_note/authors.tex
@@ -56,6 +56,8 @@
 \affiliation{McWilliams Center for Cosmology, Department of Physics, Carnegie Mellon University, Pittsburgh, PA 15213, USA}
 \author{An\v{z}e Slosar}
 \affiliation{Brookhaven National Laboratory, Physics Department, Upton, NY 11973, USA}
+\author{Tilman Tr\"oster}
+\affiliation{Institute for Astronomy, Royal Observatory Edinburgh, Edinburgh EH9 3HJ, UK}
 \author{Antonio Villarreal}
 \affiliation{Department of Physics and Astronomy, University of Pittsburgh, Pittsburgh PA 15260}
 \author{Michal Vrastil}

diff --git a/doc/0000-ccl_note/contributions.tex b/doc/0000-ccl_note/contributions.tex
@@ -25,6 +25,7 @@
 Javier Sanchez: Modified setup.py to allow pip installation and uninstall. \\
 Sukhdeep Singh: Contributed to the correlation functions code. \\
 An\v{z}e Slosar: Wrote and reviewed code. \\
+Tilman Tr\"oster: Wrote code for user-changable precision parameters, found and fixed bugs. \\
 Antonio Villarreal: Contributed to initial benchmarking, halo mass function code, and general code and issues review. \\
 Michal Vrastil: Wrote documentation and example code, reviewed code. \\
 Joe Zuntz: Wrote initial infrastructure, C testing setup, and reviewed code. \\
diff --git a/doc/0000-ccl_note/main.pdf b/doc/0000-ccl_note/main.pdf
diff --git a/doc/0000-ccl_note/main.tex b/doc/0000-ccl_note/main.tex
@@ -480,6 +480,29 @@ \subsubsection{Spline parameters \& the INI file}
 \item {\tt N$\_$K$\_$3DCOR}: Number of bins per decade for the wavenumber in the 3d correlation function calculation. 
 The default is {\tt N$\_$K$\_$3DCOR}~=~100,000.
 \end{itemize}
+
+The precision parameters of {\tt GSL} routines can also be set in {\tt ccl$\_$params.ini}. The parameters are organized in a hierarchal fashion, such that more specialized parameters take precedence over more general ones. For example, there is a general tolerance parameter, a tolerance parameter for integrals, and a tolerance parameter for the Limber integral. By default all are set to the value of the general tolerance. Changing the integration tolerance changes the tolerance for all integrals but not for root-finding routines, for example, and changing the Limber integral tolerance only affects that integral and no other.
+The parameters are:
+
+\begin{itemize}
+\item {\tt GSL$\_$EPSREL}: Relative tolerance for {\tt GSL} routines. Default is $10^{-4}$.
+\item {\tt GSL$\_$N$\_$ITERATION}: Maximum number of iterations allowed for {\tt GSL} routines. This parameter also controls the size of the integration workspaces. Default is 1000.
+\item {\tt GSL$\_$INTEGRATION$\_$GAUSS$\_$KRONROD$\_$POINTS}: Number of Gauss-Kronrod points in the {\tt QAG} integration routines. Default is {\tt GSL$\_$INTEG$\_$GAUSS41} (numerical value 4).
+\item {\tt GSL$\_$INTEGRATION$\_$EPSREL}: Relative tolerance for {\tt GSL} integration routines. Default is {\tt GSL$\_$EPSREL}.
+\item {\tt GSL$\_$INTEGRATION$\_$DISTANCE$\_$EPSREL}: Relative tolerance for {\tt GSL} integration routines for distance calculations. Default is $10^{-6}$.
+\item {\tt GSL$\_$INTEGRATION$\_$DNDZ$\_$EPSREL}: Relative tolerance for {\tt GSL} integration routines for  $\frac{dN}{dz}$ calculations. Default is $10^{-6}$.
+\item {\tt GSL$\_$INTEGRATION$\_$SIGMAR$\_$EPSREL}: Relative tolerance for {\tt GSL} integration routines for  $\sigma_R$ calculations. Default is $10^{-5}$.
+\item {\tt GSL$\_$INTEGRATION$\_$NU$\_$EPSREL}: Relative tolerance for {\tt GSL} integration routines for neutrino calculations. Default is $10^{-7}$.
+\item {\tt GSL$\_$INTEGRATION$\_$NU$\_$EPSABS}: Absolute tolerance for {\tt GSL} integration routines for neutrino calculations. Default is $10^{-7}$.
+\item {\tt GSL$\_$INTEGRATION$\_$LIMBER$\_$GAUSS$\_$KRONROD$\_$POINTS}: Number of Gauss-Kronrod points in the {\tt QAG} integration routines for the Limber integral. Default is {\tt GSL$\_$INTEGRATION$\_$GAUSS$\_$KRONROD$\_$POINTS}.
+\item {\tt GSL$\_$INTEGRATION$\_$LIMBER$\_$EPSREL}: Relative tolerance for {\tt GSL} integration routines for Limber integral. Default is {\tt GSL$\_$EPSREL}.
+\item {\tt GSL$\_$ROOT$\_$EPSREL}: Relative tolerance for {\tt GSL} root-finding routines. Default is {\tt GSL$\_$EPSREL}.
+\item {\tt GSL$\_$ROOT$\_$N$\_$ITERATION}: Maximum number of iterations allowed for {\tt GSL} root-finding routines. Default is {\tt GSL$\_$N$\_$ITERATION}.
+\item {\tt GSL$\_$ODE$\_$GROWTH$\_$EPSREL}: Relative tolerance for {\tt GSL} ODE-solving routines for growth factor calculation. Default is $10^{-6}$.
+\end{itemize}
+
+The {\tt GSL} parameters are initialized to their default values and do not need to be listed in {\tt ccl$\_$params.ini} if no adjustments are necessary.
+
 Note that a copy of {\tt ccl$\_$params.ini} is installed along with the library, so changing the version of this file inside the source directory will not have any effect unless you reinstall.
 
 For the matter power spectrum, the spline is performed in two variables: the logarithmically-spaced wavenumber and the linearly-spaced scale factor. Splining the {\tt CLASS} output leads to some precision loss (compared to direct outputs from {\tt CLASS}). We quantify this, along with the impact of extrapolation, in the following subsection.

diff --git a/include/ccl_constants.h b/include/ccl_constants.h
@@ -104,30 +104,63 @@ extern "C" {
 
 
 //Precision parameters
+/**
+ * Default relative precision if not otherwise specified
+ */
+#define GSL_EPSREL 1E-4
+
+/**
+ * Default number of iterations for integration and root-finding if not otherwise
+ * specified
+ */
+#define GSL_N_ITERATION 1000
+
+/**
+ * Default number of Gauss-Kronrod points in QAG integration if not otherwise 
+ * specified
+ */
+#define GSL_INTEGRATION_GAUSS_KRONROD_POINTS GSL_INTEG_GAUSS41
+
+/**
+ * Absolute precision in neutrino root finding
+ */
+#define GSL_EPSABS_NU 1E-7
+
+/**
+ * Relative precision in neutrino root finding
+ */
+#define GSL_EPSREL_NU 1E-7
+
+/**
+ * Number of iterations for neutrino root finding
+ */
+#define GSL_N_ITERATION_NU 1000
+
+/**
+ * Relative precision in sigma_R calculations
+ */
+#define GSL_EPSREL_SIGMAR 1E-5
+
 /**
  * Relative precision in distance calculations
  */
-#define EPSREL_DIST 1E-6
+#define GSL_EPSREL_DIST 1E-6
 
 /**
  * Relative precision in growth calculations
  */
-#define EPSREL_GROWTH 1E-6
+#define GSL_EPSREL_GROWTH 1E-6
 
 /**
  * Relative precision in dNdz calculations
  */
-#define EPSREL_DNDZ 1E-6
+#define GSL_EPSREL_DNDZ 1E-6
 
 /**
  * Absolute precision in growth calculations
  */
 #define EPS_SCALEFAC_GROWTH 1E-6
 
-//Correlation function related parameters
-#define EPSREL_CORR_FUNC 1E-3
-#define GSL_INTEGRATION_LIMIT 1000
-
 //LSST specific numbers
 #define Z_MIN_SOURCES 0.1
 #define Z_MAX_SOURCES 3.0

diff --git a/include/ccl_error.h b/include/ccl_error.h
@@ -36,6 +36,7 @@ typedef enum {
 typedef enum {
   CCL_DEBUG_MODE_OFF = 0,
   CCL_DEBUG_MODE_ON = 1,
+  CCL_DEBUG_MODE_WARNING = 2,
 } CCLDebugModePolicy;
 
 /** Raise an exception
@@ -44,6 +45,18 @@ typedef enum {
  */
 void ccl_raise_exception(int err, char* msg);
 
+/** Raise a warning
+ * Given a status, give a warning message.
+ * @return void
+ */
+void ccl_raise_warning(int err, char* msg);
+
+/** Raise a warning based on a GSL error message
+ * Given a GSL status, give a warning message.
+ * @return void
+ */
+void ccl_raise_gsl_warning(int gslstatus, char* msg);
+
 /** Set the error policy
  * @oaram error_policy the error policy
  * @return void

diff --git a/include/ccl_params.h b/include/ccl_params.h
@@ -40,6 +40,43 @@ typedef struct ccl_spline_params {
 
 extern ccl_spline_params * ccl_splines;
 
+/** 
+ * Struct that contains parameters that control the accuracy of various GSL
+ * routines.
+ */
+typedef struct ccl_gsl_params {
+  // General parameters. If not otherwise specified, those will be copied to the
+  // more specialised cases.
+  double EPSREL;
+  size_t N_ITERATION;
+
+  // Integration
+  int INTEGRATION_GAUSS_KRONROD_POINTS;
+  double INTEGRATION_EPSREL;
+  // Limber integration
+  int INTEGRATION_LIMBER_GAUSS_KRONROD_POINTS;
+  double INTEGRATION_LIMBER_EPSREL;
+  // Distance integrals
+  double INTEGRATION_DISTANCE_EPSREL;
+  // dndz integrals
+  double INTEGRATION_DNDZ_EPSREL;
+  // sigma_R integral
+  double INTEGRATION_SIGMAR_EPSREL;
+  // Neutrino integral
+  double INTEGRATION_NU_EPSREL;
+  double INTEGRATION_NU_EPSABS;
+
+  // Root finding
+  double ROOT_EPSREL;
+  int ROOT_N_ITERATION;
+
+  // ODE
+  double ODE_GROWTH_EPSREL;
+
+} ccl_gsl_params;
+
+extern ccl_gsl_params * ccl_gsl;
+
 #ifdef __cplusplus
 }
 #endif
diff --git a/pyccl/ccl.i b/pyccl/ccl.i
@@ -39,6 +39,7 @@
 %include "ccl_constants.i"
 %include "ccl_lsst_specs.i"
 %include "ccl_neutrinos.i"
+%include "ccl_params.i"
 
 %include "../include/ccl_config.h"
 %include "../include/ccl_error.h"

diff --git a/pyccl/ccl_params.i b/pyccl/ccl_params.i
@@ -0,0 +1,14 @@
+%module ccl_params
+
+%{
+#define SWIG_FILE_WITH_INIT
+#include "../include/ccl_params.h"
+%}
+
+// Automatically document arguments and output types of all functions
+%feature("autodoc", "1");
+
+// Strip the ccl_ prefix from function names
+%rename("%(strip:[ccl_])s") "";
+
+%include "../include/ccl_params.h"