Mitsuba Renderer  0.5.0
Go to the documentation of this file.
1 /*
2  This file is part of Mitsuba, a physically based rendering system.
3
4  Copyright (c) 2007-2014 by Wenzel Jakob and others.
5
6  Mitsuba is free software; you can redistribute it and/or modify
7  it under the terms of the GNU General Public License Version 3
8  as published by the Free Software Foundation.
9
10  Mitsuba is distributed in the hope that it will be useful,
11  but WITHOUT ANY WARRANTY; without even the implied warranty of
12  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13  GNU General Public License for more details.
14
15  You should have received a copy of the GNU General Public License
16  along with this program. If not, see <http://www.gnu.org/licenses/>.
17 */
18
19 #pragma once
22
23 #include <mitsuba/mitsuba.h>
24 #include <boost/function.hpp>
25
27
28 // -----------------------------------------------------------------------
29 //! \ingroup libcore
30 //! \ingroup libpython
31 //! @{ \name Basic tools for non-adaptive quadrature
32 // -----------------------------------------------------------------------
33
34 /// Evaluate the l-th Legendre polynomial using recurrence (single precision)
35 extern MTS_EXPORT_CORE float legendreP(int l, float x);
36
37 /// Evaluate the l-th Legendre polynomial using recurrence (double precision)
38 extern MTS_EXPORT_CORE double legendreP(int l, double x);
39
40 /// Evaluate an associated Legendre polynomial using recurrence (single precision)
41 extern MTS_EXPORT_CORE float legendreP(int l, int m, float x);
42
43 /// Evaluate an associated Legendre polynomial using recurrence (double precision)
44 extern MTS_EXPORT_CORE double legendreP(int l, int m, double x);
45
46 /// Evaluate the l-th Legendre polynomial and its derivative using recurrence (single precision)
47 extern MTS_EXPORT_CORE std::pair<float, float> legendrePD(int l, float x);
48
49 /// Evaluate the l-th Legendre polynomial and its derivative using recurrence (double precision)
50 extern MTS_EXPORT_CORE std::pair<double, double> legendrePD(int l, double x);
51
52 /**
53  * \brief Computes the nodes and weights of a Gauss-Legendre quadrature
54  * (aka "Gaussian quadrature") rule with the given number of evaluations.
55  *
56  * Integration is over the interval \f$[-1, 1]\f$. Gauss-Legendre quadrature
57  * maximizes the order of exactly integrable polynomials achieves this up to
58  * degree \f$2n-1\f$ (where \f$n\f$ is the number of function evaluations).
59  *
60  * This method is numerically well-behaved until about \f$n=200\f$
61  * and then becomes progressively less accurate. It is generally not a
62  * good idea to go much higher---in any case, a composite or
63  * adaptive integration scheme will be superior for large \f$n\f$.
64  *
65  * \param n
66  * Desired number of evalution points
67  * \param nodes
68  * Length-\c n array used to store the nodes of the quadrature rule
69  * \param nodes
70  * Length-\c n array used to store the weights of the quadrature rule
71  */
72 extern MTS_EXPORT_CORE void gaussLegendre(int n, Float *nodes, Float *weights);
73
74
75 /**
76  * \brief Computes the nodes and weights of a Gauss-Lobatto quadrature
77  * rule with the given number of evaluations.
78  *
79  * Integration is over the interval \f$[-1, 1]\f$. Gauss-Lobatto quadrature
80  * is preferable to Gauss-Legendre quadrature whenever the endpoints of the
81  * integration domain should explicitly be included. It maximizes the order
82  * of exactly integrable polynomials subject to this constraint and achieves
83  * this up to degree \f$2n-3\f$ (where \f$n\f$ is the number of function
84  * evaluations).
85  *
86  * This method is numerically well-behaved until about \f$n=200\f$
87  * and then becomes progressively less accurate. It is generally not a
88  * good idea to go much higher---in any case, a composite or
89  * adaptive integration scheme will be superior for large \f$n\f$.
90  *
91  * \param n
92  * Desired number of evalution points
93  * \param nodes
94  * Length-\c n array used to store the nodes of the quadrature rule
95  * \param nodes
96  * Length-\c n array used to store the weights of the quadrature rule
97  */
98 extern MTS_EXPORT_CORE void gaussLobatto(int n, Float *nodes, Float *weights);
99
100 //! @}
101 // -----------------------------------------------------------------------
102
103 // -----------------------------------------------------------------------
104 //! @{ \name Adaptive quadrature schemes
105 // -----------------------------------------------------------------------
106
107 /**
108  * \brief Computes the integral of a one-dimensional function
110  *
111  * Given a target error \f$\epsilon \f$, the integral of
112  * a function \f$f \f$ between \f$a \f$ and \f$b \f$ is
113  * calculated by means of the Gauss-Lobatto formula.
114  *
115  * References:
116  * This algorithm is a C++ implementation of the algorithm outlined in
117  *
118  * W. Gander and W. Gautschi, Adaptive Quadrature - Revisited.
119  * BIT, 40(1):84-101, March 2000. CS technical report:
120  * ftp.inf.ethz.ch/pub/publications/tech-reports/3xx/306.ps.gz
121  *
122  * The original MATLAB version can be downloaded here
124  *
125  * This particular implementation is based on code in QuantLib,
126  * a free-software/open-source library for financial quantitative
127  * analysts and developers - http://quantlib.org/
128  *
129  * \ingroup libcore
130  * \ingroup libpython
131  */
133 public:
134  typedef boost::function<Float (Float)> Integrand;
135
136  /**
137  * Initialize a Gauss-Lobatto integration scheme
138  *
139  * \param maxEvals Maximum number of function evaluations. The
140  * integrator will print a warning when this limit is
141  * exceeded. It will then stop the recursion, but a few
142  * further evaluations may still take place. Hence the limit
143  * is not a strict one.
144  *
145  * \param absError Absolute error requirement (0 to disable)
146  * \param relError Relative error requirement (0 to disable)
147  *
148  * \param useConvergenceEstimate Estimate the convergence behavior
149  * of the GL-quadrature by comparing the 4, 7 and 13-point
150  * variants and increase the absolute tolerance accordingly.
151  *
152  * \param warn Should the integrator warn when the number of
153  * function evaluations is exceeded?
154  */
155  GaussLobattoIntegrator(size_t maxEvals,
156  Float absError = 0,
157  Float relError = 0,
158  bool useConvergenceEstimate = true,
159  bool warn = true);
160
161  /**
162  * \brief Integrate the function \c f from \c a to \c b.
163  *
164  * Also returns the total number of evaluations if requested
165  */
166  Float integrate(const Integrand &f, Float a, Float b,
167  size_t *evals = NULL) const;
168 protected:
169  /**
170  * \brief Perform one step of the 4-point Gauss-Lobatto rule, then
171  * compute the same integral using a 7-point Kronrod extension and
172  * compare. If the accuracy is deemed too low, recurse.
173  *
174  * \param f Function to integrate
175  * \param a Lower integration limit
176  * \param b Upper integration limit
177  * \param fa Function evaluated at the lower limit
178  * \param fb Function evaluated at the upper limit
179  * \param is Absolute tolerance in epsilons
180  */
181  Float adaptiveGaussLobattoStep(const boost::function<Float (Float)>& f,
182  Float a, Float b, Float fa, Float fb, Float is, size_t &evals) const;
183
184  /**
185  * Compute the absolute error tolerance using a 13-point
186  * Gauss-Lobatto rule.
187  */
188  Float calculateAbsTolerance(const boost::function<Float (Float)>& f,
189  Float a, Float b, size_t &evals) const;
190 protected:
191  Float m_absError, m_relError;
192  size_t m_maxEvals;
194  bool m_warn;
195  static const Float m_alpha;
196  static const Float m_beta;
197  static const Float m_x1;
198  static const Float m_x2;
199  static const Float m_x3;
200 };
201
202 /**
203  * \brief Adaptively computes the integral of a multidimensional function using
204  * either a Gauss-Kronod (1D) or a Genz-Malik (>1D) cubature rule.
205  *
206  * This class is a C++ wrapper around the \c cubature code by Steven G. Johnson
207  * (http://ab-initio.mit.edu/wiki/index.php/Cubature)
208  *
209  * The original implementation is based on algorithms proposed in
210  *
211  * A. C. Genz and A. A. Malik, "An adaptive algorithm for numeric integration
212  * over an N-dimensional rectangular region," J. Comput. Appl. Math. 6 (4),
213  * 295–302 (1980).
214  *
215  * and
216  *
217  * J. Berntsen, T. O. Espelid, and A. Genz, "An adaptive algorithm for the
218  * approximate calculation of multiple integrals," ACM Trans. Math. Soft. 17
219  * (4), 437–451 (1991).
220  *
221  * \ingroup libcore
222  */
224 public:
225  typedef boost::function<void (const Float *, Float *)> Integrand;
226  typedef boost::function<void (size_t, const Float *, Float *)> VectorizedIntegrand;
227
228  enum EResult {
229  ESuccess = 0,
230  EFailure = 1
231  };
232
233  /**
234  * Initialize the Cubature integration scheme
235  *
236  * \param fDim Number of integrands (i.e. dimensions of the image space)
237  * \param nDim Number of integration dimensions (i.e. dimensions of the
238  * function domain)
239  * \param maxEvals Maximum number of function evaluations (0 means no
240  * limit). The error bounds will likely be exceeded when the
241  * integration is forced to stop prematurely. Note: the actual
242  * number of evaluations may somewhat exceed this value.
243  * \param absError Absolute error requirement (0 to disable)
244  * \param relError Relative error requirement (0 to disable)
245  */
246  NDIntegrator(size_t fDim, size_t dim,
247  size_t maxEvals, Float absError = 0, Float relError = 0);
248
249  /**
250  * \brief Integrate the function \c f over the rectangular domain
251  * bounded by \c min and \c max.
252  *
253  * The supplied function should have the interface
254  *
255  * <code>
256  * void integrand(const Float *in, Float *out);
257  * </code>
258  *
259  * The input array \c in consists of one set of input parameters
260  * having \c dim entries. The function is expected to store the
261  * results of the evaluation into the \c out array using \c fDim entries.
262  */
263  EResult integrate(const Integrand &f, const Float *min, const Float *max,
264  Float *result, Float *error, size_t *evals = NULL) const;
265
266  /**
267  * \brief Integrate the function \c f over the rectangular domain
268  * bounded by \c min and \c max.
269  *
270  * This function implements a vectorized version of the above
271  * integration function, which is more efficient by evaluating
272  * the integrant in `batches'. The supplied function should
273  * have the interface
274  *
275  * <code>
276  * void integrand(int numPoints, const Float *in, Float *out);
277  * </code>
278  *
279  * Note that \c in in is not a single point, but an array of \c numPoints points
280  * (length \c numPoints x \c dim), and upon return the values of all \c fDim
281  * integrands at all \c numPoints points should be stored in \c out
282  * (length \c fDim x \c numPoints). In particular, out[i*dim + j] is the j-th
283  * coordinate of the i-th point, and the k-th function evaluation (k<fDim)
284  * for the i-th point is returned in out[k*npt + i].
285  * The size of \c numPoints will vary with the dimensionality of the problem;
286  * higher-dimensional problems will have (exponentially) larger numbers,
287  * allowing for the possibility of more parallelism. Currently, \c numPoints
288  * starts at 15 in 1d, 17 in 2d, and 33 in 3d, but as the integrator
289  * calls your integrand more and more times the value will grow. e.g. if you end
290  * up requiring several thousand points in total, \c numPoints may grow to
291  * several hundred.
292  */
293  EResult integrateVectorized(const VectorizedIntegrand &f, const Float *min,
294  const Float *max, Float *result, Float *error, size_t *evals = NULL) const;
295 protected:
296  size_t m_fdim, m_dim, m_maxEvals;
297  Float m_absError, m_relError;
298 };
299
300 //! @}
301 // -----------------------------------------------------------------------
302
304
305 #endif /* __MITSUBA_CORE_QUAD_H_ */
bool m_useConvergenceEstimate
void gaussLobatto(int n, Float *nodes, Float *weights)
Computes the nodes and weights of a Gauss-Lobatto quadrature rule with the given number of evaluation...
size_t m_maxEvals
std::pair< double, double > legendrePD(int l, double x)
Evaluate the l-th Legendre polynomial and its derivative using recurrence (double precision) ...
double legendreP(int l, int m, double x)
Evaluate an associated Legendre polynomial using recurrence (double precision)
#define MTS_EXPORT_CORE
Definition: getopt.h:29
static const Float m_beta
Float m_relError
static const Float m_x1
#define MTS_NAMESPACE_BEGIN
Definition: platform.h:137
size_t m_maxEvals
boost::function< void(const Float *, Float *)> Integrand
Computes the integral of a one-dimensional function using adaptive Gauss-Lobatto quadrature.
static const Float m_x3
Adaptively computes the integral of a multidimensional function using either a Gauss-Kronod (1D) or a...
EResult
void gaussLegendre(int n, Float *nodes, Float *weights)
Computes the nodes and weights of a Gauss-Legendre quadrature (aka &quot;Gaussian quadrature&quot;) rule with t...
boost::function< Float(Float)> Integrand