Skip to content

Commit

Permalink
fix compilation
Browse files Browse the repository at this point in the history
  • Loading branch information
Matthias Wittgen committed Apr 1, 2021
1 parent 04f8c2c commit c162694
Show file tree
Hide file tree
Showing 2 changed files with 48 additions and 132 deletions.
150 changes: 19 additions & 131 deletions include/lsst/afw/math/Integrate.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -25,7 +25,7 @@
#if !defined(LSST_AFW_MATH_INTEGRATE_H)
#define LSST_AFW_MATH_INTEGRATE_H 1
/*
* Compute 1d and 2d integral
* Compute 1d
*/

#include <functional>
Expand All @@ -45,30 +45,24 @@

#include "lsst/afw/math/IntGKPData10.h"

// == The following is Mike Jarvis original comment ==
// == The following is based on Mike Jarvis original comment ==
//
// Basic Usage:
//
// First, define a function object, which should derive from
// std::unary_function<double,double>. For example, to integrate a
// First, define a function object.
// For example, to integrate a
// Gaussian, use something along the lines of this:
//
// class Gauss :
// public std::unary_function<double,double>
// {
// public :
//
// Gauss(double _mu, double _sig) :
// mu(_mu), sig(_sig), sigsq(_sig*_sig) {}
//
// double operator()(double x) const
// {
// const double SQRTTWOPI = 2.50662827463;
// return exp(-pow(x-mu,2)/2./sigsq)/SQRTTWOPI/sig;
// }
//
// private :
// double mu,sig,sigsq;
// class Gauss {
// public :
// Gauss(double _mu, double _sig) : mu(_mu), sig(_sig) {}
// double operator()(double x) const {
// constexpr double inv_sqrt_2pi = 0.3989422804014327;
// double a = (x - mu) / sig;
// return inv_sqrt_2pi / sig * std::exp(-0.5f * a * a);
// }
// private :
// double mu,sig;
// };
//
// Next, make an IntRegion object with the bounds of the integration region.
Expand All @@ -78,19 +72,15 @@
//
// For example, to integrate something from -1 to 1, use:
//
// integ::IntRegion<double> reg1(-1.,1.);
//
// (The integ:: is because everything here is in the integ namespace to
// help prevent name conflicts with other header files.)
//
// lsst::afw::math::IntRegion<double> reg1(-1.,1.);
// If a value is > 1.e10 or < -1.e10, then these values are taken to be
// infinity, rather than the actual value.
// So to integrate from 0 to infinity:
//
// integ::IntRegion<double> reg2(0.,1.e100);
//
// Or, you can use the variable integ::MOCK_INF which might be clearer.
// lsst::afw::math::IntRegion<double> reg2(0.,20);
//
// Or, you can use the variable lsst::afw::math:MOCK_INF which might be clearer.
// The integral might diverge depending on the limits chosen
//
// Finally, to perform the integral, the line would be:
//
Expand Down Expand Up @@ -138,31 +128,7 @@
//
// (It is intended to be an overestimate of the actual error,
// but it doesn't always get it completely right.)
//
//
//
// Two- and Three-Dimensional Integrals:
//
// These are slightly more complicated. The easiest case is when the
// bounds of the integral are a rectangle or 3d box. In this case,
// you can still use the regular IntRegion. The only new thing then
// is the definition of the function. For example, to integrate
// int(3x^2 + xy + y , x=0..1, y=0..1):
//
// struct Integrand :
// public std::binary_function<double,double,double>
// {
// double operator()(double x, double y) const
// { return x*(3.*x + y) + y; }
// };
//
// integ::IntRegion<double> reg3(0.,1.);
// double integ3 = int2d(Integrand(),reg3,reg3);
//
// (Which should give 1.75 as the result.)
//
//
//


