diff --git a/INSTALL.md b/INSTALL.md index 88c77b77778..4c3d552974e 100644 --- a/INSTALL.md +++ b/INSTALL.md @@ -63,17 +63,12 @@ The GalSim package also requires functionality, and can be omitted if users are happy to use JSON-style config parsing or prefer to write python scripts directly. -* Optional dependency: the scientific Python module SciPy - (http://www.scipy.org). SciPy is only required for generation of - shears due to an NFW halo. - These should installed onto your Python system so that they can be imported by: >>> import numpy >>> import pyfits >>> import yaml [ if using the galsim_yaml executable or otherwise plan to parse .yaml configuration files ] - >>> import scipy [ if generating shears according to an NFW halo ] within Python. You can test this by loading up the Python interpreter for the version of Python you will be using with the GalSim toolkit. This is usually diff --git a/galsim/__init__.py b/galsim/__init__.py index 87ad01b5a5d..5692aac51e3 100644 --- a/galsim/__init__.py +++ b/galsim/__init__.py @@ -16,3 +16,4 @@ from psfcorr import * from io import * from . import lensing +from . import integ diff --git a/galsim/integ.py b/galsim/integ.py new file mode 100644 index 00000000000..d6b372cd71b --- /dev/null +++ b/galsim/integ.py @@ -0,0 +1,36 @@ +"""@file integ.py +A Python layer version of the C++ int1d function in galim::integ +""" + +from . import _galsim + +def int1d(func, min, max, rel_err=1.e-6, abs_err=1.e-12): + """Integrate a 1-dimensional function from min to max. + + Example usage: + + >>> def func(x): return x**2 + >>> galsim.integ.int1d(func, 0, 1) + 0.33333333333333337 + >>> galsim.integ.int1d(func, 0, 2) + 2.666666666666667 + >>> galsim.integ.int1d(func, -1, 1) + 0.66666666666666674 + + @param func The function to be integrated. y = func(x) should be valid. + @param min The lower end of the integration bounds (anything < -1.e10 is treated as + negative infinity). + @param max The upper end of the integration bounds (anything > 1.e10 is treated as positive + infinity). + @param rel_err The desired relative error (default `rel_err = 1.e-6`) + @param abs_err The desired absolute error (default `abs_err = 1.e-12`) + """ + min = float(min) + max = float(max) + rel_err = float(rel_err) + abs_err = float(abs_err) + success, result = _galsim.PyInt1d(func,min,max,rel_err,abs_err) + if success: + return result + else: + raise RuntimeError(result) diff --git a/galsim/lensing.py b/galsim/lensing.py index d58133bdd89..9f5a807d118 100644 --- a/galsim/lensing.py +++ b/galsim/lensing.py @@ -414,9 +414,8 @@ def Da(self, z, z_ref=0): raise ValueError("Redshift z must not be negative") if z < z_ref: raise ValueError("Redshift z must not be smaller than the reference redshift") - # TODO: We don't want to depend on scipy, so need to move this down to c++. - from scipy.integrate import quad - d = quad(self.__angKernel, z_ref+1, z+1)[0] + + d = galsim.integ.int1d(self.__angKernel, z_ref+1, z+1) # check for curvature rk = (abs(self.omega_c))**0.5 if (rk*d > 0.01): diff --git a/include/galsim/integ/Int.h b/include/galsim/integ/Int.h index 0de6106d8ce..7944e9fa12c 100644 --- a/include/galsim/integ/Int.h +++ b/include/galsim/integ/Int.h @@ -51,7 +51,7 @@ * * There are two final arguments, which we've omitted so far, that can be used to specify * the precision required. First the relative error, then the absolute error. - * The defaults are 1.e-6 and 1.e-15 respectively, which are generally fine for most + * The defaults are 1.e-6 and 1.e-12 respectively, which are generally fine for most * purposes, but you can specify different values if you prefer. * * The absolute error only comes into play for results which are close to @@ -148,16 +148,12 @@ namespace integ { const double MOCK_INF = 1.e100; ///< May be used to indicate infinity in integration regions. const double MOCK_INF2 = 1.e10; ///< Anything larger than this is treated as infinity. const double DEFRELERR = 1.e-6; ///< The default target relative error if not specified. - const double DEFABSERR = 1.e-15; ///< The default target absolute error if not specified. + const double DEFABSERR = 1.e-12; ///< The default target absolute error if not specified. /// An exception type thrown if the integrator encounters a problem. struct IntFailure : public std::runtime_error - { - IntFailure(const std::string& s) : - std::runtime_error(s.c_str()) - {} - }; + { IntFailure(const std::string& s) : std::runtime_error(s) {} }; #ifdef NDEBUG #define integ_dbg1 if (false) (*dbgout) @@ -334,9 +330,6 @@ namespace integ { } } - // All code between the @cond and @endcond is excluded from Doxygen documentation - //! @cond - /** * @brief Add a split point to the current list to be used by the next subDivide call * @@ -366,7 +359,6 @@ namespace integ { /// Setup an fxmap for this region. void useFXMap() { _fxmap_source.reset(new std::map()); fxmap = _fxmap_source.get(); } - //! @endcond private: T _a,_b,_error,_area; @@ -380,127 +372,86 @@ namespace integ { boost::shared_ptr > _fxmap_source; }; - // All code between the @cond and @endcond is excluded from Doxygen documentation - //! @cond - - /// Rescale the error if int |f| dx or int |f-mean| dx are too large - template - inline T rescaleError( - T err, ///< The current estimate of the error - const T& int_abs, ///< An estimate of int |f| dx - const T& int_absdiff ///< An estimate of int |f-mean| dx - ) - { - const T eps = std::numeric_limits::epsilon(); - const T minrep = std::numeric_limits::min(); - - if (int_absdiff != 0. && err != 0.) { - const T scale = (200. * err / int_absdiff); - if (scale < 1.) err = int_absdiff * scale * sqrt(scale) ; - else err = int_absdiff ; - } - if (int_abs > minrep / (50. * eps)) { - const T min_err = 50. * eps * int_abs; - if (min_err > err) err = min_err; - } - return err; - } + namespace { + /// Rescale the error if int |f| dx or int |f-mean| dx are too large + template + inline T rescaleError( + T err, ///< The current estimate of the error + const T& int_abs, ///< An estimate of int |f| dx + const T& int_absdiff ///< An estimate of int |f-mean| dx + ) + { + const T eps = std::numeric_limits::epsilon(); + const T minrep = std::numeric_limits::min(); - /** - * @brief Non-adaptive GKP integration - * - * A non-adaptive integration of the function f over the region reg. - * - * The algorithm computes first a Gaussian quadrature value - * then successive Kronrod/Patterson extensions to this result. - * The functions terminates when the difference between successive - * approximations (rescaled according to rescaleError) is less than - * either abserr or relerr * I, where I is the latest estimate of the - * integral. - * - * The order of the Gauss/Kronron/Patterson scheme is determined - * by which file is included above. Currently schemes starting - * with order 1 and order 10 are calculated. There seems to be - * little practical difference in the integration times using - * the two schemes, so I haven't bothered to calculate any more. - */ - template - inline bool intGKPNA( - const UF& func, ///< The function to integrate - IntRegion& reg, ///< The region with the bounds - const typename UF::result_type relerr, ///< The target relative error - const typename UF::result_type abserr ///< The target absolute error - ) - { - typedef typename UF::result_type T; - const T a = reg.left(); - const T b = reg.right(); - - const T half_length = 0.5 * (b - a); - const T abs_half_length = std::abs(half_length); - const T center = 0.5 * (b + a); - const T f_center = func(center); - if (reg.fxmap) (*reg.fxmap)[center] = f_center; -#ifdef COUNTFEVAL - nfeval++; -#endif - const int nmax = 2*gkp_x(NGKPLEVELS-1).size()-1; - static std::vector fv1(nmax), fv2(nmax); - - fv1.clear(); - fv2.clear(); - assert(fv1.size() == 0); - assert(fv2.size() == 0); - assert(int(fv1.capacity()) == nmax); - assert(int(fv2.capacity()) == nmax); - - assert(gkp_wb(0).size() == gkp_x(0).size()+1); - T area1 = gkp_wb(0).back() * f_center; - int n0 = gkp_x(0).size(); - for (int k=0; k(0)[k]; - const T fval1 = func(center - abscissa); - const T fval2 = func(center + abscissa); - area1 += gkp_wb(0)[k] * (fval1+fval2); - fv1.push_back(fval1); - fv2.push_back(fval2); - if (reg.fxmap) { - (*reg.fxmap)[center-abscissa] = fval1; - (*reg.fxmap)[center+abscissa] = fval2; + if (int_absdiff != 0. && err != 0.) { + const T scale = (200. * err / int_absdiff); + if (scale < 1.) err = int_absdiff * scale * sqrt(scale) ; + else err = int_absdiff ; + } + if (int_abs > minrep / (50. * eps)) { + const T min_err = 50. * eps * int_abs; + if (min_err > err) err = min_err; } + return err; } - area1 *= half_length; + + /** + * @brief Non-adaptive GKP integration + * + * A non-adaptive integration of the function f over the region reg. + * + * The algorithm computes first a Gaussian quadrature value + * then successive Kronrod/Patterson extensions to this result. + * The functions terminates when the difference between successive + * approximations (rescaled according to rescaleError) is less than + * either abserr or relerr * I, where I is the latest estimate of the + * integral. + * + * The order of the Gauss/Kronron/Patterson scheme is determined + * by which file is included above. Currently schemes starting + * with order 1 and order 10 are calculated. There seems to be + * little practical difference in the integration times using + * the two schemes, so I haven't bothered to calculate any more. + */ + template + inline bool intGKPNA( + const UF& func, ///< The function to integrate + IntRegion& reg, ///< The region with the bounds + const typename UF::result_type relerr, ///< The target relative error + const typename UF::result_type abserr ///< The target absolute error + ) + { + typedef typename UF::result_type T; + const T a = reg.left(); + const T b = reg.right(); + + const T half_length = 0.5 * (b - a); + const T abs_half_length = std::abs(half_length); + const T center = 0.5 * (b + a); + const T f_center = func(center); + if (reg.fxmap) (*reg.fxmap)[center] = f_center; #ifdef COUNTFEVAL - nfeval+=gkp_x(0).size()*2; + nfeval++; #endif - - integ_dbg2<<"level 0 rule: area = "<(level).size() == fv1.size()); - assert(gkp_wa(level).size() == fv2.size()); - assert(gkp_wb(level).size() == gkp_x(level).size()+1); - T area2 = gkp_wb(level).back() * f_center; - // int_abs = approximation to integral of abs(f) - if (calc_int_abs) int_abs = std::abs(area2); - for (size_t k=0; k(level)[k] * (fv1[k]+fv2[k]); - if (calc_int_abs) - int_abs += gkp_wa(level)[k] * - (std::abs(fv1[k]) + std::abs(fv2[k])); - } - int nl = gkp_x(level).size(); - for (int k=0; k(level)[k]; + const int nmax = 2*gkp_x(NGKPLEVELS-1).size()-1; + static std::vector fv1(nmax), fv2(nmax); + + fv1.clear(); + fv2.clear(); + assert(fv1.size() == 0); + assert(fv2.size() == 0); + assert(int(fv1.capacity()) == nmax); + assert(int(fv2.capacity()) == nmax); + + assert(gkp_wb(0).size() == gkp_x(0).size()+1); + T area1 = gkp_wb(0).back() * f_center; + int n0 = gkp_x(0).size(); + for (int k=0; k(0)[k]; const T fval1 = func(center - abscissa); const T fval2 = func(center + abscissa); - const T fval = fval1 + fval2; - area2 += gkp_wb(level)[k] * fval; - if (calc_int_abs) - int_abs += gkp_wb(level)[k] * (std::abs(fval1) + std::abs(fval2)); + area1 += gkp_wb(0)[k] * (fval1+fval2); fv1.push_back(fval1); fv2.push_back(fval2); if (reg.fxmap) { @@ -508,271 +459,304 @@ namespace integ { (*reg.fxmap)[center+abscissa] = fval2; } } + area1 *= half_length; #ifdef COUNTFEVAL - nfeval+=gkp_x(level).size()*2; + nfeval+=gkp_x(0).size()*2; #endif - if (calc_int_abs) { - const T mean = area1*T(0.5); - // int_absdiff = approximation to the integral of abs(f-mean) - int_absdiff = gkp_wb(level).back() * std::abs(f_center-mean); - for (size_t k=0; k(level).size(); k++) { - int_absdiff += gkp_wa(level)[k] * - (std::abs(fv1[k]-mean) + std::abs(fv2[k]-mean)); + + integ_dbg2<<"level 0 rule: area = "<(level).size() == fv1.size()); + assert(gkp_wa(level).size() == fv2.size()); + assert(gkp_wb(level).size() == gkp_x(level).size()+1); + T area2 = gkp_wb(level).back() * f_center; + // int_abs = approximation to integral of abs(f) + if (calc_int_abs) int_abs = std::abs(area2); + for (size_t k=0; k(level)[k] * (fv1[k]+fv2[k]); + if (calc_int_abs) + int_abs += gkp_wa(level)[k] * + (std::abs(fv1[k]) + std::abs(fv2[k])); } - for (size_t k=0; k(level).size(); k++) { - int_absdiff += gkp_wb(level)[k] * - (std::abs(fv1[k]-mean) + std::abs(fv2[k]-mean)); + int nl = gkp_x(level).size(); + for (int k=0; k(level)[k]; + const T fval1 = func(center - abscissa); + const T fval2 = func(center + abscissa); + const T fval = fval1 + fval2; + area2 += gkp_wb(level)[k] * fval; + if (calc_int_abs) + int_abs += gkp_wb(level)[k] * (std::abs(fval1) + std::abs(fval2)); + fv1.push_back(fval1); + fv2.push_back(fval2); + if (reg.fxmap) { + (*reg.fxmap)[center-abscissa] = fval1; + (*reg.fxmap)[center+abscissa] = fval2; + } } - int_absdiff *= abs_half_length ; - int_abs *= abs_half_length; - } - area2 *= half_length; - err = rescaleError(std::abs(area2-area1), int_abs, int_absdiff) ; - if (err < int_absdiff) calc_int_abs = false; - - integ_dbg2<<"at level "< - inline void intGKP( - const UF& func, IntRegion& reg, - const typename UF::result_type relerr, - const typename UF::result_type abserr) - { - typedef typename UF::result_type T; - const T eps = std::numeric_limits::epsilon(); - - integ_dbg2<<"Start intGKP\n"; - - assert(abserr >= 0.); - assert(relerr > 0.); - - // Check for early exit: - if (reg.left() == reg.right()) { - integ_dbg2<<"left == right, so integral is trivially 0.\n"; - reg.setArea(0.,0.); - return; + return false; } - // Perform the first integration - bool done = intGKPNA(func, reg, relerr, abserr); - if (done) { - integ_dbg2<<"GKPNA suceeded, so we're done.\n"; - return; - } + /** + * @brief Adaptive GKP integration + * + * An adaptive integration algorithm which computes the integral of f + * over the region reg. + * + * First the non-adaptive GKP algorithm is tried. + * + * If that is not accurate enough (according to the absolute and + * relative accuracies, abserr and relerr), the region is split in half, + * and each new region is integrated. + * + * The routine continues by successively splitting the subregion + * which gave the largest absolute error until the integral converges. + * + * The area and estimated error are returned as reg.getArea() and reg.getErr() + */ + template + inline void intGKP( + const UF& func, IntRegion& reg, + const typename UF::result_type relerr, + const typename UF::result_type abserr) + { + typedef typename UF::result_type T; + const T eps = std::numeric_limits::epsilon(); - integ_dbg2<<"In adaptive GKP, failed first pass... subdividing\n"; - integ_dbg2<<"Intial range = "< for IntRegoins such that the "largest" is the one - // with the largest current error estimate. This is the next one to be split - // if we need to split further. - std::priority_queue,std::vector > > allregions; - allregions.push(reg); - T finalarea = reg.getArea(); - T finalerr = reg.getErr(); - T tolerance= std::max(abserr, relerr * std::abs(finalarea)); - assert(finalerr > tolerance); - - while(!error_type && finalerr > tolerance) { - // Bisect the subinterval with the largest error estimate - integ_dbg2<<"Current answer = "< parent = allregions.top(); - allregions.pop(); - integ_dbg2<<"Subdividing largest error region "; - integ_dbg2< > children; - parent.subDivide(children); - // For "GKP", there are only two, but for GKPOSC, there is one - // for each oscillation in region - - // Try to do at least 3x better with the children - T factor = 3*children.size()*finalerr/tolerance; - T newabserr = std::abs(parent.getErr()/factor); - T newrelerr = newabserr/std::abs(parent.getArea()); - integ_dbg2<<"New abserr,rel = "<& child = children[i]; - integ_dbg2<<"Integrating child "<= 10 && newerror > parent.getErr() && - std::abs(delta) <= newerror-parent.getErr()) { - integ_dbg2<<"roundoff type 2: newerror/error = "; - integ_dbg2< tolerance) { - if (roundoff_type1 >= 200) { - error_type = 1; // round off error - integ_dbg2<<"GKP: Round off error 1\n"; + integ_dbg2<<"In adaptive GKP, failed first pass... subdividing\n"; + integ_dbg2<<"Intial range = "< for IntRegoins such that the "largest" is the one + // with the largest current error estimate. This is the next one to be split + // if we need to split further. + std::priority_queue,std::vector > > allregions; + allregions.push(reg); + T finalarea = reg.getArea(); + T finalerr = reg.getErr(); + T tolerance= std::max(abserr, relerr * std::abs(finalarea)); + assert(finalerr > tolerance); + + while(!error_type && finalerr > tolerance) { + // Bisect the subinterval with the largest error estimate + integ_dbg2<<"Current answer = "< parent = allregions.top(); + allregions.pop(); + integ_dbg2<<"Subdividing largest error region "; + integ_dbg2< > children; + parent.subDivide(children); + // For "GKP", there are only two, but for GKPOSC, there is one + // for each oscillation in region + + // Try to do at least 3x better with the children + T factor = 3*children.size()*finalerr/tolerance; + T newabserr = std::abs(parent.getErr()/factor); + T newrelerr = newabserr/std::abs(parent.getArea()); + integ_dbg2<<"New abserr,rel = "<& child = children[i]; + integ_dbg2<<"Integrating child "<= 200.) { - error_type = 2; // round off error - integ_dbg2<<"GKP: Round off error 2\n"; + integ_dbg2<<"Compare: newerr = "<= 10 && newerror > parent.getErr() && + std::abs(delta) <= newerror-parent.getErr()) { + integ_dbg2<<"roundoff type 2: newerror/error = "; + integ_dbg2< tolerance) { + if (roundoff_type1 >= 200) { + error_type = 1; // round off error + integ_dbg2<<"GKP: Round off error 1\n"; + } + if (roundoff_type2 >= 200.) { + error_type = 2; // round off error + integ_dbg2<<"GKP: Round off error 2\n"; + } + const double parent_size = parent.right()-parent.left(); + const double reg_size = reg.right()-parent.left(); + if (std::abs(parent_size / reg_size) < eps) { + error_type = 3; // found singularity + integ_dbg2<<"GKP: Probable singularity\n"; + } + } + for(size_t i=0;i& r=allregions.top(); + finalarea += r.getArea(); + finalerr += r.getErr(); + allregions.pop(); + } + reg.setArea(finalarea,finalerr); + + if (error_type == 1) { + integ_dbg2<<"Type 1 roundoff = "<& r=allregions.top(); - finalarea += r.getArea(); - finalerr += r.getErr(); - allregions.pop(); - } - reg.setArea(finalarea,finalerr); - - if (error_type == 1) { - std::ostringstream s; - s << "Type 1 roundoff's = "< - struct AuxFunc1 : // f(1/x-1) for int(a..infinity) - public std::unary_function - { - public: - AuxFunc1(const UF& _f) : f(_f) {} - typename UF::result_type operator()( - typename UF::argument_type x) const - { return f(1./x-1.)/(x*x); } - private: - const UF& f; - }; - - template - AuxFunc1 inline Aux1(const UF& uf) - { return AuxFunc1(uf); } - template - struct AuxFunc2 : // f(1/x+1) for int(-infinity..b) - public std::unary_function - { - public: - AuxFunc2(const UF& _f) : f(_f) {} - typename UF::result_type operator()( - typename UF::argument_type x) const - { return f(1./x+1.)/(x*x); } - private: - const UF& f; - }; - - template AuxFunc2 - inline Aux2(const UF& uf) - { return AuxFunc2(uf); } - //! @endcond + template + struct AuxFunc1 : // f(1/x-1) for int(a..infinity) + public std::unary_function + { + public: + AuxFunc1(const UF& _f) : f(_f) {} + typename UF::result_type operator()( + typename UF::argument_type x) const + { return f(1./x-1.)/(x*x); } + private: + const UF& f; + }; + + template + AuxFunc1 inline Aux1(const UF& uf) + { return AuxFunc1(uf); } + + template + struct AuxFunc2 : // f(1/x+1) for int(-infinity..b) + public std::unary_function + { + public: + AuxFunc2(const UF& _f) : f(_f) {} + typename UF::result_type operator()( + typename UF::argument_type x) const + { return f(1./x+1.)/(x*x); } + private: + const UF& f; + }; + + template AuxFunc2 + inline Aux2(const UF& uf) + { return AuxFunc2(uf); } + } // anonymous namespace /// Perform a 1-dimensional integral using an IntRegion template @@ -851,36 +835,35 @@ namespace integ { return int1d(func,reg,relerr,abserr); } - // All code between the @cond and @endcond is excluded from Doxygen documentation - //! @cond - template - class Int2DAuxType : - public std::unary_function - { - public: - Int2DAuxType(const BF& _func, const YREG& _yreg, - const typename BF::result_type& _relerr, - const typename BF::result_type& _abserr) : - func(_func),yreg(_yreg),relerr(_relerr),abserr(_abserr) - {} - - typename BF::result_type operator()( - typename BF::first_argument_type x) const + namespace { + template + class Int2DAuxType : + public std::unary_function { - typename YREG::result_type tempreg = yreg(x); - typename BF::result_type result = - int1d(bind21(func,x),tempreg,relerr,abserr); - integ_dbg3<<"Evaluated int2dAux at x = "< @@ -901,37 +884,36 @@ namespace integ { return answer; } - // All code between the @cond and @endcond is excluded from Doxygen documentation - //! @cond - template - class Int3DAuxType : - public std::unary_function - { - public: - Int3DAuxType(const TF& _func, const YREG& _yreg, const ZREG& _zreg, - const typename TF::result_type& _relerr, - const typename TF::result_type& _abserr) : - func(_func),yreg(_yreg),zreg(_zreg),relerr(_relerr),abserr(_abserr) - {} - - typename TF::result_type operator()( - typename TF::firstof3_argument_type x) const + namespace { + template + class Int3DAuxType : + public std::unary_function { - typename YREG::result_type tempreg = yreg(x); - typename TF::result_type result = - int2d(bind31(func,x),tempreg,bind21(zreg,x),relerr,abserr); - integ_dbg3<<"Evaluated int3dAux at x = "< @@ -956,28 +938,27 @@ namespace integ { // Helpers for constant regions for int2d, int3d: - // All code between the @cond and @endcond is excluded from Doxygen documentation - //! @cond - template - struct ConstantReg1 : - public std::unary_function > - { - ConstantReg1(T a,T b) : ir(a,b) {} - ConstantReg1(const IntRegion& r) : ir(r) {} - IntRegion operator()(T x) const { return ir; } - IntRegion ir; - }; - - template - struct ConstantReg2 : - public std::binary_function > - { - ConstantReg2(T a,T b) : ir(a,b) {} - ConstantReg2(const IntRegion& r) : ir(r) {} - IntRegion operator()(T x, T y) const { return ir; } - IntRegion ir; - }; - //! @endcond + namespace { + template + struct ConstantReg1 : + public std::unary_function > + { + ConstantReg1(T a,T b) : ir(a,b) {} + ConstantReg1(const IntRegion& r) : ir(r) {} + IntRegion operator()(T x) const { return ir; } + IntRegion ir; + }; + + template + struct ConstantReg2 : + public std::binary_function > + { + ConstantReg2(T a,T b) : ir(a,b) {} + ConstantReg2(const IntRegion& r) : ir(r) {} + IntRegion operator()(T x, T y) const { return ir; } + IntRegion ir; + }; + } // anonymous namespace /// Perform a 3-dimensional integral using constant IntRegions for both regions /// (i.e. the integral is over a square) diff --git a/pysrc/Integ.cpp b/pysrc/Integ.cpp new file mode 100644 index 00000000000..53d8233e36e --- /dev/null +++ b/pysrc/Integ.cpp @@ -0,0 +1,51 @@ +#include "boost/python.hpp" +#include "integ/Int.h" + +#include + +namespace bp = boost::python; + +namespace galsim { +namespace integ { +namespace { + + // A C++ function object that just calls a python function. + class PyFunc : + public std::unary_function + { + public: + PyFunc(const bp::object& func) : _func(func) {} + double operator()(double x) const + { return bp::extract(_func(x)); } + private: + const bp::object& _func; + }; + + // Integrate a python function using int1d. + bp::tuple PyInt1d(const bp::object& func, double min, double max, + double rel_err=DEFRELERR, double abs_err=DEFABSERR) + { + PyFunc pyfunc(func); + try { + double res = int1d(pyfunc, min, max, rel_err, abs_err); + return bp::make_tuple(true, res); + } catch (std::exception& e) { + return bp::make_tuple(false, e.what()); + } + } + +} // anonymous + + +void pyExportInteg() { + + bp::def("PyInt1d", + &PyInt1d, (bp::args("func", "min", "max"), + bp::arg("rel_err")=DEFRELERR, bp::arg("abs_err")=DEFABSERR), + "Calculate the integral of the given 1-d function from min to max."); + +} + +} // namespace integ +} // namespace galsim + diff --git a/pysrc/files.txt b/pysrc/files.txt index 5b9e5df6a79..fe946d26895 100644 --- a/pysrc/files.txt +++ b/pysrc/files.txt @@ -20,3 +20,4 @@ SBKolmogorov.cpp PhotonArray.cpp Random.cpp PSFCorr.cpp +Integ.cpp diff --git a/pysrc/module.cpp b/pysrc/module.cpp index 5ea813ddf7e..88a1a4ce7a8 100644 --- a/pysrc/module.cpp +++ b/pysrc/module.cpp @@ -26,10 +26,14 @@ namespace galsim { void pyExportSBKolmogorov(); void pyExportRandom(); - namespace hsm{ + namespace hsm { void pyExportPSFCorr(); } // namespace hsm + namespace integ { + void pyExportInteg(); + } // namespace hsm + } // namespace galsim BOOST_PYTHON_MODULE(_galsim) { @@ -55,4 +59,5 @@ BOOST_PYTHON_MODULE(_galsim) { galsim::pyExportSBKolmogorov(); galsim::pyExportRandom(); galsim::hsm::pyExportPSFCorr(); + galsim::integ::pyExportInteg(); } diff --git a/tests/files.txt b/tests/files.txt index f799ba6710b..d77dc263b98 100644 --- a/tests/files.txt +++ b/tests/files.txt @@ -1 +1,2 @@ test_Image.cpp +test_integ.cpp diff --git a/tests/run_all_tests b/tests/run_all_tests index f69ac5795c3..4dabdf340de 100755 --- a/tests/run_all_tests +++ b/tests/run_all_tests @@ -13,3 +13,4 @@ python test_shear.py || exit python test_spline.py || exit python test_utilities.py || exit python test_lensing.py || exit +python test_integ.py || exit diff --git a/tests/test_integ.cpp b/tests/test_integ.cpp new file mode 100644 index 00000000000..9811c2915c1 --- /dev/null +++ b/tests/test_integ.cpp @@ -0,0 +1,186 @@ + +#include "galsim/integ/Int.h" + +#define BOOST_TEST_DYN_LINK + +// icpc pretends to be GNUC, since it thinks it's compliant, but it's not. +// It doesn't understand "pragma GCC" +#ifndef __INTEL_COMPILER + +// The boost unit tests have some unused variables, so suppress the warnings about that. +// I think pragma GCC was introduced in gcc 4.2, so guard for >= that version +#if defined(__GNUC__) && __GNUC__ >= 4 && (__GNUC__ >= 5 || __GNUC_MINOR__ >= 2) +#pragma GCC diagnostic ignored "-Wunused-variable" +#endif + +// Not sure when this was added. Currently check for it for versions >= 4.3 +#if defined(__GNUC__) && __GNUC__ >= 4 && (__GNUC__ >= 5 || __GNUC_MINOR__ >= 3) +#pragma GCC diagnostic ignored "-Warray-bounds" +#endif + +#endif + +#include +#include +#include + +const double test_sigma = 7.; // test value of Gaussian sigma for integral tests +const double test_rel_err = 1.e-7; // the relative accuracy at which to test +const double test_abs_err = 1.e-13; // the absolute accuracy at which to test +const double test_mock_inf = 2.e10; // number large enough to get interpreted as infinity by + // integration routines +// Note: all "true" answers in the below tests are found using Wolfram Alpha. + +// A simple Gaussian that works as a functional object +class Gauss : public std::unary_function +{ +public : + Gauss(double sig) : _sig(sig) {} + + double operator()(double x) const + { return exp(-0.5*pow(x/_sig,2)); } + +private : + double _sig; +}; + +// A simple power law +class Power : public std::unary_function +{ +public : + Power(double expon) : _expon(expon) {} + + double operator()(double x) const + { return pow(x,_expon); } + +private : + double _expon; +}; + +// A straight function, rather than a functional class: +double osc_func(double x) +{ return sin(pow(x,2)) * exp(-std::abs(x)); } + +// A simple function: +// f(x,y) = x*(3*x+y) + y +double twod_func(double x, double y) +{ return x * (3.*x + y) + y; } + +BOOST_AUTO_TEST_SUITE(integ_tests); + +BOOST_AUTO_TEST_CASE( TestGaussian ) +{ + Gauss gauss(test_sigma); + + BOOST_CHECK_CLOSE( + galsim::integ::int1d(gauss, -1., 1., test_rel_err, test_abs_err), + 1.99321805307377285009, + // Note: it seems like this should be test_rel_err, but boost has this last + // parameter as a _percentage_ error. (!?!) So need to multiply by 100. + 100 * test_rel_err); + + BOOST_CHECK_CLOSE( + galsim::integ::int1d(gauss, 0., 20., test_rel_err, test_abs_err), + 8.73569586966967345835, + 100 * test_rel_err); + + BOOST_CHECK_CLOSE( + galsim::integ::int1d(gauss, -50., -40., test_rel_err, test_abs_err), + 9.66426031085587421984e-8, + // Difference is allowed to be test_abs_err, which corresponds to a relative error of + // test_abs_err / 9.66e-8 + 100 * test_abs_err / 9.66e-8); + + BOOST_CHECK_CLOSE( + galsim::integ::int1d(gauss, 0., test_mock_inf, test_rel_err, test_abs_err), + 8.77319896120850210849, + 100 * test_rel_err); + + BOOST_CHECK_CLOSE( + galsim::integ::int1d(gauss, -test_mock_inf, 5.4, test_rel_err, test_abs_err), + 13.68221660030048620971, + 100 * test_rel_err); + + BOOST_CHECK_CLOSE( + galsim::integ::int1d(gauss, -test_mock_inf, test_mock_inf, test_rel_err, test_abs_err), + 17.54639792241700421699, + 100 * test_rel_err); +} + +BOOST_AUTO_TEST_CASE( TestOscillatory ) +{ + BOOST_CHECK_CLOSE( + galsim::integ::int1d(std::ptr_fun(osc_func), -1., 1., test_rel_err, test_abs_err), + 0.30182513444548879567, + 100 * test_rel_err); + + BOOST_CHECK_CLOSE( + galsim::integ::int1d(std::ptr_fun(osc_func), 0., 20., test_rel_err, test_abs_err), + 0.27051358019041255485, + 100 * test_rel_err); + + BOOST_CHECK_CLOSE( + galsim::integ::int1d(std::ptr_fun(osc_func), -15., -14., test_rel_err, test_abs_err), + 7.81648378350593176887e-9, + 100 * test_abs_err / 7.82e-8); + + BOOST_CHECK_CLOSE( + galsim::integ::int1d(std::ptr_fun(osc_func), 0., test_mock_inf, test_rel_err, + test_abs_err), + 0.27051358016221414426, + 100 * test_rel_err); + + BOOST_CHECK_CLOSE( + galsim::integ::int1d(std::ptr_fun(osc_func), -test_mock_inf, 5.4, test_rel_err, + test_abs_err), + 0.5413229824941895221, + 100 * test_rel_err); + + BOOST_CHECK_CLOSE( + galsim::integ::int1d(std::ptr_fun(osc_func), -test_mock_inf, test_mock_inf, test_rel_err, + test_abs_err), + 0.54102716032442828852, + 100 * test_rel_err); +} + +BOOST_AUTO_TEST_CASE( TestPole ) +{ + Power powm05(-0.5); + + BOOST_CHECK_CLOSE( + galsim::integ::int1d(powm05, 0., 1., test_rel_err, test_abs_err), + 2., + 100 * test_rel_err); + + BOOST_CHECK_CLOSE( + galsim::integ::int1d(powm05, 0., 300., test_rel_err, test_abs_err), + 34.64101615137754587055, + 100 * test_rel_err); + + Power powm2(-2.); + + BOOST_CHECK_CLOSE( + galsim::integ::int1d(powm2, 1., 2., test_rel_err, test_abs_err), + 0.5, + 100 * test_rel_err); + + BOOST_CHECK_CLOSE( + galsim::integ::int1d(powm2, 1., test_mock_inf, test_rel_err, test_abs_err), + 1., + 100 * test_rel_err); + + BOOST_CHECK_THROW( + galsim::integ::int1d(powm2, 0., 1., test_rel_err, test_abs_err), + galsim::integ::IntFailure); +} + +BOOST_AUTO_TEST_CASE( Test2d ) +{ + BOOST_CHECK_CLOSE( + galsim::integ::int2d(std::ptr_fun(twod_func),0.,1.,0.,1., test_rel_err, test_abs_err), + 1.75, + 100 * test_rel_err); +} + +BOOST_AUTO_TEST_SUITE_END(); + diff --git a/tests/test_integ.py b/tests/test_integ.py new file mode 100644 index 00000000000..61ddd9e6900 --- /dev/null +++ b/tests/test_integ.py @@ -0,0 +1,215 @@ +"""Unit tests for integration routines at the Python layer. +""" + +import numpy as np +try: + import galsim +except ImportError: + import os + import sys + path, filename = os.path.split(__file__) + sys.path.append(os.path.abspath(os.path.join(path, ".."))) + import galsim + +test_sigma = 7. # test value of Gaussian sigma for integral tests +test_rel_err = 1.e-7 # the relative accuracy at which to test +test_abs_err = 1.e-13 # the absolute accuracy at which to test +test_mock_inf = 2.e10 # number large enough to get interpreted as infinity by + # integration routines +test_decimal = 7 + +def funcname(): + import inspect + return inspect.stack()[1][3] + +def test_gaussian_finite_limits(): + """Test the integration of a 1D zero-mean Gaussian across intervals of [-1, 1], [0, 20] + and [-50, -40]. + """ + import time + t1 = time.time() + + # Define our test function + def test_func(x): return np.exp(-.5 * x**2 / test_sigma**2) + + test_integral = galsim.integ.int1d(test_func, -1., 1., test_rel_err, test_abs_err) + # test results easily calculated using Wolfram alpha + true_result = 1.99321805307377285009 + np.testing.assert_almost_equal( + test_integral, true_result, decimal=test_decimal, verbose=True, + err_msg="Gaussian integral failed across interval [-1, 1].") + + test_integral = galsim.integ.int1d(test_func, 0., 20., test_rel_err, test_abs_err) + true_result = 8.73569586966967345835 + np.testing.assert_almost_equal( + test_integral, true_result, decimal=test_decimal, verbose=True, + err_msg="Gaussian integral failed across interval [0, 20].") + + test_integral = galsim.integ.int1d(test_func, -50., -40., test_rel_err, test_abs_err) + true_result = 9.66426031085587421984e-8 + np.testing.assert_almost_equal( + test_integral, true_result, decimal=test_decimal, verbose=True, + err_msg="Gaussian integral failed across interval [-50, -40].") + + t2 = time.time() + print 'time for %s = %.2f'%(funcname(),t2-t1) + +def test_gaussian_infinite_limits(): + """Test the integration of a 1D zero-mean Gaussian across intervals of [0, inf], [-inf, 5.4] + and [-inf, inf]. + """ + import time + t1 = time.time() + + # Define our test function + def test_func(x): return np.exp(-.5 * x**2 / test_sigma**2) + + test_integral = galsim.integ.int1d(test_func, 0., test_mock_inf, test_rel_err, test_abs_err) + # test results easily calculated using Wolfram alpha + true_result = 8.77319896120850210849 + np.testing.assert_almost_equal( + test_integral, true_result, decimal=test_decimal, verbose=True, + err_msg="Gaussian integral failed across interval [0, inf].") + + test_integral = galsim.integ.int1d(test_func, -test_mock_inf, 5.4, test_rel_err, test_abs_err) + true_result = 13.68221660030048620971 + np.testing.assert_almost_equal( + test_integral, true_result, decimal=test_decimal, verbose=True, + err_msg="Gaussian integral failed across interval [-inf, 5.4].") + + test_integral = galsim.integ.int1d( + test_func, -test_mock_inf, test_mock_inf, test_rel_err, test_abs_err) + true_result = 17.54639792241700421699 + np.testing.assert_almost_equal( + test_integral, true_result, decimal=test_decimal, verbose=True, + err_msg="Gaussian integral failed across interval [-inf, inf].") + + t2 = time.time() + print 'time for %s = %.2f'%(funcname(),t2-t1) + +def test_sinxsqexpabsx_finite_limits(): + """Test the integration of a slightly tricky oscillating sin(x^2) * exp(-|x|) function across + finite intervals [-1, 1], [0, 20], [-15, 14]. + """ + import time + t1 = time.time() + + # Define our test function + def test_func(x): return np.sin(x**2) * np.exp(-np.abs(x)) + + test_integral = galsim.integ.int1d(test_func, -1., 1., test_rel_err, test_abs_err) + # test results easily calculated using Wolfram alpha + true_result = 0.30182513444548879567 + np.testing.assert_almost_equal( + test_integral, true_result, decimal=test_decimal, verbose=True, + err_msg="Sin(x^2) * exp(-|x|) integral failed across interval [-1, 1].") + + test_integral = galsim.integ.int1d(test_func, 0., 20., test_rel_err, test_abs_err) + true_result = 0.27051358019041255485 + np.testing.assert_almost_equal( + test_integral, true_result, decimal=test_decimal, verbose=True, + err_msg="Sin(x^2) * exp(-|x|) integral failed across interval [0, 20].") + + test_integral = galsim.integ.int1d(test_func, -15., -14., test_rel_err, test_abs_err) + true_result = 7.81648378350593176887e-9 + np.testing.assert_almost_equal( + (test_integral - true_result) / true_result, 0., decimal=test_decimal, verbose=True, + err_msg="Sin(x^2) * exp(-|x|) integral failed across interval [-15, -14].") + + t2 = time.time() + print 'time for %s = %.2f'%(funcname(),t2-t1) + +def test_sinxsqexpabsx_infinite_limits(): + """Test the integration of a slightly tricky oscillating sin(x^2) * exp(-|x|) function across + infinite intervals [0, inf], [-inf, 5.4], [-inf, inf]. + """ + import time + t1 = time.time() + + # Define our test function + def test_func(x): return np.sin(x**2) * np.exp(-np.abs(x)) + + test_integral = galsim.integ.int1d(test_func, 0., test_mock_inf, test_rel_err, test_abs_err) + # test results easily calculated using Wolfram alpha + true_result = 0.27051358016221414426 + np.testing.assert_almost_equal( + test_integral, true_result, decimal=test_decimal, verbose=True, + err_msg="Sin(x^2) * exp(-|x|) integral failed across interval [0, inf].") + + test_integral = galsim.integ.int1d(test_func, -test_mock_inf, 5.4, test_rel_err, test_abs_err) + true_result = 0.5413229824941895221 + np.testing.assert_almost_equal( + test_integral, true_result, decimal=test_decimal, verbose=True, + err_msg="Sin(x^2) * exp(-|x|) integral failed across interval [-inf, 5.4].") + + test_integral = galsim.integ.int1d( + test_func, -test_mock_inf, test_mock_inf, test_rel_err, test_abs_err) + true_result = 0.54102716032442828852 + np.testing.assert_almost_equal( + test_integral, true_result, decimal=test_decimal, verbose=True, + err_msg="Sin(x^2) * exp(-|x|) integral failed across interval [-inf, inf].") + + t2 = time.time() + print 'time for %s = %.2f'%(funcname(),t2-t1) + +def test_invroot_finite_limits(): + """Test the integration of |x|^(-1/2) across intervals [0,1], [0,300] (integrable pole at x=0). + """ + import time + t1 = time.time() + + # Define our test function + def test_func(x): return 1. / np.sqrt(np.abs(x)) + test_integral = galsim.integ.int1d(test_func, 0, 1., test_rel_err, test_abs_err) + # test results easily calculated using Wolfram alpha + true_result = 2. + np.testing.assert_almost_equal( + test_integral, true_result, decimal=test_decimal, verbose=True, + err_msg="|x|^(-1/2) integral failed across interval [0, 1].") + + test_integral = galsim.integ.int1d(test_func, 0., 300., test_rel_err, test_abs_err) + true_result = 34.64101615137754587055 + np.testing.assert_almost_equal( + test_integral, true_result, decimal=test_decimal, verbose=True, + err_msg="|x|^(-1/2) integral failed across interval [0, 300].") + + t2 = time.time() + print 'time for %s = %.2f'%(funcname(),t2-t1) + +def test_invroot_infinite_limits(): + """Test the integration of |x|^(-2) across intervals [1,2], [1,inf]. + Also check that [0,1] raises an exception. + """ + import time + t1 = time.time() + + # Define our test function + def test_func(x): return x**-2 + test_integral = galsim.integ.int1d(test_func, 1., 2., test_rel_err, test_abs_err) + true_result = 0.5 + np.testing.assert_almost_equal( + test_integral, true_result, decimal=test_decimal, verbose=True, + err_msg="x^(-2) integral failed across interval [1, 2].") + + test_integral = galsim.integ.int1d(test_func, 1., test_mock_inf, test_rel_err, test_abs_err) + true_result = 1.0 + np.testing.assert_almost_equal( + test_integral, true_result, decimal=test_decimal, verbose=True, + err_msg="x^(-2) integral failed across interval [1, inf].") + + np.testing.assert_raises( + RuntimeError, + galsim.integ.int1d, test_func, 0., 1., test_rel_err, test_abs_err) + + t2 = time.time() + print 'time for %s = %.2f'%(funcname(),t2-t1) + + +if __name__ == "__main__": + test_gaussian_finite_limits() + test_gaussian_infinite_limits() + test_sinxsqexpabsx_finite_limits() + test_sinxsqexpabsx_infinite_limits() + test_invroot_finite_limits() + test_invroot_infinite_limits() + diff --git a/tests/test_lensing.py b/tests/test_lensing.py index e2585f5dc96..62fcdb53f0a 100644 --- a/tests/test_lensing.py +++ b/tests/test_lensing.py @@ -28,49 +28,45 @@ def test_nfwhalo(): # distance go from 1 .. 599 arcsec ref = np.loadtxt(refdir + '/nfw_lens.dat') + # set up the same halo + halo = galsim.lensing.NFWHalo(mass=1e15, conc=4, z=1, pos_x=0, pos_y=0) + pos_x = np.arange(1,600) + pos_y = np.zeros_like(pos_x) + z_s = 2 + kappa = halo.getConvergence(pos_x, pos_y, z_s) + gamma1, gamma2 = halo.getShear(pos_x, pos_y, z_s, reduced=False) + g1, g2 = halo.getShear(pos_x, pos_y, z_s, reduced=True) + + # check internal correctness: + # g1 = gamma1/(1-kappa), and g2 = 0 + np.testing.assert_array_equal(g1, gamma1/(1-kappa), + err_msg="Computation of reduced shear g incorrect.") + np.testing.assert_array_equal(g2, np.zeros_like(g2), + err_msg="Computation of reduced shear g2 incorrect.") + + # comparison to reference: + # tangential shear in x-direction is purely negative in g1 try: - # set up the same halo - halo = galsim.lensing.NFWHalo(mass=1e15, conc=4, z=1, pos_x=0, pos_y=0) - pos_x = np.arange(1,600) - pos_y = np.zeros_like(pos_x) - z_s = 2 - kappa = halo.getConvergence(pos_x, pos_y, z_s) - gamma1, gamma2 = halo.getShear(pos_x, pos_y, z_s, reduced=False) - g1, g2 = halo.getShear(pos_x, pos_y, z_s, reduced=True) - - # check internal correctness: - # g1 = gamma1/(1-kappa), and g2 = 0 - np.testing.assert_array_equal(g1, gamma1/(1-kappa), - err_msg="Computation of reduced shear g incorrect.") - np.testing.assert_array_equal(g2, np.zeros_like(g2), - err_msg="Computation of reduced shear g2 incorrect.") - - # comparison to reference: - # tangential shear in x-direction is purely negative in g1 - try: - np.testing.assert_allclose( - -ref[:,2], gamma1, rtol=1e-4, - err_msg="Computation of shear deviates from reference.") - np.testing.assert_allclose( - -ref[:,3], g1, rtol=1e-4, - err_msg="Computation of reduced shear deviates from reference.") - np.testing.assert_allclose( - ref[:,4], kappa, rtol=1e-4, - err_msg="Computation of convergence deviates from reference.") - except AttributeError: - # Older numpy versions don't have assert_allclose, so use this instead: - np.testing.assert_array_almost_equal( - -ref[:,2], gamma1, decimal=4, - err_msg="Computation of shear deviates from reference.") - np.testing.assert_array_almost_equal( - -ref[:,3], g1, decimal=4, - err_msg="Computation of reduced shear deviates from reference.") - np.testing.assert_array_almost_equal( - ref[:,4], kappa, decimal=4, - err_msg="Computation of convergence deviates from reference.") - except ImportError: - import warnings - warnings.warn("scipy not found! Integrator required for angular diameter distances") + np.testing.assert_allclose( + -ref[:,2], gamma1, rtol=1e-4, + err_msg="Computation of shear deviates from reference.") + np.testing.assert_allclose( + -ref[:,3], g1, rtol=1e-4, + err_msg="Computation of reduced shear deviates from reference.") + np.testing.assert_allclose( + ref[:,4], kappa, rtol=1e-4, + err_msg="Computation of convergence deviates from reference.") + except AttributeError: + # Older numpy versions don't have assert_allclose, so use this instead: + np.testing.assert_array_almost_equal( + -ref[:,2], gamma1, decimal=4, + err_msg="Computation of shear deviates from reference.") + np.testing.assert_array_almost_equal( + -ref[:,3], g1, decimal=4, + err_msg="Computation of reduced shear deviates from reference.") + np.testing.assert_array_almost_equal( + ref[:,4], kappa, decimal=4, + err_msg="Computation of convergence deviates from reference.") t2 = time.time() print 'time for %s = %.2f'%(funcname(),t2-t1)