326 lines
9.5 KiB
C++
326 lines
9.5 KiB
C++
|
// Copyright 2005, Ondra Kamenik
|
||
|
|
||
|
// Quadrature.
|
||
|
|
||
|
/* This file defines an interface for one dimensional (non-nested) quadrature
|
||
|
|OneDQuadrature|, and a parent for all multi-dimensional
|
||
|
quadratures. This parent class |Quadrature| presents a general concept of
|
||
|
quadrature, this is
|
||
|
$$\int f(x){\rm d}x \approx\sum_{i=1}^N w_ix_i$$
|
||
|
The class |Quadrature| just declares this concept. The concept is
|
||
|
implemented by class |QuadratureImpl| which paralelizes the
|
||
|
summation. All implementations therefore wishing to use the parallel
|
||
|
implementation should inherit from |QuadratureImpl| and integration is
|
||
|
done.
|
||
|
|
||
|
The integration concept relies on a point iterator, which goes through
|
||
|
all $x_i$ and $w_i$ for $i=1,\ldots,N$. All the iterators must be able
|
||
|
to go through only a portion of the set $i=1,\ldots,N$. This enables
|
||
|
us to implement paralelism, for two threads for example, one iterator
|
||
|
goes from the beginning to the (approximately) half, and the other
|
||
|
goes from the half to the end.
|
||
|
|
||
|
Besides this concept of the general quadrature, this file defines also
|
||
|
one dimensional quadrature, which is basically a scheme of points and
|
||
|
weights for different levels. The class |OneDQuadrature| is a parent
|
||
|
of all such objects, the classes |GaussHermite| and |GaussLegendre|
|
||
|
are specific implementations for Gauss--Hermite and Gauss--Legendre
|
||
|
quadratures resp. */
|
||
|
|
||
|
#ifndef QUADRATURE_H
|
||
|
#define QUADRATURE_H
|
||
|
|
||
|
#include <cstdlib>
|
||
|
#include "vector_function.hh"
|
||
|
#include "int_sequence.h"
|
||
|
#include "sthread.h"
|
||
|
|
||
|
/* This pure virtual class represents a concept of one-dimensional
|
||
|
(non-nested) quadrature. So, one dimensional quadrature must return
|
||
|
number of levels, number of points in a given level, and then a point
|
||
|
and a weight in a given level and given order. */
|
||
|
|
||
|
class OneDQuadrature
|
||
|
{
|
||
|
public:
|
||
|
virtual ~OneDQuadrature()
|
||
|
{
|
||
|
}
|
||
|
virtual int numLevels() const = 0;
|
||
|
virtual int numPoints(int level) const = 0;
|
||
|
virtual double point(int level, int i) const = 0;
|
||
|
virtual double weight(int lelel, int i) const = 0;
|
||
|
};
|
||
|
|
||
|
/* This is a general concept of multidimensional quadrature. at this
|
||
|
general level, we maintain only a dimension, and declare virtual
|
||
|
functions for integration. The function take two forms; first takes a
|
||
|
constant |VectorFunction| as an argument, creates locally
|
||
|
|VectorFunctionSet| and do calculation, second one takes as an
|
||
|
argument |VectorFunctionSet|.
|
||
|
|
||
|
Part of the interface is a method returning a number of evaluations
|
||
|
for a specific level. Note two things: this assumes that the number of
|
||
|
evaluations is known apriori and thus it is not applicable for
|
||
|
adaptive quadratures, second for Monte Carlo type of quadrature, the
|
||
|
level is a number of evaluations. */
|
||
|
|
||
|
class Quadrature
|
||
|
{
|
||
|
protected:
|
||
|
int dim;
|
||
|
public:
|
||
|
Quadrature(int d) : dim(d)
|
||
|
{
|
||
|
}
|
||
|
virtual ~Quadrature()
|
||
|
{
|
||
|
}
|
||
|
int
|
||
|
dimen() const
|
||
|
{
|
||
|
return dim;
|
||
|
}
|
||
|
virtual void integrate(const VectorFunction &func, int level,
|
||
|
int tn, Vector &out) const = 0;
|
||
|
virtual void integrate(VectorFunctionSet &fs, int level, Vector &out) const = 0;
|
||
|
virtual int numEvals(int level) const = 0;
|
||
|
};
|
||
|
|
||
|
/* This is just an integration worker, which works over a given
|
||
|
|QuadratureImpl|. It also needs the function, level, a specification
|
||
|
of the subgroup of points, and output vector.
|
||
|
|
||
|
See |@<|QuadratureImpl| class declaration@>| for details. */
|
||
|
|
||
|
template <typename _Tpit>
|
||
|
class QuadratureImpl;
|
||
|
|
||
|
template <typename _Tpit>
|
||
|
class IntegrationWorker : public THREAD
|
||
|
{
|
||
|
const QuadratureImpl<_Tpit> &quad;
|
||
|
VectorFunction &func;
|
||
|
int level;
|
||
|
int ti;
|
||
|
int tn;
|
||
|
Vector &outvec;
|
||
|
public:
|
||
|
IntegrationWorker(const QuadratureImpl<_Tpit> &q, VectorFunction &f, int l,
|
||
|
int tii, int tnn, Vector &out)
|
||
|
: quad(q), func(f), level(l), ti(tii), tn(tnn), outvec(out)
|
||
|
{
|
||
|
}
|
||
|
|
||
|
/* This integrates the given portion of the integral. We obtain first
|
||
|
and last iterators for the portion (|beg| and |end|). Then we iterate
|
||
|
through the portion. and finally we add the intermediate result to the
|
||
|
result |outvec|.
|
||
|
|
||
|
This method just everything up as it is coming. This might be imply
|
||
|
large numerical errors, perhaps in future I will implement something
|
||
|
smarter. */
|
||
|
|
||
|
void
|
||
|
operator()()
|
||
|
{
|
||
|
_Tpit beg = quad.begin(ti, tn, level);
|
||
|
_Tpit end = quad.begin(ti+1, tn, level);
|
||
|
Vector tmpall(outvec.length());
|
||
|
tmpall.zeros();
|
||
|
Vector tmp(outvec.length());
|
||
|
|
||
|
// note that since beg came from begin, it has empty signal
|
||
|
// and first evaluation gets no signal
|
||
|
for (_Tpit run = beg; run != end; ++run)
|
||
|
{
|
||
|
func.eval(run.point(), run.signal(), tmp);
|
||
|
tmpall.add(run.weight(), tmp);
|
||
|
}
|
||
|
|
||
|
{
|
||
|
SYNCHRO syn(&outvec, "IntegrationWorker");
|
||
|
outvec.add(1.0, tmpall);
|
||
|
}
|
||
|
}
|
||
|
};
|
||
|
|
||
|
/* This is the class which implements the integration. The class is
|
||
|
templated by the iterator type. We declare a method |begin| returning
|
||
|
an iterator to the beginnning of the |ti|-th portion out of total |tn|
|
||
|
portions for a given level.
|
||
|
|
||
|
In addition, we define a method which saves all the points to a given
|
||
|
file. Only for debugging purposes. */
|
||
|
|
||
|
template <typename _Tpit>
|
||
|
class QuadratureImpl : public Quadrature
|
||
|
{
|
||
|
friend class IntegrationWorker<_Tpit>;
|
||
|
public:
|
||
|
QuadratureImpl(int d) : Quadrature(d)
|
||
|
{
|
||
|
}
|
||
|
|
||
|
/* Just fill a thread group with workes and run it. */
|
||
|
void
|
||
|
integrate(VectorFunctionSet &fs, int level, Vector &out) const
|
||
|
{
|
||
|
// todo: out.length()==func.outdim()
|
||
|
// todo: dim == func.indim()
|
||
|
out.zeros();
|
||
|
THREAD_GROUP gr;
|
||
|
for (int ti = 0; ti < fs.getNum(); ti++)
|
||
|
{
|
||
|
gr.insert(new IntegrationWorker<_Tpit>(*this, fs.getFunc(ti),
|
||
|
level, ti, fs.getNum(), out));
|
||
|
}
|
||
|
gr.run();
|
||
|
}
|
||
|
void
|
||
|
integrate(const VectorFunction &func,
|
||
|
int level, int tn, Vector &out) const
|
||
|
{
|
||
|
VectorFunctionSet fs(func, tn);
|
||
|
integrate(fs, level, out);
|
||
|
}
|
||
|
|
||
|
/* Just for debugging. */
|
||
|
void
|
||
|
savePoints(const char *fname, int level) const
|
||
|
{
|
||
|
FILE *fd;
|
||
|
if (NULL == (fd = fopen(fname, "w")))
|
||
|
{
|
||
|
// todo: raise
|
||
|
fprintf(stderr, "Cannot open file %s for writing.\n", fname);
|
||
|
exit(1);
|
||
|
}
|
||
|
_Tpit beg = begin(0, 1, level);
|
||
|
_Tpit end = begin(1, 1, level);
|
||
|
for (_Tpit run = beg; run != end; ++run)
|
||
|
{
|
||
|
fprintf(fd, "%16.12g", run.weight());
|
||
|
for (int i = 0; i < dimen(); i++)
|
||
|
fprintf(fd, "\t%16.12g", run.point()[i]);
|
||
|
fprintf(fd, "\n");
|
||
|
}
|
||
|
fclose(fd);
|
||
|
}
|
||
|
|
||
|
_Tpit
|
||
|
start(int level) const
|
||
|
{
|
||
|
return begin(0, 1, level);
|
||
|
}
|
||
|
_Tpit
|
||
|
end(int level) const
|
||
|
{
|
||
|
return begin(1, 1, level);
|
||
|
}
|
||
|
protected:
|
||
|
virtual _Tpit begin(int ti, int tn, int level) const = 0;
|
||
|
};
|
||
|
|
||
|
/* This is only an interface to a precalculated data in file {\tt
|
||
|
precalc\_quadrature.dat} which is basically C coded static data. It
|
||
|
implements |OneDQuadrature|. The data file is supposed to define the
|
||
|
following data: number of levels, array of number of points at each
|
||
|
level, an array of weights and array of points. The both latter array
|
||
|
store data level by level. An offset for a specific level is stored in
|
||
|
|offsets| integer sequence.
|
||
|
|
||
|
The implementing subclasses just fill the necessary data from the
|
||
|
file, the rest is calculated here. */
|
||
|
|
||
|
class OneDPrecalcQuadrature : public OneDQuadrature
|
||
|
{
|
||
|
int num_levels;
|
||
|
const int *num_points;
|
||
|
const double *weights;
|
||
|
const double *points;
|
||
|
IntSequence offsets;
|
||
|
public:
|
||
|
OneDPrecalcQuadrature(int nlevels, const int *npoints,
|
||
|
const double *wts, const double *pts)
|
||
|
: num_levels(nlevels), num_points(npoints),
|
||
|
weights(wts), points(pts), offsets(num_levels)
|
||
|
{
|
||
|
calcOffsets();
|
||
|
}
|
||
|
virtual ~OneDPrecalcQuadrature()
|
||
|
{
|
||
|
}
|
||
|
int
|
||
|
numLevels() const
|
||
|
{
|
||
|
return num_levels;
|
||
|
}
|
||
|
int
|
||
|
numPoints(int level) const
|
||
|
{
|
||
|
return num_points[level-1];
|
||
|
}
|
||
|
double
|
||
|
point(int level, int i) const
|
||
|
{
|
||
|
return points[offsets[level-1]+i];
|
||
|
}
|
||
|
double
|
||
|
weight(int level, int i) const
|
||
|
{
|
||
|
return weights[offsets[level-1]+i];
|
||
|
}
|
||
|
protected:
|
||
|
void calcOffsets();
|
||
|
};
|
||
|
|
||
|
/* Just precalculated Gauss--Hermite quadrature. This quadrature integrates exactly integrals
|
||
|
$$\int_{-\infty}^{\infty} x^ke^{-x^2}{\rm d}x$$
|
||
|
for level $k$.
|
||
|
|
||
|
Note that if pluging this one-dimensional quadrature to product or
|
||
|
Smolyak rule in order to integrate a function $f$ through normally
|
||
|
distributed inputs, one has to wrap $f$ to
|
||
|
|GaussConverterFunction| and apply the product or Smolyak rule to the
|
||
|
new function.
|
||
|
|
||
|
Check {\tt precalc\_quadrature.dat} for available levels. */
|
||
|
|
||
|
class GaussHermite : public OneDPrecalcQuadrature
|
||
|
{
|
||
|
public:
|
||
|
GaussHermite();
|
||
|
};
|
||
|
|
||
|
/* Just precalculated Gauss--Legendre quadrature. This quadrature integrates exactly integrals
|
||
|
$$\int_0^1x^k{\rm d}x$$
|
||
|
for level $k$.
|
||
|
|
||
|
Check {\tt precalc\_quadrature.dat} for available levels. */
|
||
|
|
||
|
class GaussLegendre : public OneDPrecalcQuadrature
|
||
|
{
|
||
|
public:
|
||
|
GaussLegendre();
|
||
|
};
|
||
|
|
||
|
/* This is just an inverse cummulative density function of normal
|
||
|
distribution. Its only method |get| returns for a given number
|
||
|
$x\in(0,1)$ a number $y$ such that $P(z<y)=x$, where the probability
|
||
|
is taken over normal distribution $N(0,1)$.
|
||
|
|
||
|
Currently, the implementation is done by a table lookup which implies
|
||
|
that the tails had to be chopped off. This further implies that Monte
|
||
|
Carlo quadratures using this transformation to draw points from
|
||
|
multidimensional $N(0,I)$ fail to integrate with satisfactory
|
||
|
precision polynomial functions, for which the tails matter. */
|
||
|
|
||
|
class NormalICDF
|
||
|
{
|
||
|
public:
|
||
|
static double get(double x);
|
||
|
};
|
||
|
|
||
|
#endif
|