namespace lsst {
namespace afw {
Expand Down Expand Up @@ -643,25 +609,6 @@ AuxFunc2<UF> inline Aux2(UF uf) {
return AuxFunc2<UF>(uf);
}

/**
* Helpers for constant regions for int2d, int3d:
*
*/
template <class T>
struct ConstantReg1 : public std::unary_function<T, IntRegion<T> > {
ConstantReg1(T a, T b) : ir(a, b) {}
ConstantReg1(IntRegion<T> const &r) : ir(r) {}
IntRegion<T> operator()(T) const { return ir; }
IntRegion<T> ir;
};

template <class T>
struct ConstantReg2 : public std::binary_function<T, T, IntRegion<T> > {
ConstantReg2(T a, T b) : ir(a, b) {}
ConstantReg2(IntRegion<T> const &r) : ir(r) {}
IntRegion<T> operator()(T x, T y) const { return ir; }
IntRegion<T> ir;
};

// pulled from MoreFunctional.h. Needed in class Int2DAuxType and Int3DAuxType
template <class BF>
Expand Down Expand Up @@ -808,65 +755,6 @@ inline Arg int1d(UnaryFunctionT func, IntRegion<Arg> &reg,
}
}

/**
* Front end for the 2d integrator
*/
template <class BF, class YREG>
inline typename BF::result_type int2d(BF const &func, IntRegion<typename BF::result_type> &reg,
YREG const &yreg, typename BF::result_type const &abserr = DEFABSERR,
typename BF::result_type const &relerr = DEFRELERR) {
using namespace details;
integ_dbg2 << "Starting int2d: range = ";
integ_dbg2 << reg.Left() << ".." << reg.Right() << std::endl;
Int2DAuxType<BF, YREG> faux(func, yreg, abserr * 1.0e-3, relerr * 1.0e-3);
typename BF::result_type answer = int1d(faux, reg, abserr, relerr);
integ_dbg2 << "done int2d answer = " << answer << " +- " << reg.Err() << std::endl;
return answer;
}

/**
* Front end for the 3d integrator
*/
template <class TF, class YREG, class ZREG>
inline typename TF::result_type int3d(TF const &func, IntRegion<typename TF::result_type> &reg,
YREG const &yreg, ZREG const &zreg,
typename TF::result_type const &abserr = DEFABSERR,
typename TF::result_type const &relerr = DEFRELERR) {
using namespace details;
integ_dbg2 << "Starting int3d: range = ";
integ_dbg2 << reg.Left() << ".." << reg.Right() << std::endl;
Int3DAuxType<TF, YREG, ZREG> faux(func, yreg, zreg, abserr * 1.e-3, relerr * 1.e-3);
typename TF::result_type answer = int1d(faux, reg, abserr, relerr);
integ_dbg2 << "done int3d answer = " << answer << " +- " << reg.Err() << std::endl;
return answer;
}

/**
* Front end for the 2d integrator
*/
template <class BF>
inline typename BF::result_type int2d(BF const &func, IntRegion<typename BF::result_type> &reg,
IntRegion<typename BF::result_type> &yreg,
typename BF::result_type const &abserr = DEFABSERR,
typename BF::result_type const &relerr = DEFRELERR) {
using namespace details;
return int2d(func, reg, ConstantReg1<typename BF::result_type>(yreg), abserr, relerr);
}

/**
* Front end for the 3d integrator
*/
template <class TF>
inline typename TF::result_type int3d(TF const &func, IntRegion<typename TF::result_type> &reg,
IntRegion<typename TF::result_type> &yreg,
IntRegion<typename TF::result_type> &zreg,
typename TF::result_type const &abserr = DEFABSERR,
typename TF::result_type const &relerr = DEFRELERR) {
using namespace details;
return int3d(func, reg, ConstantReg1<typename TF::result_type>(yreg),
ConstantReg2<typename TF::result_type>(zreg), abserr, relerr);
}

// =============================================================
/**
* The 1D integrator
Expand Down
30 changes: 29 additions & 1 deletion tests/integrate.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -129,7 +129,7 @@ BOOST_AUTO_TEST_CASE(
double parab_area_integrate = math::integrate(parab1d, x1, x2);
double parab_area_analytic = parab1d.getAnalyticArea(x1, x2);

double parab_area_integrate_function = math::integrate(std::function(parabola1d), x1, x2);
double parab_area_integrate_function = math::integrate(parabola1d, x1, x2);

BOOST_CHECK_CLOSE(parab_area_integrate, parab_area_analytic, 1e-6);
BOOST_CHECK_CLOSE(parab_area_integrate_function, parab_area_analytic, 1e-6);
Expand Down Expand Up @@ -158,3 +158,31 @@ BOOST_AUTO_TEST_CASE(
BOOST_CHECK_CLOSE(parab_volume_integrate, parab_volume_analytic, 1e-6);
BOOST_CHECK_CLOSE(parab_volume_integrate_function, parab_volume_analytic, 1e-6);
}

/*
* Test int1d function on Gaussian
*/
BOOST_AUTO_TEST_CASE(Integrate) {
constexpr double tolerance = 1e-6;
constexpr double result1 = 0.68268949;
constexpr double result2 = 0.5;
class Gauss {
public :
Gauss(double _mu, double _sig) : mu(_mu), sig(_sig) {}

double operator()(double x) const {
constexpr double inv_sqrt_2pi = 0.3989422804014327;
double a = (x - mu) / sig;
return inv_sqrt_2pi / sig * std::exp(-0.5d * a * a);
}

private :
double mu, sig;
};
lsst::afw::math::IntRegion<double> reg1(-1., 1.);
lsst::afw::math::IntRegion<double> reg2(0, 20);
double integ1 = lsst::afw::math::int1d(Gauss(0., 1.), reg1, 1.e-10, 1.e-4);
double integ2 = lsst::afw::math::int1d(Gauss(0., 2.), reg2, 1.e-10, 1.e-4);
BOOST_CHECK_CLOSE(integ1, result1, tolerance);
BOOST_CHECK_CLOSE(integ2, result2, tolerance);
}

0 comments on commit c162694

Please sign in to comment.