quadrature.hh 9.62 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34
// 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"
35 36
#include "int_sequence.hh"
#include "sthread.hh"
37 38 39 40 41 42 43 44 45 46

/* 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()
47
  = default;
48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75
  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()
76
  = default;
77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122
  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
123
  operator()() override
124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164
  {
    _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
165
  integrate(VectorFunctionSet &fs, int level, Vector &out) const override
166 167 168 169 170 171 172 173 174 175 176 177 178 179
  {
    // 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,
180
            int level, int tn, Vector &out) const override
181 182 183 184 185 186 187 188 189 190
  {
    VectorFunctionSet fs(func, tn);
    integrate(fs, level, out);
  }

  /* Just for debugging. */
  void
  savePoints(const char *fname, int level) const
  {
    FILE *fd;
191
    if (nullptr == (fd = fopen(fname, "w")))
192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248
      {
        // 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();
  }
249 250
  ~OneDPrecalcQuadrature()
  override = default;
251
  int
252
  numLevels() const override
253 254 255 256
  {
    return num_levels;
  }
  int
257
  numPoints(int level) const override
258 259 260 261
  {
    return num_points[level-1];
  }
  double
262
  point(int level, int i) const override
263 264 265 266
  {
    return points[offsets[level-1]+i];
  }
  double
267
  weight(int level, int i) const override
268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322
  {
    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