dynare++/integ/cc/product.hh

-// Copyright 2005, Ondra Kamenik
-
-// Product quadrature.
-
-/* This file defines a product multidimensional quadrature. If $Q_k$
-   denotes the one dimensional quadrature, then the product quadrature
-   $Q$ of $k$ level and dimension $d$ takes the form
-   $$Qf=\sum_{i_1=1}^{n_k}\ldots\sum_{i_d=1}^{n^k}w_{i_1}\cdot\ldots\cdot w_{i_d}
-   f(x_{i_1},\ldots,x_{i_d})$$
-   which can be written in terms of the one dimensional quadrature $Q_k$ as
-   $$Qf=(Q_k\otimes\ldots\otimes Q_k)f$$
-
-   Here we define the product quadrature iterator |prodpit| and plug it
-   into |QuadratureImpl| to obtains |ProductQuadrature|. */
-
-#ifndef PRODUCT_H
-#define PRODUCT_H
-
-#include "int_sequence.hh"
-#include "vector_function.hh"
-#include "quadrature.hh"
-
-/* This defines a product point iterator. We have to maintain the
-   following: a pointer to product quadrature in order to know the
-   dimension and the underlying one dimensional quadrature, then level,
-   number of points in the level, integer sequence of indices, signal,
-   the coordinates of the point and the weight.
-
-   The point indices, signal, and point coordinates are implmented as
-   pointers in order to allow for empty constructor.
-
-   The constructor |prodpit(const ProductQuadrature& q, int j0, int l)|
-   constructs an iterator pointing to $(j0,0,\ldots,0)$, which is used by
-   |begin| dictated by |QuadratureImpl|. */
-
-class ProductQuadrature;
-
-class prodpit
-{
-protected:
-  const ProductQuadrature &prodq;
-  int level{0};
-  int npoints{0};
-  IntSequence jseq;
-  bool end_flag{true};
-  ParameterSignal sig;
-  Vector p;
-  double w;
-public:
-  prodpit() = default;
-  prodpit(const ProductQuadrature &q, int j0, int l);
-  prodpit(const prodpit &ppit) = default;
-  ~prodpit() = default;
-  bool operator==(const prodpit &ppit) const;
-  bool
-  operator!=(const prodpit &ppit) const
-  {
-    return !operator==(ppit);
-  }
-  prodpit &operator=(const prodpit &spit) = delete;
-  prodpit &operator++();
-  const ParameterSignal &
-  signal() const
-  {
-    return sig;
-  }
-  const Vector &
-  point() const
-  {
-    return p;
-  }
-  double
-  weight() const
-  {
-    return w;
-  }
-  void print() const;
-protected:
-  void setPointAndWeight();
-};
-
-/* The product quadrature is just |QuadratureImpl| with the product
-   iterator plugged in. The object is constructed by just giving the
-   underlying one dimensional quadrature, and the dimension. The only
-   extra method is |designLevelForEvals| which for the given maximum
-   number of evaluations (and dimension and underlying quadrature from
-   the object) returns a maximum level yeilding number of evaluations
-   less than the given number. */
-
-class ProductQuadrature : public QuadratureImpl<prodpit>
-{
-  friend class prodpit;
-  const OneDQuadrature &uquad;
-public:
-  ProductQuadrature(int d, const OneDQuadrature &uq);
-  ~ProductQuadrature() override = default;
-  int
-  numEvals(int l) const override
-  {
-    int res = 1;
-    for (int i = 0; i < dimen(); i++)
-      res *= uquad.numPoints(l);
-    return res;
-  }
-  void designLevelForEvals(int max_eval, int &lev, int &evals) const;
-protected:
-  prodpit begin(int ti, int tn, int level) const override;
-};
-
-#endif

dynare++/integ/cc/quadrature.cc

-// Copyright 2005, Ondra Kamenik
-
-#include "quadrature.hh"
-#include "precalc_quadrature.hh"
-
-#include <cmath>
-
-void
-OneDPrecalcQuadrature::calcOffsets()
-{
-  offsets[0] = 0;
-  for (int i = 1; i < num_levels; i++)
-    offsets[i] = offsets[i-1] + num_points[i-1];
-}
-
-GaussHermite::GaussHermite()
-  : OneDPrecalcQuadrature(gh_num_levels, gh_num_points, gh_weights, gh_points)
-{
-}
-
-GaussLegendre::GaussLegendre()
-  : OneDPrecalcQuadrature(gl_num_levels, gl_num_points, gl_weights, gl_points)
-{
-}

dynare++/integ/cc/quadrature.hh

-// 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 <iostream>
-#include <fstream>
-#include <iomanip>
-
-#include "vector_function.hh"
-#include "int_sequence.hh"
-#include "sthread.hh"
-
-/* 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() = default;
-  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() = default;
-  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 sthread::detach_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()(std::mutex &mut) override
-  {
-    _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);
-      }
-
-    {
-      std::unique_lock<std::mutex> lk{mut};
-      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 override
-  {
-    // todo: out.length()==func.outdim()
-    // todo: dim == func.indim()
-    out.zeros();
-    sthread::detach_thread_group gr;
-    for (int ti = 0; ti < fs.getNum(); ti++)
-      gr.insert(std::make_unique<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 override
-  {
-    VectorFunctionSet fs(func, tn);
-    integrate(fs, level, out);
-  }
-
-  /* Just for debugging. */
-  void
-  savePoints(const std::string &fname, int level) const
-  {
-    std::ofstream fd{fname, std::ios::out | std::ios::trunc};
-    if (fd.fail())
-      {
-        // todo: raise
-        std::cerr << "Cannot open file " << fname << " for writing." << std::endl;
-        exit(EXIT_FAILURE);
-      }
-    _Tpit beg = begin(0, 1, level);
-    _Tpit end = begin(1, 1, level);
-    fd << std::setprecision(12);
-    for (_Tpit run = beg; run != end; ++run)
-      {
-        fd << std::setw(16) << run.weight();
-        for (int i = 0; i < dimen(); i++)
-          fd << '\t' << std::setw(16) << run.point()[i];
-        fd << std::endl;
-      }
-    fd.close();
-  }
-
-  _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();
-  }
-  ~OneDPrecalcQuadrature() override = default;
-  int
-  numLevels() const override
-  {
-    return num_levels;
-  }
-  int
-  numPoints(int level) const override
-  {
-    return num_points[level-1];
-  }
-  double
-  point(int level, int i) const override
-  {
-    return points[offsets[level-1]+i];
-  }
-  double
-  weight(int level, int i) const override
-  {
-    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();
-};
-
-#endif

dynare++/integ/cc/quasi_mcarlo.cc

-// Copyright 2005, Ondra Kamenik
-
-#include "quasi_mcarlo.hh"
-
-#include <cmath>
-#include <iostream>
-#include <iomanip>
-
-/* Here in the constructor, we have to calculate a maximum length of
-   |coeff| array for a given |base| and given maximum |maxn|. After
-   allocation, we calculate the coefficients. */
-
-RadicalInverse::RadicalInverse(int n, int b, int mxn)
-  : num(n), base(b), maxn(mxn),
-    coeff(static_cast<int>(floor(log(static_cast<double>(maxn))/log(static_cast<double>(b)))+2), 0)
-{
-  int nr = num;
-  j = -1;
-  do
-    {
-      j++;
-      coeff[j] = nr % base;
-      nr = nr / base;
-    }
-  while (nr > 0);
-}
-
-/* This evaluates the radical inverse. If there was no permutation, we have to calculate
-   $$
-   {c_0\over b}+{c_1\over b^2}+\ldots+{c_j\over b^{j+1}}
-   $$
-   which is evaluated as
-   $$
-   \left(\ldots\left(\left({c_j\over b}\cdot{1\over b}+{c_{j-1}\over b}\right)
-   \cdot{1\over b}+{c_{j-2}\over b}\right)
-   \ldots\right)\cdot{1\over b}+{c_0\over b}
-   $$
-   Now with permutation $\pi$, we have
-   $$
-   \left(\ldots\left(\left({\pi(c_j)\over b}\cdot{1\over b}+
-   {\pi(c_{j-1})\over b}\right)\cdot{1\over b}+
-   {\pi(c_{j-2})\over b}\right)
-   \ldots\right)\cdot{1\over b}+{\pi(c_0)\over b}
-   $$ */
-
-double
-RadicalInverse::eval(const PermutationScheme &p) const
-{
-  double res = 0;
-  for (int i = j; i >= 0; i--)
-    {
-      int cper = p.permute(i, base, coeff[i]);
-      res = (cper + res)/base;
-    }
-  return res;
-}
-
-/* We just add 1 to the lowest coefficient and check for overflow with respect
-   to the base. */
-void
-RadicalInverse::increase()
-{
-  // todo: raise if |num+1 > maxn|
-  num++;
-  int i = 0;
-  coeff[i]++;
-  while (coeff[i] == base)
-    {
-      coeff[i] = 0;
-      coeff[++i]++;
-    }
-  if (i > j)
-    j = i;
-}
-
-/* Debug print. */
-void
-RadicalInverse::print() const
-{
-  std::cout << "n=" << num << " b=" << base << " c=";
-  coeff.print();
-}
-
-/* Here we have the first 170 primes. This means that we are not able
-   to integrate dimensions greater than 170. */
-
-std::array<int, 170> HaltonSequence::primes = {
-  2,     3,     5,     7,    11,    13,    17,    19,    23,    29,
-  31,    37,    41,    43,    47,    53,    59,    61,    67,    71,
-  73,    79,    83,    89,    97,   101,   103,   107,   109,   113,
-  127,   131,   137,   139,   149,   151,   157,   163,   167,   173,
-  179,   181,   191,   193,   197,   199,   211,   223,   227,   229,
-  233,   239,   241,   251,   257,   263,   269,   271,   277,   281,
-  283,   293,   307,   311,   313,   317,   331,   337,   347,   349,
-  353,   359,   367,   373,   379,   383,   389,   397,   401,   409,
-  419,   421,   431,   433,   439,   443,   449,   457,   461,   463,
-  467,   479,   487,   491,   499,   503,   509,   521,   523,   541,
-  547,   557,   563,   569,   571,   577,   587,   593,   599,   601,
-  607,   613,   617,   619,   631,   641,   643,   647,   653,   659,
-  661,   673,   677,   683,   691,   701,   709,   719,   727,   733,
-  739,   743,   751,   757,   761,   769,   773,   787,   797,   809,
-  811,   821,   823,   827,   829,   839,   853,   857,   859,   863,
-  877,   881,   883,   887,   907,   911,   919,   929,   937,   941,
-  947,   953,   967,   971,   977,   983,   991,   997,  1009,  1013
-};
-
-/* This takes first |dim| primes and constructs |dim| radical inverses
-   and calls |eval|. */
-
-HaltonSequence::HaltonSequence(int n, int mxn, int dim, const PermutationScheme &p)
-  : num(n), maxn(mxn), per(p), pt(dim)
-{
-  // todo: raise if |dim > num_primes|
-  // todo: raise if |n > mxn|
-  for (int i = 0; i < dim; i++)
-    ri.emplace_back(num, primes[i], maxn);
-  eval();
-}
-
-/* This calls |RadicalInverse::increase| for all radical inverses and
-   calls |eval|. */
-
-void
-HaltonSequence::increase()
-{
-  for (auto & i : ri)
-    i.increase();
-  num++;
-  if (num <= maxn)
-    eval();
-}
-
-/* This sets point |pt| to radical inverse evaluations in each dimension. */
-void
-HaltonSequence::eval()
-{
-  for (unsigned int i = 0; i < ri.size(); i++)
-    pt[i] = ri[i].eval(per);
-}
-
-/* Debug print. */
-void
-HaltonSequence::print() const
-{
-  auto ff = std::cout.flags();
-  for (const auto & i : ri)
-    i.print();
-  std::cout << "point=[ "
-            << std::fixed << std::setprecision(6);
-  for (unsigned int i = 0; i < ri.size(); i++)
-    std::cout << std::setw(7) << pt[i] << ' ';
-  std::cout << ']' << std::endl;
-  std::cout.flags(ff);
-}
-
-qmcpit::qmcpit(const QMCSpecification &s, int n)
-  : spec(s), halton{n, s.level(), s.dimen(), s.getPerScheme()},
-    sig{s.dimen()}
-{
-}
-
-bool
-qmcpit::operator==(const qmcpit &qpit) const
-{
-  return &spec == &qpit.spec && halton.getNum() == qpit.halton.getNum();
-}
-
-qmcpit &
-qmcpit::operator++()
-{
-  // todo: raise if |halton == null || qmcq == NULL|
-  halton.increase();
-  return *this;
-}
-
-double
-qmcpit::weight() const
-{
-  return 1.0/spec.level();
-}
-
-/* Clear from code. */
-int
-WarnockPerScheme::permute(int i, int base, int c) const
-{
-  return (c+i) % base;
-}
-
-/* Clear from code. */
-int
-ReversePerScheme::permute(int i, int base, int c) const
-{
-  return (base-c) % base;
-}

dynare++/integ/cc/quasi_mcarlo.hh

-// Copyright 2005, Ondra Kamenik
-
-// Quasi Monte Carlo quadrature.
-
-/* This defines quasi Monte Carlo quadratures for cube and for a function
-   multiplied by normal density. The quadrature for a cube is named
-   |QMCarloCubeQuadrature| and integrates:
-   $$\int_{\langle 0,1\rangle^n}f(x){\rm d}x$$
-   The quadrature for a function of normally distributed parameters is
-   named |QMCarloNormalQuadrature| and integrates:
-   $${1\over\sqrt{(2\pi)^n}}\int_{(-\infty,\infty)^n}f(x)e^{-{1\over 2}x^Tx}{\rm d}x$$
-
-   For a cube we define |qmcpit| as iterator of |QMCarloCubeQuadrature|,
-   and for the normal density multiplied function we define |qmcnpit| as
-   iterator of |QMCarloNormalQuadrature|.
-
-   The quasi Monte Carlo method generates low discrepancy points with
-   equal weights. The one dimensional low discrepancy sequences are
-   generated by |RadicalInverse| class, the sequences are combined for
-   higher dimensions by |HaltonSequence| class. The Halton sequence can
-   use a permutation scheme; |PermutattionScheme| is an abstract class
-   for all permutaton schemes. We have three implementations:
-   |WarnockPerScheme|, |ReversePerScheme|, and |IdentityPerScheme|. */
-
-#ifndef QUASI_MCARLO_H
-#define QUASI_MCARLO_H
-
-#include "int_sequence.hh"
-#include "quadrature.hh"
-
-#include "Vector.hh"
-
-#include <vector>
-
-/* This abstract class declares |permute| method which permutes
-   coefficient |c| having index of |i| fro the base |base| and returns
-   the permuted coefficient which must be in $0,\ldots,base-1$. */
-
-class PermutationScheme
-{
-public:
-  PermutationScheme() = default;
-  virtual ~PermutationScheme() = default;
-  virtual int permute(int i, int base, int c) const = 0;
-};
-
-/* This class represents an integer number |num| as
-   $c_0+c_1b+c_2b^2+\ldots+c_jb^j$, where $b$ is |base| and
-   $c_0,\ldots,c_j$ is stored in |coeff|. The size of |IntSequence| coeff
-   does not grow with growing |num|, but is fixed from the very beginning
-   and is set according to supplied maximum |maxn|.
-
-   The basic method is |eval| which evaluates the |RadicalInverse| with a
-   given permutation scheme and returns the point, and |increase| which
-   increases |num| and recalculates the coefficients. */
-
-class RadicalInverse
-{
-  int num;
-  int base;
-  int maxn;
-  int j;
-  IntSequence coeff;
-public:
-  RadicalInverse(int n, int b, int mxn);
-  RadicalInverse(const RadicalInverse &ri) = default;
-  RadicalInverse &operator=(const RadicalInverse &radi) = default;
-  double eval(const PermutationScheme &p) const;
-  void increase();
-  void print() const;
-};
-
-/* This is a vector of |RadicalInverse|s, each |RadicalInverse| has a
-   different prime as its base. The static members |primes| and
-   |num_primes| define a precalculated array of primes. The |increase|
-   method of the class increases indices in all |RadicalInverse|s and
-   sets point |pt| to contain the points in each dimension. */
-
-class HaltonSequence
-{
-private:
-  static std::array<int, 170> primes;
-protected:
-  int num;
-  int maxn;
-  std::vector<RadicalInverse> ri;
-  const PermutationScheme &per;
-  Vector pt;
-public:
-  HaltonSequence(int n, int mxn, int dim, const PermutationScheme &p);
-  HaltonSequence(const HaltonSequence &hs) = default;
-  HaltonSequence &operator=(const HaltonSequence &hs) = delete;
-  void increase();
-  const Vector &
-  point() const
-  {
-    return pt;
-  }
-  const int
-  getNum() const
-  {
-    return num;
-  }
-  void print() const;
-protected:
-  void eval();
-};
-
-/* This is a specification of quasi Monte Carlo quadrature. It consists
-   of dimension |dim|, number of points (or level) |lev|, and the
-   permutation scheme. This class is common to all quasi Monte Carlo
-   classes. */
-
-class QMCSpecification
-{
-protected:
-  int dim;
-  int lev;
-  const PermutationScheme &per_scheme;
-public:
-  QMCSpecification(int d, int l, const PermutationScheme &p)
-    : dim(d), lev(l), per_scheme(p)
-  {
-  }
-  virtual ~QMCSpecification() = default;
-  int
-  dimen() const
-  {
-    return dim;
-  }
-  int
-  level() const
-  {
-    return lev;
-  }
-  const PermutationScheme &
-  getPerScheme() const
-  {
-    return per_scheme;
-  }
-};
-
-/* This is an iterator for quasi Monte Carlo over a cube
-   |QMCarloCubeQuadrature|. The iterator maintains |HaltonSequence| of
-   the same dimension as given by the specification. An iterator can be
-   constructed from a given number |n|, or by a copy constructor. */
-
-class qmcpit
-{
-protected:
-  const QMCSpecification &spec;
-  HaltonSequence halton;
-  ParameterSignal sig;
-public:
-  qmcpit(const QMCSpecification &s, int n);
-  qmcpit(const qmcpit &qpit) = default;
-  virtual ~qmcpit() = default;
-  bool operator==(const qmcpit &qpit) const;
-  bool
-  operator!=(const qmcpit &qpit) const
-  {
-    return !operator==(qpit);
-  }
-  qmcpit &operator=(const qmcpit &qpit) = delete;
-  qmcpit &operator++();
-  const ParameterSignal &
-  signal() const
-  {
-    return sig;
-  }
-  const Vector &
-  point() const
-  {
-    return halton.point();
-  }
-  double weight() const;
-  void
-  print() const
-  {
-    halton.print();
-  }
-};
-
-/* This is an easy declaration of quasi Monte Carlo quadrature for a
-   cube. Everything important has been done in its iterator |qmcpit|, so
-   we only inherit from general |Quadrature| and reimplement |begin| and
-   |numEvals|. */
-
-class QMCarloCubeQuadrature : public QuadratureImpl<qmcpit>, public QMCSpecification
-{
-public:
-  QMCarloCubeQuadrature(int d, int l, const PermutationScheme &p)
-    : QuadratureImpl<qmcpit>(d), QMCSpecification(d, l, p)
-  {
-  }
-  ~QMCarloCubeQuadrature() override = default;
-  int
-  numEvals(int l) const override
-  {
-    return l;
-  }
-protected:
-  qmcpit
-  begin(int ti, int tn, int lev) const override
-  {
-    return qmcpit(*this, ti*level()/tn + 1);
-  }
-};
-
-/* Declares Warnock permutation scheme. */
-class WarnockPerScheme : public PermutationScheme
-{
-public:
-  int permute(int i, int base, int c) const override;
-};
-
-/* Declares reverse permutation scheme. */
-class ReversePerScheme : public PermutationScheme
-{
-public:
-  int permute(int i, int base, int c) const override;
-};
-
-/* Declares no permutation (identity) scheme. */
-class IdentityPerScheme : public PermutationScheme
-{
-public:
-  int
-  permute(int i, int base, int c) const override
-  {
-    return c;
-  }
-};
-
-#endif

dynare++/integ/cc/smolyak.cc

-// Copyright 2005, Ondra Kamenik
-
-#include "smolyak.hh"
-#include "symmetry.hh"
-
-#include <iostream>
-#include <iomanip>
-
-/* This constructs a beginning of |isum| summand in |smolq|. We must be
-   careful here, since |isum| can be past-the-end, so no reference to
-   vectors in |smolq| by |isum| must be done in this case. */
-
-smolpit::smolpit(const SmolyakQuadrature &q, unsigned int isum)
-  : smolq(q), isummand(isum),
-    jseq{q.dimen(), 0},
-    sig{q.dimen()},
-    p{q.dimen()}
-{
-  if (isummand < q.numSummands())
-    setPointAndWeight();
-}
-
-bool
-smolpit::operator==(const smolpit &spit) const
-{
-  return &smolq == &spit.smolq && isummand == spit.isummand && jseq == spit.jseq;
-}
-
-/* We first try to increase index within the current summand. If we are
-   at maximum, we go to a subsequent summand. Note that in this case all
-   indices in |jseq| will be zero, so no change is needed. */
-
-smolpit &
-smolpit::operator++()
-{
-  // todo: throw if |smolq==NULL| or |jseq==NULL| or |sig==NULL|
-  const IntSequence &levpts = smolq.levpoints[isummand];
-  int i = smolq.dimen()-1;
-  jseq[i]++;
-  while (i >= 0 && jseq[i] == levpts[i])
-    {
-      jseq[i] = 0;
-      i--;
-      if (i >= 0)
-        jseq[i]++;
-    }
-  sig.signalAfter(std::max(i, 0));
-
-  if (i < 0)
-    isummand++;
-
-  if (isummand < smolq.numSummands())
-    setPointAndWeight();
-
-  return *this;
-}
-
-/* Here we set the point coordinates according to |jseq| and
-   |isummand|. Also the weight is set here. */
-
-void
-smolpit::setPointAndWeight()
-{
-  // todo: raise if |smolq==NULL| or |jseq==NULL| or |sig==NULL| or
-  // |p==NULL| or |isummand>=smolq.numSummands()|
-  int l = smolq.level;
-  int d = smolq.dimen();
-  int sumk = (smolq.levels[isummand]).sum();
-  int m1exp = l + d - sumk - 1;
-  w = (2*(m1exp/2) == m1exp) ? 1.0 : -1.0;
-  w *= smolq.psc.noverk(d-1, sumk-l);
-  for (int i = 0; i < d; i++)
-    {
-      int ki = (smolq.levels[isummand])[i];
-      p[i] = (smolq.uquad).point(ki, jseq[i]);
-      w *= (smolq.uquad).weight(ki, jseq[i]);
-    }
-}
-
-/* Debug print. */
-void
-smolpit::print() const
-{
-  auto ff = std::cout.flags();
-  std::cout << "isum=" << std::left << std::setw(3) << isummand << std::right << ": [";
-  for (int i = 0; i < smolq.dimen(); i++)
-    std::cout << std::setw(2) << (smolq.levels[isummand])[i] << ' ';
-  std::cout << "] j=[";
-  for (int i = 0; i < smolq.dimen(); i++)
-    std::cout << std::setw(2) << jseq[i] << ' ';
-  std::cout << std::showpos << std::fixed << std::setprecision(3)
-            << "] " << std::setw(4) << w << "*(";
-  for (int i = 0; i < smolq.dimen()-1; i++)
-    std::cout << std::setw(4) << p[i] << ' ';
-  std::cout << std::setw(4) << p[smolq.dimen()-1] << ')' << std::endl;
-  std::cout.flags(ff);
-}
-
-/* Here is the constructor of |SmolyakQuadrature|. We have to setup
-   |levels|, |levpoints| and |cumevals|. We have to go through all
-   $d$-dimensional sequences $k$, such that $l\leq \vert k\vert\leq
-   l+d-1$ and all $k_i$ are positive integers. This is equivalent to
-   going through all $k$ such that $l-d\leq\vert k\vert\leq l-1$ and all
-   $k_i$ are non-negative integers. This is equivalent to going through
-   $d+1$ dimensional sequences $(k,x)$ such that $\vert(k,x)\vert =l-1$
-   and $x=0,\ldots,d-1$. The resulting sequence of positive integers is
-   obtained by adding $1$ to all $k_i$. */
-
-SmolyakQuadrature::SmolyakQuadrature(int d, int l, const OneDQuadrature &uq)
-  : QuadratureImpl<smolpit>(d), level(l), uquad(uq), psc(d-1, d-1)
-{
-  // todo: check |l>1|, |l>=d|
-  // todo: check |l>=uquad.miLevel()|, |l<=uquad.maxLevel()|
-  int cum = 0;
-  SymmetrySet ss(l-1, d+1);
-  for (symiterator si(ss); !si.isEnd(); ++si)
-    {
-      if ((*si)[d] <= d-1)
-        {
-          IntSequence lev((const IntSequence &)*si, 0, d);
-          lev.add(1);
-          levels.push_back(lev);
-          IntSequence levpts(d);
-          for (int i = 0; i < d; i++)
-            levpts[i] = uquad.numPoints(lev[i]);
-          levpoints.push_back(levpts);
-          cum += levpts.mult();
-          cumevals.push_back(cum);
-        }
-    }
-}
-
-/* Here we return a number of evalutions of the quadrature for the
-   given level. If the given level is the current one, we simply return
-   the maximum cumulative number of evaluations. Otherwise we call costly
-   |calcNumEvaluations| method. */
-
-int
-SmolyakQuadrature::numEvals(int l) const
-{
-  if (l != level)
-    return calcNumEvaluations(l);
-  else
-    return cumevals[numSummands()-1];
-}
-
-/* This divides all the evaluations to |tn| approximately equal groups,
-   and returns the beginning of the specified group |ti|. The granularity
-   of divisions are summands as listed by |levels|. */
-
-smolpit
-SmolyakQuadrature::begin(int ti, int tn, int l) const
-{
-  // todo: raise is |level!=l|
-  if (ti == tn)
-    return smolpit(*this, numSummands());
-
-  int totevals = cumevals[numSummands()-1];
-  int evals = (totevals*ti)/tn;
-  unsigned int isum = 0;
-  while (isum+1 < numSummands() && cumevals[isum+1] < evals)
-    isum++;
-  return smolpit(*this, isum);
-}
-
-/* This is the same in a structure as |@<|SmolyakQuadrature| constructor@>|.
-   We have to go through all summands and calculate
-   a number of evaluations in each summand. */
-
-int
-SmolyakQuadrature::calcNumEvaluations(int lev) const
-{
-  int cum = 0;
-  SymmetrySet ss(lev-1, dim+1);
-  for (symiterator si(ss); !si.isEnd(); ++si)
-    {
-      if ((*si)[dim] <= dim-1)
-        {
-          IntSequence lev((const IntSequence &)*si, 0, dim);
-          lev.add(1);
-          IntSequence levpts(dim);
-          for (int i = 0; i < dim; i++)
-            levpts[i] = uquad.numPoints(lev[i]);
-          cum += levpts.mult();
-        }
-    }
-  return cum;
-}
-
-/* This returns a maximum level such that the number of evaluations is
-   less than the given number. */
-
-void
-SmolyakQuadrature::designLevelForEvals(int max_evals, int &lev, int &evals) const
-{
-  int last_evals;
-  evals = 1;
-  lev = 1;
-  do
-    {
-      lev++;
-      last_evals = evals;
-      evals = calcNumEvaluations(lev);
-    }
-  while (lev < uquad.numLevels() && evals <= max_evals);
-  lev--;
-  evals = last_evals;
-}

dynare++/integ/cc/smolyak.hh

-// Copyright 2005, Ondra Kamenik
-
-// Smolyak quadrature.
-
-/* This file defines Smolyak (sparse grid) multidimensional quadrature
-   for non-nested underlying one dimensional quadrature. Let $Q^1_l$ denote
-   the one dimensional quadrature of $l$ level. Let $n_l$ denote a
-   number of points in the $l$ level. Than the Smolyak quadrature can be
-   defined as
-   $$Q^df=\sum_{l\leq\vert k\vert\leq l+d-1}(-1)^{l+d-\vert k\vert-1}\left(\matrix{d-1\cr
-   \vert k\vert-l}\right)(Q_{k_1}^1\otimes\ldots\otimes Q_{k_d}^1)f,$$
-   where $d$ is the dimension, $k$ is $d$-dimensional sequence of
-   integers, and $\vert k\vert$ denotes a sum of the sequence.
-
-   Here we define |smolpit| as Smolyak iterator and |SmolyakQuadrature|. */
-
-#ifndef SMOLYAK_H
-#define SMOLYAK_H
-
-#include "int_sequence.hh"
-#include "tl_static.hh"
-#include "vector_function.hh"
-#include "quadrature.hh"
-
-/* Here we define the Smolyak point iterator. The Smolyak formula can
-   be broken to a sum of product quadratures with various combinations of
-   levels. The iterator follows this pattern. It maintains an index to a
-   summand and then a point coordinates within the summand (product
-   quadrature). The array of summands to which the |isummand| points is
-   maintained by the |SmolyakQuadrature| class to which the object knows
-   the pointer |smolq|.
-
-   We provide a constructor which points to the beginning of the given
-   summand. This constructor is used in |SmolyakQuadrature::begin| method
-   which approximately divideds all the iterators to subsets of equal
-   size. */
-
-class SmolyakQuadrature;
-
-class smolpit
-{
-protected:
-  const SmolyakQuadrature &smolq;
-  unsigned int isummand{0};
-  IntSequence jseq;
-  ParameterSignal sig;
-  Vector p;
-  double w;
-public:
-  smolpit(const SmolyakQuadrature &q, unsigned int isum);
-  smolpit(const smolpit &spit) = default;
-  ~smolpit() = default;
-  bool operator==(const smolpit &spit) const;
-  bool
-  operator!=(const smolpit &spit) const
-  {
-    return !operator==(spit);
-  }
-  smolpit &operator=(const smolpit &spit) = delete;
-  smolpit &operator++();
-  const ParameterSignal &
-  signal() const
-  {
-    return sig;
-  }
-  const Vector &
-  point() const
-  {
-    return p;
-  }
-  double
-  weight() const
-  {
-    return w;
-  }
-  void print() const;
-protected:
-  void setPointAndWeight();
-};
-
-/* Here we define the class |SmolyakQuadrature|. It maintains an array
-   of summands of the Smolyak quadrature formula:
-   $$\sum_{l\leq\vert k\vert\leq l+d-1}(-1)^{l+d-\vert
-   k\vert-1}\left(\matrix{d-1\cr
-   \vert k\vert-l}\right)(Q_{k_1}^1\otimes\ldots\otimes Q_{k_d}^1)f$$
-   Each summand is fully specified by sequence $k$. The summands are here
-   represented (besides $k$) also by sequence of number of points in each
-   level selected by $k$, and also by a cummulative number of
-   evaluations. The latter two are added only for conveniency.
-
-   The summands in the code are given by |levels|, which is a vector of
-   $k$ sequences, further by |levpoints| which is a vector of sequences
-   of nuber of points in each level, and by |cumevals| which is the
-   cumulative number of points, this is $\sum_k\prod_{i=1}^dn_{k_i}$,
-   where the sum is done through all $k$ before the current.
-
-   The |levels| and |levpoints| vectors are used by |smolpit|. */
-
-class SmolyakQuadrature : public QuadratureImpl<smolpit>
-{
-  friend class smolpit;
-  int level;
-  const OneDQuadrature &uquad;
-  std::vector<IntSequence> levels;
-  std::vector<IntSequence> levpoints;
-  std::vector<int> cumevals;
-  PascalTriangle psc;
-public:
-  SmolyakQuadrature(int d, int l, const OneDQuadrature &uq);
-  ~SmolyakQuadrature() override = default;
-  int numEvals(int level) const override;
-  void designLevelForEvals(int max_eval, int &lev, int &evals) const;
-protected:
-  smolpit begin(int ti, int tn, int level) const override;
-  unsigned int
-  numSummands() const
-  {
-    return levels.size();
-  }
-private:
-  int calcNumEvaluations(int level) const;
-};
-
-#endif

dynare++/integ/cc/vector_function.cc

-// Copyright 2005, Ondra Kamenik
-
-#include "vector_function.hh"
-
-#include <dynlapack.h>
-
-#include <cmath>
-#include <algorithm>
-
-/* Just an easy constructor of sequence of booleans defaulting to
-   change everywhere. */
-
-ParameterSignal::ParameterSignal(int n)
-  : data(n, true)
-{
-}
-
-/* This sets |false| (no change) before a given parameter, and |true|
-   (change) after the given parameter (including). */
-
-void
-ParameterSignal::signalAfter(int l)
-{
-  for (size_t i = 0; i < std::min((size_t) l, data.size()); i++)
-    data[i] = false;
-  for (size_t i = l; i < data.size(); i++)
-    data[i] = true;
-}
-
-/* This constructs a function set hardcopying also the first. */
-VectorFunctionSet::VectorFunctionSet(const VectorFunction &f, int n)
-  : funcs(n)
-{
-  for (int i = 0; i < n; i++)
-    {
-      func_copies.push_back(f.clone());
-      funcs[i] = func_copies.back().get();
-    }
-}
-
-/* This constructs a function set with shallow copy in the first and
-   hard copies in others. */
-VectorFunctionSet::VectorFunctionSet(VectorFunction &f, int n)
-  : funcs(n)
-{
-  if (n > 0)
-    funcs[0] = &f;
-  for (int i = 1; i < n; i++)
-    {
-      func_copies.push_back(f.clone());
-      funcs[i] = func_copies.back().get();
-    }
-}
-
-
-/* Here we construct the object from the given function $f$ and given
-   variance-covariance matrix $\Sigma=$|vcov|. The matrix $A$ is
-   calculated as lower triangular and yields $\Sigma=AA^T$. */
-
-GaussConverterFunction::GaussConverterFunction(VectorFunction &f, const GeneralMatrix &vcov)
-  : VectorFunction(f), func(&f), A(vcov.numRows(), vcov.numRows()),
-    multiplier(calcMultiplier())
-{
-  // todo: raise if |A.numRows() != indim()|
-  calcCholeskyFactor(vcov);
-}
-
-/* Here we construct the object in the same way, however we mark the
-   function as to be deleted. */
-GaussConverterFunction::GaussConverterFunction(std::unique_ptr<VectorFunction> f, const GeneralMatrix &vcov)
-  : VectorFunction(*f), func_storage{move(f)}, func{func_storage.get()}, A(vcov.numRows(), vcov.numRows()),
-    multiplier(calcMultiplier())
-{
-  // todo: raise if |A.numRows() != indim()|
-  calcCholeskyFactor(vcov);
-}
-
-GaussConverterFunction::GaussConverterFunction(const GaussConverterFunction &f)
-  : VectorFunction(f), func_storage{f.func->clone()}, func{func_storage.get()}, A(f.A),
-    multiplier(f.multiplier)
-{
-}
-
-/* Here we evaluate the function
-   $g(y)={1\over\sqrt{\pi^n}}f\left(\sqrt{2}Ay\right)$. Since the matrix $A$ is lower
-   triangular, the change signal for the function $f$ will look like
-   $(0,\ldots,0,1,\ldots,1)$ where the first $1$ is in the same position
-   as the first change in the given signal |sig| of the input
-   $y=$|point|. */
-
-void
-GaussConverterFunction::eval(const Vector &point, const ParameterSignal &sig, Vector &out)
-{
-  ParameterSignal s(sig);
-  int i = 0;
-  while (i < indim() && !sig[i])
-    i++;
-  s.signalAfter(i);
-
-  Vector x(indim());
-  x.zeros();
-  A.multaVec(x, point);
-  x.mult(sqrt(2.0));
-
-  func->eval(x, s, out);
-
-  out.mult(multiplier);
-}
-
-/* This returns $1\over\sqrt{\pi^n}$. */
-
-double
-GaussConverterFunction::calcMultiplier() const
-{
-  return sqrt(pow(M_PI, -1*indim()));
-}
-
-void
-GaussConverterFunction::calcCholeskyFactor(const GeneralMatrix &vcov)
-{
-  A = vcov;
-
-  lapack_int rows = A.numRows(), lda = A.getLD();
-  for (int i = 0; i < rows; i++)
-    for (int j = i+1; j < rows; j++)
-      A.get(i, j) = 0.0;
-
-  lapack_int info;
-  dpotrf("L", &rows, A.base(), &lda, &info);
-  // todo: raise if |info!=1|
-}

dynare++/integ/cc/vector_function.hh

-// Copyright 2005, Ondra Kamenik
-
-// Vector function.
-
-/* This file defines interface for functions taking a vector as an input
-   and returning a vector (with a different size) as an output. We are
-   also introducing a parameter signalling; it is a boolean vector which
-   tracks parameters which were changed from the previous call. The
-   |VectorFunction| implementation can exploit this information and
-   evaluate the function more efficiently. The information can be
-   completely ignored.
-
-   From the signalling reason, and from other reasons, the function
-   evaluation is not |const|. */
-
-#ifndef VECTOR_FUNCTION_H
-#define VECTOR_FUNCTION_H
-
-#include "Vector.hh"
-#include "GeneralMatrix.hh"
-
-#include <vector>
-#include <memory>
-
-/* This is a simple class representing a vector of booleans. The items
-   night be retrieved or changed, or can be set |true| after some
-   point. This is useful when we multiply the vector with lower
-   triangular matrix.
-
-   |true| means that a parameter was changed. */
-
-class ParameterSignal
-{
-protected:
-  std::vector<bool> data;
-public:
-  ParameterSignal(int n);
-  ParameterSignal(const ParameterSignal &sig) = default;
-  ~ParameterSignal() = default;
-  void signalAfter(int l);
-  bool
-  operator[](int i) const
-  {
-    return data[i];
-  }
-  std::vector<bool>::reference
-  operator[](int i)
-  {
-    return data[i];
-  }
-};
-
-/* This is the abstract class for vector function. At this level of
-   abstraction we only need to know size of input vector and a size of
-   output vector.
-
-   The important thing here is a clone method, we will need to make hard
-   copies of vector functions since the evaluations are not |const|. The
-   hardcopies apply for parallelization. */
-
-class VectorFunction
-{
-protected:
-  int in_dim;
-  int out_dim;
-public:
-  VectorFunction(int idim, int odim)
-    : in_dim(idim), out_dim(odim)
-  {
-  }
-  VectorFunction(const VectorFunction &func) = default;
-  virtual ~VectorFunction() = default;
-  virtual std::unique_ptr<VectorFunction> clone() const = 0;
-  virtual void eval(const Vector &point, const ParameterSignal &sig, Vector &out) = 0;
-  int
-  indim() const
-  {
-    return in_dim;
-  }
-  int
-  outdim() const
-  {
-    return out_dim;
-  }
-};
-
-/* This makes |n| copies of |VectorFunction|. The first constructor
-   make exactly |n| new copies, the second constructor copies only the
-   pointer to the first and others are hard (real) copies.
-
-   The class is useful for making a given number of copies at once, and
-   this set can be reused many times if we need mupliple copis of the
-   function (for example for paralelizing the code). */
-
-class VectorFunctionSet
-{
-private:
-  // Stores the hard copies made by the class
-  std::vector<std::unique_ptr<VectorFunction>> func_copies;
-protected:
-  std::vector<VectorFunction *> funcs;
-public:
-  VectorFunctionSet(const VectorFunction &f, int n);
-  VectorFunctionSet(VectorFunction &f, int n);
-  ~VectorFunctionSet() = default;
-  VectorFunction &
-  getFunc(int i)
-  {
-    return *(funcs[i]);
-  }
-  int
-  getNum() const
-  {
-    return funcs.size();
-  }
-};
-
-/* This class wraps another |VectorFunction| to allow integration of a
-   function through normally distributed inputs. Namely, if one wants to
-   integrate
-   $${1\over\sqrt{(2\pi)^n\vert\Sigma\vert}}\int f(x)e^{-{1\over2}x^T\Sigma^{-1}x}{\rm d}x$$
-   then if we write $\Sigma=AA^T$ and $x=\sqrt{2}Ay$, we get integral
-   $${1\over\sqrt{(2\pi)^n\vert\Sigma\vert}}
-   \int f\left(\sqrt{2}Ay\right)e^{-y^Ty}\sqrt{2^n}\vert A\vert{\rm d}y=
-   {1\over\sqrt{\pi^n}}\int f\left(\sqrt{2}Ay\right)e^{-y^Ty}{\rm d}y,$$
-   which means that a given function $f$ we have to wrap to yield a function
-   $$g(y)={1\over\sqrt{\pi^n}}f\left(\sqrt{2}Ay\right).$$
-   This is exactly what this class is doing. This transformation is
-   useful since the Gauss--Hermite points and weights are defined for
-   weighting function $e^{-y^2}$, so this transformation allows using
-   Gauss--Hermite quadratures seemlessly in a context of integration through
-   normally distributed inputs.
-
-   The class maintains a pointer to the function $f$. When the object is
-   constructed by the first constructor, the $f$ is assumed to be owned by the
-   caller. If the object of this class is copied, then $f$ is copied and hence
-   stored in a std::unique_ptr. The second constructor takes a smart pointer to
-   the function and in that case the class takes ownership of $f$. */
-
-class GaussConverterFunction : public VectorFunction
-{
-private:
-  std::unique_ptr<VectorFunction> func_storage;
-protected:
-  VectorFunction *func;
-  GeneralMatrix A;
-  double multiplier;
-public:
-  GaussConverterFunction(VectorFunction &f, const GeneralMatrix &vcov);
-  GaussConverterFunction(std::unique_ptr<VectorFunction> f, const GeneralMatrix &vcov);
-  GaussConverterFunction(const GaussConverterFunction &f);
-  ~GaussConverterFunction() override = default;
-  std::unique_ptr<VectorFunction>
-  clone() const override
-  {
-    return std::make_unique<GaussConverterFunction>(*this);
-  }
-  void eval(const Vector &point, const ParameterSignal &sig, Vector &out) override;
-private:
-  double calcMultiplier() const;
-  void calcCholeskyFactor(const GeneralMatrix &vcov);
-};
-
-#endif

dynare++/integ/src/Makefile.am

-noinst_PROGRAMS = quadrature-points
-
-quadrature_points_SOURCES = quadrature-points.cc
-quadrature_points_CPPFLAGS = -I../.. -I../../sylv/cc -I../../integ/cc -I../../tl/cc
-quadrature_points_CXXFLAGS = $(AM_CXXFLAGS) $(THREAD_CXXFLAGS)
-quadrature_points_LDADD = ../cc/libinteg.a ../../tl/cc/libtl.a ../../parser/cc/libparser.a ../../sylv/cc/libsylv.a ../../utils/cc/libutils.a $(LAPACK_LIBS) $(BLAS_LIBS) $(LIBS) $(FLIBS)

dynare++/integ/src/quadrature-points.cc

-// Copyright (C) 2008-2011, Ondra Kamenik
-
-#include "parser/cc/matrix_parser.hh"
-#include "utils/cc/exception.hh"
-#include "sylv/cc/GeneralMatrix.hh"
-#include "sylv/cc/Vector.hh"
-#include "sylv/cc/SymSchurDecomp.hh"
-#include "sylv/cc/SylvException.hh"
-#include "integ/cc/quadrature.hh"
-#include "integ/cc/smolyak.hh"
-#include "integ/cc/product.hh"
-
-#include <getopt.h>
-
-#include <cmath>
-#include <cstdlib>
-#include <iostream>
-#include <fstream>
-#include <sstream>
-#include <memory>
-#include <string>
-
-struct QuadParams
-{
-  std::string outname;
-  std::string vcovname;
-  int max_level{3};
-  double discard_weight{0.0};
-  QuadParams(int argc, char **argv);
-  void check_consistency() const;
-private:
-  enum {opt_max_level, opt_discard_weight, opt_vcov};
-};
-
-QuadParams::QuadParams(int argc, char **argv)
-{
-  if (argc == 1)
-    {
-      // print the help and exit
-      std::cerr << "Usage: " << argv[0] << " [--max-level INTEGER] [--discard-weight FLOAT] [--vcov FILENAME] OUTPUT_FILENAME" << std::endl;
-      std::exit(EXIT_FAILURE);
-    }
-
-  outname = argv[argc-1];
-  argc--;
-
-  struct option const opts [] = {
-    {"max-level", required_argument, nullptr, opt_max_level},
-    {"discard-weight", required_argument, nullptr, opt_discard_weight},
-    {"vcov", required_argument, nullptr, opt_vcov},
-    {nullptr, 0, nullptr, 0}
-  };
-
-  int ret;
-  int index;
-  while (-1 != (ret = getopt_long(argc, argv, "", opts, &index)))
-    {
-      switch (ret)
-        {
-        case opt_max_level:
-          try
-            {
-              max_level = std::stoi(std::string{optarg});
-            }
-          catch (const std::invalid_argument &e)
-            {
-              std::cerr << "Couldn't parse integer " << optarg << ", ignored" << std::endl;
-            }
-          break;
-        case opt_discard_weight:
-          try
-            {
-              discard_weight = std::stod(std::string{optarg});
-            }
-          catch (const std::invalid_argument &e)
-            {
-              std::cerr << "Couldn't parse float " << optarg << ", ignored" << std::endl;
-            }
-          break;
-        case opt_vcov:
-          vcovname = optarg;
-          break;
-        }
-    }
-
-  check_consistency();
-}
-
-void
-QuadParams::check_consistency() const
-{
-  if (outname.empty())
-    {
-      std::cerr << "Error: output name not set" << std::endl;
-      std::exit(EXIT_FAILURE);
-    }
-
-  if (vcovname.empty())
-    {
-      std::cerr << "Error: vcov file name not set" << std::endl;
-      std::exit(EXIT_FAILURE);
-    }
-}
-
-int
-main(int argc, char **argv)
-{
-  QuadParams params(argc, argv);
-
-  // open output file for writing
-  std::ofstream fout{params.outname, std::ios::out | std::ios::trunc};
-  if (fout.fail())
-    {
-      std::cerr << "Could not open " << params.outname << " for writing" << std::endl;
-      std::exit(EXIT_FAILURE);
-    }
-
-  try
-    {
-      std::ifstream f{params.vcovname};
-      std::ostringstream buffer;
-      buffer << f.rdbuf();
-      std::string contents{buffer.str()};
-
-      // parse the vcov matrix
-      ogp::MatrixParser mp;
-      mp.parse(contents.length(), contents.c_str());
-      if (mp.nrows() != mp.ncols())
-        throw ogu::Exception(__FILE__, __LINE__,
-                             "VCOV matrix not square");
-      // and put to the GeneralMatrix
-      GeneralMatrix vcov(mp.nrows(), mp.ncols());
-      vcov.zeros();
-      for (ogp::MPIterator it = mp.begin(); it != mp.end(); ++it)
-        vcov.get(it.row(), it.col()) = *it;
-
-      // calculate the factor A of vcov, so that A*A^T=VCOV
-      GeneralMatrix A(vcov.numRows(), vcov.numRows());
-      SymSchurDecomp ssd(vcov);
-      ssd.getFactor(A);
-
-      // construct Gauss-Hermite quadrature
-      GaussHermite ghq;
-      // construct Smolyak quadrature
-      int level = params.max_level;
-      SmolyakQuadrature sq(vcov.numRows(), level, ghq);
-
-      std::cout << "Dimension:                " << vcov.numRows() << std::endl
-                << "Maximum level:            " << level << std::endl
-                << "Total number of nodes:    " << sq.numEvals(level) << std::endl;
-
-      // put the points to the vector
-      std::vector<std::unique_ptr<Vector>> points;
-      for (smolpit qit = sq.start(level); qit != sq.end(level); ++qit)
-        points.push_back(std::make_unique<Vector>((const Vector &) qit.point()));
-      // sort and uniq
-      std::sort(points.begin(), points.end(), [](auto &a, auto &b) { return a.get() < b.get(); });
-      auto new_end = std::unique(points.begin(), points.end());
-      points.erase(new_end, points.end());
-
-      std::cout << "Duplicit nodes removed:   " << (unsigned long) (sq.numEvals(level)-points.size())
-                << std::endl;
-
-      // calculate weights and mass
-      double mass = 0.0;
-      std::vector<double> weights;
-      for (auto & point : points)
-        {
-          weights.push_back(std::exp(-point->dot(*point)));
-          mass += weights.back();
-        }
-
-      // calculate discarded mass
-      double discard_mass = 0.0;
-      for (double weight : weights)
-        if (weight/mass < params.discard_weight)
-          discard_mass += weight;
-
-      std::cout << "Total mass discarded:     " << std::fixed << discard_mass/mass << std::endl;
-
-      // dump the results
-      int npoints = 0;
-      double upscale_weight = 1/(mass-discard_mass);
-      Vector x(vcov.numRows());
-      fout << std::setprecision(16);
-      for (int i = 0; i < (int) weights.size(); i++)
-        if (weights[i]/mass >= params.discard_weight)
-          {
-            // print the upscaled weight
-            fout << std::setw(20) << upscale_weight*weights[i];
-            // multiply point with the factor A and sqrt(2)
-            A.multVec(0.0, x, std::sqrt(2.), *(points[i]));
-            // print the coordinates
-            for (int j = 0; j < x.length(); j++)
-              fout << ' ' << std::setw(20) << x[j];
-            fout << std::endl;
-            npoints++;
-          }
-
-      std::cout << "Final number of points:   " << npoints << std::endl;
-
-      fout.close();
-    }
-  catch (const SylvException &e)
-    {
-      e.printMessage();
-      return EXIT_FAILURE;
-    }
-  catch (const ogu::Exception &e)
-    {
-      e.print();
-      return EXIT_FAILURE;
-    }
-
-  return EXIT_SUCCESS;
-}

dynare++/integ/testing/Makefile.am

-check_PROGRAMS = tests
-
-tests_SOURCES = tests.cc
-tests_CPPFLAGS = -I../cc -I../../tl/cc -I../../sylv/cc -I$(top_srcdir)/mex/sources
-tests_CXXFLAGS = $(AM_CXXFLAGS) $(THREAD_CXXFLAGS)
-tests_LDFLAGS = $(AM_LDFLAGS) $(LDFLAGS_MATIO)
-tests_LDADD = ../../tl/cc/libtl.a ../../sylv/cc/libsylv.a ../cc/libinteg.a $(LAPACK_LIBS) $(BLAS_LIBS) $(LIBS) $(FLIBS) $(LIBADD_MATIO)
-
-check-local:
-	./tests

dynare++/integ/testing/tests.cc

-/* $Id: tests.cpp 431 2005-08-16 15:41:01Z kamenik $ */
-/* Copyright 2005, Ondra Kamenik */
-
-#include "GeneralMatrix.hh"
-#include <dynlapack.h>
-#include "SylvException.hh"
-
-#include "rfs_tensor.hh"
-#include "normal_moments.hh"
-
-#include "vector_function.hh"
-#include "quadrature.hh"
-#include "smolyak.hh"
-#include "product.hh"
-#include "quasi_mcarlo.hh"
-
-#include <iomanip>
-#include <chrono>
-#include <cmath>
-#include <iostream>
-#include <utility>
-#include <array>
-#include <memory>
-#include <cstdlib>
-
-// evaluates unfolded (Dx)^k power, where x is a vector, D is a
-// Cholesky factor (lower triangular)
-class MomentFunction : public VectorFunction
-{
-  GeneralMatrix D;
-  int k;
-public:
-  MomentFunction(const GeneralMatrix &inD, int kk)
-    : VectorFunction(inD.numRows(), UFSTensor::calcMaxOffset(inD.numRows(), kk)),
-      D(inD), k(kk)
-  {
-  }
-  MomentFunction(const MomentFunction &func) = default;
-
-  std::unique_ptr<VectorFunction>
-  clone() const override
-  {
-    return std::make_unique<MomentFunction>(*this);
-  }
-  void eval(const Vector &point, const ParameterSignal &sig, Vector &out) override;
-};
-
-void
-MomentFunction::eval(const Vector &point, const ParameterSignal &sig, Vector &out)
-{
-  if (point.length() != indim() || out.length() != outdim())
-    {
-      std::cerr << "Wrong length of vectors in MomentFunction::eval" << std::endl;
-      std::exit(EXIT_FAILURE);
-    }
-  Vector y(point);
-  y.zeros();
-  D.multaVec(y, point);
-  URSingleTensor ypow(y, k);
-  out.zeros();
-  out.add(1.0, ypow.getData());
-}
-
-class TensorPower : public VectorFunction
-{
-  int k;
-public:
-  TensorPower(int nvar, int kk)
-    : VectorFunction(nvar, UFSTensor::calcMaxOffset(nvar, kk)), k(kk)
-  {
-  }
-  TensorPower(const TensorPower &func) = default;
-
-  std::unique_ptr<VectorFunction>
-  clone() const override
-  {
-    return std::make_unique<TensorPower>(*this);
-  }
-  void eval(const Vector &point, const ParameterSignal &sig, Vector &out) override;
-};
-
-void
-TensorPower::eval(const Vector &point, const ParameterSignal &sig, Vector &out)
-{
-  if (point.length() != indim() || out.length() != outdim())
-    {
-      std::cerr << "Wrong length of vectors in TensorPower::eval" << std::endl;
-      std::exit(EXIT_FAILURE);
-    }
-  URSingleTensor ypow(point, k);
-  out.zeros();
-  out.add(1.0, ypow.getData());
-}
-
-// evaluates (1+1/d)^d*(x_1*...*x_d)^(1/d), its integral over <0,1>^d
-// is 1.0, and its variation grows exponetially
-// d = dim
-class Function1 : public VectorFunction
-{
-  int dim;
-public:
-  Function1(int d)
-    : VectorFunction(d, 1), dim(d)
-  {
-  }
-  Function1(const Function1 &f)
-    : VectorFunction(f.indim(), f.outdim()), dim(f.dim)
-  {
-  }
-  std::unique_ptr<VectorFunction>
-  clone() const override
-  {
-    return std::make_unique<Function1>(*this);
-  }
-  void eval(const Vector &point, const ParameterSignal &sig, Vector &out) override;
-};
-
-void
-Function1::eval(const Vector &point, const ParameterSignal &sig, Vector &out)
-{
-  if (point.length() != dim || out.length() != 1)
-    {
-      std::cerr << "Wrong length of vectors in Function1::eval" << std::endl;
-      std::exit(EXIT_FAILURE);
-    }
-  double r = 1;
-  for (int i = 0; i < dim; i++)
-    r *= point[i];
-  r = pow(r, 1.0/dim);
-  r *= pow(1.0 + 1.0/dim, (double) dim);
-  out[0] = r;
-}
-
-// evaluates Function1 but with transformation x_i=0.5(y_i+1)
-// this makes the new function integrate over <-1,1>^d to 1.0
-class Function1Trans : public Function1
-{
-public:
-  Function1Trans(int d)
-    : Function1(d)
-  {
-  }
-  Function1Trans(const Function1Trans &func) = default;
-
-  std::unique_ptr<VectorFunction>
-  clone() const override
-  {
-    return std::make_unique<Function1Trans>(*this);
-  }
-  void eval(const Vector &point, const ParameterSignal &sig, Vector &out) override;
-};
-
-void
-Function1Trans::eval(const Vector &point, const ParameterSignal &sig, Vector &out)
-{
-  Vector p(point.length());
-  for (int i = 0; i < p.length(); i++)
-    p[i] = 0.5*(point[i]+1);
-  Function1::eval(p, sig, out);
-  out.mult(pow(0.5, indim()));
-}
-
-// WallTimer class. Constructor saves the wall time, destructor
-// cancels the current time from the saved, and prints the message
-// with time information
-class WallTimer
-{
-  std::string mes;
-  std::chrono::time_point<std::chrono::high_resolution_clock> start;
-  bool new_line;
-public:
-  WallTimer(std::string m, bool nl = true)
-    : mes{m}, start{std::chrono::high_resolution_clock::now()}, new_line{nl}
-  {
-  }
-  ~WallTimer()
-  {
-    auto end = std::chrono::high_resolution_clock::now();
-    std::chrono::duration<double> duration = end - start;
-    std::cout << mes << std::setw(8) << std::setprecision(4) << duration.count();
-    if (new_line)
-      std::cout << std::endl;
-  }
-};
-
-/****************************************************/
-/*     declaration of TestRunnable class            */
-/****************************************************/
-class TestRunnable
-{
-public:
-  const std::string name;
-  int dim; // dimension of the solved problem
-  int nvar; // number of variable of the solved problem
-  TestRunnable(std::string name_arg, int d, int nv)
-    : name{move(name_arg)}, dim(d), nvar(nv)
-  {
-  }
-  virtual ~TestRunnable() = default;
-  bool test() const;
-  virtual bool run() const = 0;
-protected:
-  static bool smolyak_normal_moments(const GeneralMatrix &m, int imom, int level);
-  static bool product_normal_moments(const GeneralMatrix &m, int imom, int level);
-  static bool qmc_normal_moments(const GeneralMatrix &m, int imom, int level);
-  static bool smolyak_product_cube(const VectorFunction &func, const Vector &res,
-                                   double tol, int level);
-  static bool qmc_cube(const VectorFunction &func, double res, double tol, int level);
-};
-
-bool
-TestRunnable::test() const
-{
-  std::cout << "Running test <" << name << ">" << std::endl;
-  bool passed;
-  {
-    WallTimer tim("Wall clock time ", false);
-    passed = run();
-  }
-  if (passed)
-    {
-      std::cout << "............................ passed" << std::endl << std::endl;
-      return passed;
-    }
-  else
-    {
-      std::cout << "............................ FAILED" << std::endl << std::endl;
-      return passed;
-    }
-}
-
-/****************************************************/
-/*     definition of TestRunnable static methods    */
-/****************************************************/
-bool
-TestRunnable::smolyak_normal_moments(const GeneralMatrix &m, int imom, int level)
-{
-  // first make m*m' and then Cholesky factor
-  GeneralMatrix mtr(m, "transpose");
-  GeneralMatrix msq(m, mtr);
-
-  // make vector function
-  int dim = m.numRows();
-  TensorPower tp(dim, imom);
-  GaussConverterFunction func(tp, msq);
-
-  // smolyak quadrature
-  Vector smol_out(UFSTensor::calcMaxOffset(dim, imom));
-  {
-    WallTimer tim("\tSmolyak quadrature time:         ");
-    GaussHermite gs;
-    SmolyakQuadrature quad(dim, level, gs);
-    quad.integrate(func, level, sthread::detach_thread_group::max_parallel_threads, smol_out);
-    std::cout << "\tNumber of Smolyak evaluations:    " << quad.numEvals(level) << std::endl;
-  }
-
-  // check against theoretical moments
-  UNormalMoments moments(imom, msq);
-  smol_out.add(-1.0, (moments.get(Symmetry(imom)))->getData());
-  std::cout << "\tError:                         " << std::setw(16) << std::setprecision(12) << smol_out.getMax() << std::endl;
-  return smol_out.getMax() < 1.e-7;
-}
-
-bool
-TestRunnable::product_normal_moments(const GeneralMatrix &m, int imom, int level)
-{
-  // first make m*m' and then Cholesky factor
-  GeneralMatrix mtr(m, "transpose");
-  GeneralMatrix msq(m, mtr);
-
-  // make vector function
-  int dim = m.numRows();
-  TensorPower tp(dim, imom);
-  GaussConverterFunction func(tp, msq);
-
-  // product quadrature
-  Vector prod_out(UFSTensor::calcMaxOffset(dim, imom));
-  {
-    WallTimer tim("\tProduct quadrature time:         ");
-    GaussHermite gs;
-    ProductQuadrature quad(dim, gs);
-    quad.integrate(func, level, sthread::detach_thread_group::max_parallel_threads, prod_out);
-    std::cout << "\tNumber of product evaluations:    " << quad.numEvals(level) << std::endl;
-  }
-
-  // check against theoretical moments
-  UNormalMoments moments(imom, msq);
-  prod_out.add(-1.0, (moments.get(Symmetry(imom)))->getData());
-  std::cout << "\tError:                         " << std::setw(16) << std::setprecision(12) << prod_out.getMax() << std::endl;
-  return prod_out.getMax() < 1.e-7;
-}
-
-bool
-TestRunnable::smolyak_product_cube(const VectorFunction &func, const Vector &res,
-                                   double tol, int level)
-{
-  if (res.length() != func.outdim())
-    {
-      std::cerr << "Incompatible dimensions of check value and function." << std::endl;
-      std::exit(EXIT_FAILURE);
-    }
-
-  GaussLegendre glq;
-  Vector out(func.outdim());
-  double smol_error;
-  double prod_error;
-  {
-    WallTimer tim("\tSmolyak quadrature time:         ");
-    SmolyakQuadrature quad(func.indim(), level, glq);
-    quad.integrate(func, level, sthread::detach_thread_group::max_parallel_threads, out);
-    out.add(-1.0, res);
-    smol_error = out.getMax();
-    std::cout << "\tNumber of Smolyak evaluations:    " << quad.numEvals(level) << std::endl;
-    std::cout << "\tError:                            " << std::setw(16) << std::setprecision(12) << smol_error << std::endl;
-  }
-  {
-    WallTimer tim("\tProduct quadrature time:         ");
-    ProductQuadrature quad(func.indim(), glq);
-    quad.integrate(func, level, sthread::detach_thread_group::max_parallel_threads, out);
-    out.add(-1.0, res);
-    prod_error = out.getMax();
-    std::cout << "\tNumber of product evaluations:    " << quad.numEvals(level) << std::endl;
-    std::cout << "\tError:                            " << std::setw(16) << std::setprecision(12) << prod_error << std::endl;
-  }
-
-  return smol_error < tol && prod_error < tol;
-}
-
-bool
-TestRunnable::qmc_cube(const VectorFunction &func, double res, double tol, int level)
-{
-  Vector r(1);
-  double error1;
-  {
-    WallTimer tim("\tQuasi-Monte Carlo (Warnock scrambling) time:  ");
-    WarnockPerScheme wps;
-    QMCarloCubeQuadrature qmc(func.indim(), level, wps);
-    //		qmc.savePoints("warnock.txt", level);
-    qmc.integrate(func, level, sthread::detach_thread_group::max_parallel_threads, r);
-    error1 = std::max(res - r[0], r[0] - res);
-    std::cout << "\tQuasi-Monte Carlo (Warnock scrambling) error: " << std::setw(16) << std::setprecision(12) << error1 << std::endl;
-  }
-  double error2;
-  {
-    WallTimer tim("\tQuasi-Monte Carlo (reverse scrambling) time:  ");
-    ReversePerScheme rps;
-    QMCarloCubeQuadrature qmc(func.indim(), level, rps);
-    //		qmc.savePoints("reverse.txt", level);
-    qmc.integrate(func, level, sthread::detach_thread_group::max_parallel_threads, r);
-    error2 = std::max(res - r[0], r[0] - res);
-    std::cout << "\tQuasi-Monte Carlo (reverse scrambling) error: " << std::setw(16) << std::setprecision(12) << error2 << std::endl;
-  }
-  double error3;
-  {
-    WallTimer tim("\tQuasi-Monte Carlo (no scrambling) time:       ");
-    IdentityPerScheme ips;
-    QMCarloCubeQuadrature qmc(func.indim(), level, ips);
-    //		qmc.savePoints("identity.txt", level);
-    qmc.integrate(func, level, sthread::detach_thread_group::max_parallel_threads, r);
-    error3 = std::max(res - r[0], r[0] - res);
-    std::cout << "\tQuasi-Monte Carlo (no scrambling) error:      " << std::setw(16) << std::setprecision(12) << error3 << std::endl;
-  }
-
-  return error1 < tol && error2 < tol && error3 < tol;
-}
-
-/****************************************************/
-/*     definition of TestRunnable subclasses        */
-/****************************************************/
-class SmolyakNormalMom1 : public TestRunnable
-{
-public:
-  SmolyakNormalMom1()
-    : TestRunnable("Smolyak normal moments (dim=2, level=4, order=4)", 4, 2)
-  {
-  }
-
-  bool
-  run() const override
-  {
-    GeneralMatrix m(2, 2);
-    m.zeros(); m.get(0, 0) = 1; m.get(1, 1) = 1;
-    return smolyak_normal_moments(m, 4, 4);
-  }
-};
-
-class SmolyakNormalMom2 : public TestRunnable
-{
-public:
-  SmolyakNormalMom2()
-    : TestRunnable("Smolyak normal moments (dim=3, level=8, order=8)", 8, 3)
-  {
-  }
-
-  bool
-  run() const override
-  {
-    GeneralMatrix m(3, 3);
-    m.zeros();
-    m.get(0, 0) = 1; m.get(0, 2) = 0.5; m.get(1, 1) = 1;
-    m.get(1, 0) = 0.5; m.get(2, 2) = 2; m.get(2, 1) = 4;
-    return smolyak_normal_moments(m, 8, 8);
-  }
-};
-
-class ProductNormalMom1 : public TestRunnable
-{
-public:
-  ProductNormalMom1()
-    : TestRunnable("Product normal moments (dim=2, level=4, order=4)", 4, 2)
-  {
-  }
-
-  bool
-  run() const override
-  {
-    GeneralMatrix m(2, 2);
-    m.zeros(); m.get(0, 0) = 1; m.get(1, 1) = 1;
-    return product_normal_moments(m, 4, 4);
-  }
-};
-
-class ProductNormalMom2 : public TestRunnable
-{
-public:
-  ProductNormalMom2()
-    : TestRunnable("Product normal moments (dim=3, level=8, order=8)", 8, 3)
-  {
-  }
-
-  bool
-  run() const override
-  {
-    GeneralMatrix m(3, 3);
-    m.zeros();
-    m.get(0, 0) = 1; m.get(0, 2) = 0.5; m.get(1, 1) = 1;
-    m.get(1, 0) = 0.5; m.get(2, 2) = 2; m.get(2, 1) = 4;
-    return product_normal_moments(m, 8, 8);
-  }
-};
-
-// note that here we pass 1,1 to tls since smolyak has its own PascalTriangle
-class F1GaussLegendre : public TestRunnable
-{
-public:
-  F1GaussLegendre()
-    : TestRunnable("Function1 Gauss-Legendre (dim=6, level=13", 1, 1)
-  {
-  }
-
-  bool
-  run() const override
-  {
-    Function1Trans f1(6);
-    Vector res(1); res[0] = 1.0;
-    return smolyak_product_cube(f1, res, 1e-2, 13);
-  }
-};
-
-class F1QuasiMCarlo : public TestRunnable
-{
-public:
-  F1QuasiMCarlo()
-    : TestRunnable("Function1 Quasi-Monte Carlo (dim=6, level=1000000)", 1, 1)
-  {
-  }
-
-  bool
-  run() const override
-  {
-    Function1 f1(6);
-    return qmc_cube(f1, 1.0, 1.e-4, 1000000);
-  }
-};
-
-int
-main()
-{
-  std::vector<std::unique_ptr<TestRunnable>> all_tests;
-  // fill in vector of all tests
-  all_tests.push_back(std::make_unique<SmolyakNormalMom1>());
-  all_tests.push_back(std::make_unique<SmolyakNormalMom2>());
-  all_tests.push_back(std::make_unique<ProductNormalMom1>());
-  all_tests.push_back(std::make_unique<ProductNormalMom2>());
-  all_tests.push_back(std::make_unique<F1GaussLegendre>());
-  all_tests.push_back(std::make_unique<F1QuasiMCarlo>());
-
-  // find maximum dimension and maximum nvar
-  int dmax = 0;
-  int nvmax = 0;
-  for (const auto &test : all_tests)
-    {
-      if (dmax < test->dim)
-        dmax = test->dim;
-      if (nvmax < test->nvar)
-        nvmax = test->nvar;
-    }
-  tls.init(dmax, nvmax); // initialize library
-
-  // launch the tests
-  int success = 0;
-  for (const auto &test : all_tests)
-    {
-      try
-        {
-          if (test->test())
-            success++;
-        }
-      catch (const TLException &e)
-        {
-          std::cout << "Caught TL exception in <" << test->name << ">:" << std::endl;
-          e.print();
-        }
-      catch (SylvException &e)
-        {
-          std::cout << "Caught Sylv exception in <" << test->name << ">:" << std::endl;
-          e.printMessage();
-        }
-    }
-
-  int nfailed = all_tests.size() - success;
-  std::cout << "There were " << nfailed << " tests that failed out of "
-            << all_tests.size() << " tests run." << std::endl;
-
-  if (nfailed)
-    return EXIT_FAILURE;
-  else
-    return EXIT_SUCCESS;
-}

dynare++/kord/Makefile.am

-noinst_LIBRARIES = libkord.a
-
-libkord_a_SOURCES = \
-	approximation.cc \
-	approximation.hh \
-	decision_rule.cc \
-	decision_rule.hh \
-	dynamic_model.cc \
-	dynamic_model.hh \
-	faa_di_bruno.cc \
-	faa_di_bruno.hh \
-	first_order.cc \
-	first_order.hh \
-	global_check.cc \
-	global_check.hh \
-	kord_exception.hh \
-	korder.cc \
-	korder.hh \
-	korder_stoch.cc \
-	korder_stoch.hh \
-	journal.cc \
-	journal.hh \
-	mersenne_twister.hh \
-	normal_conjugate.cc \
-	normal_conjugate.hh \
-	random.cc \
-	random.hh
-
-libkord_a_CPPFLAGS = -I../sylv/cc -I../tl/cc -I../integ/cc -I$(top_srcdir)/mex/sources $(CPPFLAGS_MATIO)
-libkord_a_CXXFLAGS = $(AM_CXXFLAGS) $(THREAD_CXXFLAGS)
-
-check_PROGRAMS = tests
-
-tests_SOURCES = tests.cc
-tests_CPPFLAGS = -I../sylv/cc -I../tl/cc -I../integ/cc -I$(top_srcdir)/mex/sources
-tests_CXXFLAGS = $(AM_CXXFLAGS) $(THREAD_CXXFLAGS)
-tests_LDFLAGS = $(AM_LDFLAGS) $(LDFLAGS_MATIO)
-tests_LDADD = libkord.a ../tl/cc/libtl.a ../sylv/cc/libsylv.a $(LAPACK_LIBS) $(BLAS_LIBS) $(LIBS) $(FLIBS) $(LIBADD_MATIO)
-
-check-local:
-	./tests
-
-CLEANFILES = out.txt

dynare++/kord/approximation.hh

-// Copyright 2005, Ondra Kamenik
-
-// Approximating model solution
-
-/* The class |Approximation| in this file is a main interface to the
-   algorithms calculating approximations to the decision rule about
-   deterministic and stochastic steady states.
-
-   The approximation about a deterministic steady state is solved by
-   classes |@<|FirstOrder| class declaration@>| and |@<|KOrder| class
-   declaration@>|. The approximation about the stochastic steady state is
-   solved by class |@<|KOrderStoch| class declaration@>| together with a
-   method of |Approximation| class |@<|Approximation::walkStochSteady|
-   code@>|.
-
-   The approximation about the stochastic steady state is done with
-   explicit expression of forward derivatives of $g^{**}$. More formally,
-   we have to solve the decision rule $g$ from the implicit system:
-   $$E_t(f(g^{**}(g^*(y^*,u_t,\sigma),u_{t+1},\sigma),g(y^*,u_t,\sigma),y_t,u_t))=0$$
-   The term within the expectations can be Taylor expanded, and the
-   expectation can be driven into the formula. However, when doing this
-   at $\sigma\not=0$, the term $g^{**}$ at $\sigma\not=0$ is dependent on
-   $u_{t+1}$ and thus the integral of its approximation includes all
-   derivatives wrt. $u$ of $g^{**}$. Note that for $\sigma=0$, the
-   derivatives of $g^{**}$ in this context are constant. This is the main
-   difference between the approximation at deterministic steady
-   ($\sigma=0$), and stochastic steady ($\sigma\not=0$). This means that
-   $k$-order derivative of the above equation at $\sigma\not=0$ depends of
-   all derivatives of $g^**$ (including those with order greater than
-   $k$).
-
-   The explicit expression of the forward $g^{**}$ means that the
-   derivatives of $g$ are not solved simultaneously, but that the forward
-   derivatives of $g^{**}$ are calculated as an extrapolation based on
-   the approximation at lower $\sigma$. This is exactly what does the
-   |@<|Approximation::walkStochSteady| code@>|. It starts at the
-   deterministic steady state, and in a few steps it adds to $\sigma$
-   explicitly expressing forward $g^{**}$ from a previous step.
-
-   Further details on the both solution methods are given in (todo: put
-   references here when they exist).
-
-   Very important note: all classes here used for calculation of decision
-   rule approximation are folded. For the time being, it seems that faa
-   Di Bruno formula is quicker for folded tensors, and that is why we
-   stick to folded tensors here. However, when the calcs are done, we
-   calculate also its unfolded versions, to be available for simulations
-   and so on. */
-
-#ifndef APPROXIMATION_H
-#define APPROXIMATION_H
-
-#include "dynamic_model.hh"
-#include "decision_rule.hh"
-#include "korder.hh"
-#include "journal.hh"
-
-/* This class is used to calculate derivatives by faa Di Bruno of the
-   $$f(g^{**}(g^*(y^*,u,\sigma),u',\sigma),g(y^*,u,\sigma),y^*,u)$$ with
-   respect $u'$. In order to keep it as simple as possible, the class
-   represents an equivalent (with respect to $u'$) container for
-   $f(g^{**}(y^*,u',\sigma),0,0,0)$. The class is used only for
-   evaluation of approximation error in |Approximation| class, which is
-   calculated in |Approximation::calcStochShift| method.
-
-   Since it is a folded version, we inherit from
-   |StackContainer<FGSTensor>| and |FoldedStackContainer|. To construct
-   it, we need only the $g^{**}$ container and size of stacks. */
-
-class ZAuxContainer : public StackContainer<FGSTensor>, public FoldedStackContainer
-{
-public:
-  using _Ctype = StackContainer<FGSTensor>::_Ctype;
-  using itype = StackContainer<FGSTensor>::itype;
-  ZAuxContainer(const _Ctype *gss, int ngss, int ng, int ny, int nu);
-  itype getType(int i, const Symmetry &s) const override;
-};
-
-/* This class provides an interface to approximation algorithms. The
-   core method is |walkStochSteady| which calculates the approximation
-   about stochastic steady state in a given number of steps. The number
-   is given as a parameter |ns| of the constructor. If the number is
-   equal to zero, the resulted approximation is about the deterministic
-   steady state.
-
-   An object is constructed from the |DynamicModel|, and the number of
-   steps |ns|. Also, we pass a reference to journal. That's all. The
-   result of the core method |walkStochSteady| is a decision rule |dr|
-   and a matrix |ss| whose columns are steady states for increasing
-   $\sigma$ during the walk. Both can be retrived by public methods. The
-   first column of the matrix is the deterministic steady state, the last
-   is the stochastic steady state for the full size shocks.
-
-   The method |walkStochSteady| calls the following methods:
-   |approxAtSteady| calculates an initial approximation about the
-   deterministic steady, |saveRuleDerivs| saves derivatives of a rule for
-   the following step in |rule_ders| and |rule_ders_ss| (see
-   |@<|Approximation::saveRuleDerivs| code@>| for their description),
-   |check| reports an error of the current approximation and
-   |calcStochShift| (called from |check|) calculates a shift of the
-   system equations due to uncertainity.
-
-   dr\_centralize is a new option. dynare++ was automatically expressing
-   results around the fixed point instead of the deterministic steady
-   state. dr\_centralize controls this behavior. */
-
-class Approximation
-{
-  DynamicModel &model;
-  Journal &journal;
-  FGSContainer *rule_ders;
-  FGSContainer *rule_ders_ss;
-  FoldDecisionRule *fdr;
-  UnfoldDecisionRule *udr;
-  const PartitionY ypart;
-  const FNormalMoments mom;
-  IntSequence nvs;
-  int steps;
-  bool dr_centralize;
-  double qz_criterium;
-  TwoDMatrix ss;
-public:
-  Approximation(DynamicModel &m, Journal &j, int ns, bool dr_centr, double qz_crit);
-  virtual
-  ~Approximation();
-
-  const FoldDecisionRule&getFoldDecisionRule() const;
-  const UnfoldDecisionRule&getUnfoldDecisionRule() const;
-  const TwoDMatrix &
-  getSS() const
-  {
-    return ss;
-  }
-  const DynamicModel &
-  getModel() const
-  {
-    return model;
-  }
-
-  void walkStochSteady();
-  TwoDMatrix *calcYCov() const;
-  const FGSContainer *
-  get_rule_ders() const
-  {
-    return rule_ders;
-  }
-  const FGSContainer *
-  get_rule_ders_ss() const
-  {
-    return rule_ders;
-  }
-protected:
-  void approxAtSteady();
-  void calcStochShift(Vector &out, double at_sigma) const;
-  void saveRuleDerivs(const FGSContainer &g);
-  void check(double at_sigma) const;
-};
-
-#endif

dynare++/kord/decision_rule.cc

-// Copyright 2004, Ondra Kamenik
-
-#include "kord_exception.hh"
-#include "decision_rule.hh"
-#include "dynamic_model.hh"
-
-#include "SymSchurDecomp.hh"
-
-#include <dynlapack.h>
-
-#include <limits>
-#include <utility>
-#include <memory>
-
-template <>
-int DRFixPoint<KOrder::fold>::max_iter = 10000;
-template <>
-int DRFixPoint<KOrder::unfold>::max_iter = 10000;
-template <>
-double DRFixPoint<KOrder::fold>::tol = 1.e-10;
-template <>
-double DRFixPoint<KOrder::unfold>::tol = 1.e-10;
-template <>
-int DRFixPoint<KOrder::fold>::max_newton_iter = 50;
-template <>
-int DRFixPoint<KOrder::unfold>::max_newton_iter = 50;
-template <>
-int DRFixPoint<KOrder::fold>::newton_pause = 100;
-template <>
-int DRFixPoint<KOrder::unfold>::newton_pause = 100;
-
-// |FoldDecisionRule| conversion from |UnfoldDecisionRule|
-FoldDecisionRule::FoldDecisionRule(const UnfoldDecisionRule &udr)
-  : DecisionRuleImpl<KOrder::fold>(ctraits<KOrder::fold>::Tpol(udr.nrows(), udr.nvars()),
-                                   udr.ypart, udr.nu, udr.ysteady)
-{
-  for (const auto & it : udr)
-    {
-      insert(new ctraits<KOrder::fold>::Ttensym(*(it.second)));
-    }
-}
-
-// |UnfoldDecisionRule| conversion from |FoldDecisionRule|
-UnfoldDecisionRule::UnfoldDecisionRule(const FoldDecisionRule &fdr)
-  : DecisionRuleImpl<KOrder::unfold>(ctraits<KOrder::unfold>::Tpol(fdr.nrows(), fdr.nvars()),
-                                     fdr.ypart, fdr.nu, fdr.ysteady)
-{
-  for (const auto & it : fdr)
-    {
-      insert(new ctraits<KOrder::unfold>::Ttensym(*(it.second)));
-    }
-}
-
-SimResults::~SimResults()
-{
-  for (int i = 0; i < getNumSets(); i++)
-    {
-      delete data[i];
-      delete shocks[i];
-    }
-}
-
-/* This runs simulations with an output to journal file. Note that we
-   report how many simulations had to be thrown out due to Nan or Inf. */
-
-void
-SimResults::simulate(int num_sim, const DecisionRule &dr, const Vector &start,
-                     const TwoDMatrix &vcov, Journal &journal)
-{
-  JournalRecordPair paa(journal);
-  paa << "Performing " << num_sim << " stochastic simulations for "
-      << num_per << " periods burning " << num_burn << " initial periods"  << endrec;
-  simulate(num_sim, dr, start, vcov);
-  int thrown = num_sim - data.size();
-  if (thrown > 0)
-    {
-      JournalRecord rec(journal);
-      rec << "I had to throw " << thrown << " simulations away due to Nan or Inf" << endrec;
-    }
-}
-
-/* This runs a given number of simulations by creating
-   |SimulationWorker| for each simulation and inserting them to the
-   thread group. */
-
-void
-SimResults::simulate(int num_sim, const DecisionRule &dr, const Vector &start,
-                     const TwoDMatrix &vcov)
-{
-  std::vector<RandomShockRealization> rsrs;
-  rsrs.reserve(num_sim);
-
-  sthread::detach_thread_group gr;
-  for (int i = 0; i < num_sim; i++)
-    {
-      RandomShockRealization sr(vcov, system_random_generator.int_uniform());
-      rsrs.push_back(sr);
-      gr.insert(std::make_unique<SimulationWorker>(*this, dr, DecisionRule::horner,
-                                                   num_per+num_burn, start, rsrs.back()));
-    }
-  gr.run();
-}
-
-/* This adds the data with the realized shocks. It takes only periods
-   which are not to be burnt. If the data is not finite, the both data
-   and shocks are thrown away. */
-
-bool
-SimResults::addDataSet(TwoDMatrix *d, ExplicitShockRealization *sr, const ConstVector &st)
-{
-  KORD_RAISE_IF(d->nrows() != num_y,
-                "Incompatible number of rows for SimResults::addDataSets");
-  KORD_RAISE_IF(d->ncols() != num_per+num_burn,
-                "Incompatible number of cols for SimResults::addDataSets");
-  bool ret = false;
-  if (d->isFinite())
-    {
-      data.push_back(new TwoDMatrix((const TwoDMatrix &) (*d), num_burn, num_per));
-      shocks.push_back(new ExplicitShockRealization(
-                                                    ConstTwoDMatrix(sr->getShocks(), num_burn, num_per)));
-      if (num_burn == 0)
-        start.emplace_back(st);
-      else
-        start.emplace_back(d->getCol(num_burn-1));
-      ret = true;
-    }
-
-  delete d;
-  delete sr;
-  return ret;
-}
-
-void
-SimResults::writeMat(const char *base, const char *lname) const
-{
-  char matfile_name[100];
-  sprintf(matfile_name, "%s.mat", base);
-  mat_t *matfd = Mat_Create(matfile_name, nullptr);
-  if (matfd != nullptr)
-    {
-      writeMat(matfd, lname);
-      Mat_Close(matfd);
-    }
-}
-
-/* This save the results as matrices with given prefix and with index
-   appended. If there is only one matrix, the index is not appended. */
-
-void
-SimResults::writeMat(mat_t *fd, const char *lname) const
-{
-  char tmp[100];
-  for (int i = 0; i < getNumSets(); i++)
-    {
-      if (getNumSets() > 1)
-        sprintf(tmp, "%s_data%d", lname, i+1);
-      else
-        sprintf(tmp, "%s_data", lname);
-      ConstTwoDMatrix m(*(data[i]));
-      m.writeMat(fd, tmp);
-    }
-}
-
-void
-SimResultsStats::simulate(int num_sim, const DecisionRule &dr,
-                          const Vector &start,
-                          const TwoDMatrix &vcov, Journal &journal)
-{
-  SimResults::simulate(num_sim, dr, start, vcov, journal);
-  {
-    JournalRecordPair paa(journal);
-    paa << "Calculating means from the simulations." << endrec;
-    calcMean();
-  }
-  {
-    JournalRecordPair paa(journal);
-    paa << "Calculating covariances from the simulations." << endrec;
-    calcVcov();
-  }
-}
-
-/* Here we do not save the data itself, we save only mean and vcov. */
-void
-SimResultsStats::writeMat(mat_t *fd, const char *lname) const
-{
-  char tmp[100];
-  sprintf(tmp, "%s_mean", lname);
-  ConstTwoDMatrix m(num_y, 1, mean);
-  m.writeMat(fd, tmp);
-  sprintf(tmp, "%s_vcov", lname);
-  ConstTwoDMatrix(vcov).writeMat(fd, tmp);
-}
-
-void
-SimResultsStats::calcMean()
-{
-  mean.zeros();
-  if (data.size()*num_per > 0)
-    {
-      double mult = 1.0/data.size()/num_per;
-      for (auto & i : data)
-        {
-          for (int j = 0; j < num_per; j++)
-            {
-              ConstVector col{i->getCol(j)};
-              mean.add(mult, col);
-            }
-        }
-    }
-}
-
-void
-SimResultsStats::calcVcov()
-{
-  if (data.size()*num_per > 1)
-    {
-      vcov.zeros();
-      double mult = 1.0/(data.size()*num_per - 1);
-      for (auto & i : data)
-        {
-          const TwoDMatrix &d = *i;
-          for (int j = 0; j < num_per; j++)
-            {
-              for (int m = 0; m < num_y; m++)
-                {
-                  for (int n = m; n < num_y; n++)
-                    {
-                      double s = (d.get(m, j)-mean[m])*(d.get(n, j)-mean[n]);
-                      vcov.get(m, n) += mult*s;
-                      if (m != n)
-                        vcov.get(n, m) += mult*s;
-                    }
-                }
-            }
-        }
-    }
-  else
-    {
-      vcov.infs();
-    }
-}
-
-void
-SimResultsDynamicStats::simulate(int num_sim, const DecisionRule &dr,
-                                 const Vector &start,
-                                 const TwoDMatrix &vcov, Journal &journal)
-{
-  SimResults::simulate(num_sim, dr, start, vcov, journal);
-  {
-    JournalRecordPair paa(journal);
-    paa << "Calculating means of the conditional simulations." << endrec;
-    calcMean();
-  }
-  {
-    JournalRecordPair paa(journal);
-    paa << "Calculating variances of the conditional simulations." << endrec;
-    calcVariance();
-  }
-}
-
-void
-SimResultsDynamicStats::writeMat(mat_t *fd, const char *lname) const
-{
-  char tmp[100];
-  sprintf(tmp, "%s_cond_mean", lname);
-  ConstTwoDMatrix(mean).writeMat(fd, tmp);
-  sprintf(tmp, "%s_cond_variance", lname);
-  ConstTwoDMatrix(variance).writeMat(fd, tmp);
-}
-
-void
-SimResultsDynamicStats::calcMean()
-{
-  mean.zeros();
-  if (data.size() > 0)
-    {
-      double mult = 1.0/data.size();
-      for (int j = 0; j < num_per; j++)
-        {
-          Vector meanj{mean.getCol(j)};
-          for (auto & i : data)
-            {
-              ConstVector col{i->getCol(j)};
-              meanj.add(mult, col);
-            }
-        }
-    }
-}
-
-void
-SimResultsDynamicStats::calcVariance()
-{
-  if (data.size() > 1)
-    {
-      variance.zeros();
-      double mult = 1.0/(data.size()-1);
-      for (int j = 0; j < num_per; j++)
-        {
-          ConstVector meanj{mean.getCol(j)};
-          Vector varj{variance.getCol(j)};
-          for (auto & i : data)
-            {
-              Vector col{i->getCol(j)};
-              col.add(-1.0, meanj);
-              for (int k = 0; k < col.length(); k++)
-                col[k] = col[k]*col[k];
-              varj.add(mult, col);
-            }
-        }
-    }
-  else
-    {
-      variance.infs();
-    }
-}
-
-void
-SimResultsIRF::simulate(const DecisionRule &dr, Journal &journal)
-{
-  JournalRecordPair paa(journal);
-  paa << "Performing " << control.getNumSets() << " IRF simulations for "
-      << num_per << " periods; shock=" << ishock << ", impulse=" << imp << endrec;
-  simulate(dr);
-  int thrown = control.getNumSets() - data.size();
-  if (thrown > 0)
-    {
-      JournalRecord rec(journal);
-      rec << "I had to throw " << thrown
-          << " simulations away due to Nan or Inf" << endrec;
-    }
-  calcMeans();
-  calcVariances();
-}
-
-void
-SimResultsIRF::simulate(const DecisionRule &dr)
-{
-  sthread::detach_thread_group gr;
-  for (int idata = 0; idata < control.getNumSets(); idata++)
-    gr.insert(std::make_unique<SimulationIRFWorker>(*this, dr, DecisionRule::horner,
-                                                    num_per, idata, ishock, imp));
-  gr.run();
-}
-
-void
-SimResultsIRF::calcMeans()
-{
-  means.zeros();
-  if (data.size() > 0)
-    {
-      for (auto & i : data)
-        means.add(1.0, *i);
-      means.mult(1.0/data.size());
-    }
-}
-
-void
-SimResultsIRF::calcVariances()
-{
-  if (data.size() > 1)
-    {
-      variances.zeros();
-      for (auto & i : data)
-        {
-          TwoDMatrix d((const TwoDMatrix &)(*i));
-          d.add(-1.0, means);
-          for (int j = 0; j < d.nrows(); j++)
-            for (int k = 0; k < d.ncols(); k++)
-              variances.get(j, k) += d.get(j, k)*d.get(j, k);
-          d.mult(1.0/(data.size()-1));
-        }
-    }
-  else
-    {
-      variances.infs();
-    }
-}
-
-void
-SimResultsIRF::writeMat(mat_t *fd, const char *lname) const
-{
-  char tmp[100];
-  sprintf(tmp, "%s_mean", lname);
-  means.writeMat(fd, tmp);
-  sprintf(tmp, "%s_var", lname);
-  variances.writeMat(fd, tmp);
-}
-
-void
-RTSimResultsStats::simulate(int num_sim, const DecisionRule &dr, const Vector &start,
-                            const TwoDMatrix &v, Journal &journal)
-{
-  JournalRecordPair paa(journal);
-  paa << "Performing " << num_sim << " real-time stochastic simulations for "
-      << num_per << " periods" << endrec;
-  simulate(num_sim, dr, start, v);
-  mean = nc.getMean();
-  mean.add(1.0, dr.getSteady());
-  nc.getVariance(vcov);
-  if (thrown_periods > 0)
-    {
-      JournalRecord rec(journal);
-      rec << "I had to throw " << thrown_periods << " periods away due to Nan or Inf" << endrec;
-      JournalRecord rec1(journal);
-      rec1 << "This affected " << incomplete_simulations << " out of "
-           << num_sim << " simulations" << endrec;
-    }
-}
-
-void
-RTSimResultsStats::simulate(int num_sim, const DecisionRule &dr, const Vector &start,
-                            const TwoDMatrix &vcov)
-{
-  std::vector<RandomShockRealization> rsrs;
-  rsrs.reserve(num_sim);
-
-  sthread::detach_thread_group gr;
-  for (int i = 0; i < num_sim; i++)
-    {
-      RandomShockRealization sr(vcov, system_random_generator.int_uniform());
-      rsrs.push_back(sr);
-      gr.insert(std::make_unique<RTSimulationWorker>(*this, dr, DecisionRule::horner,
-                                                     num_per, start, rsrs.back()));
-    }
-  gr.run();
-}
-
-void
-RTSimResultsStats::writeMat(mat_t *fd, const char *lname)
-{
-  char tmp[100];
-  sprintf(tmp, "%s_rt_mean", lname);
-  ConstTwoDMatrix m(nc.getDim(), 1, mean);
-  m.writeMat(fd, tmp);
-  sprintf(tmp, "%s_rt_vcov", lname);
-  ConstTwoDMatrix(vcov).writeMat(fd, tmp);
-}
-
-IRFResults::IRFResults(const DynamicModel &mod, const DecisionRule &dr,
-                       const SimResults &control, std::vector<int> ili,
-                       Journal &journal)
-  : model(mod), irf_list_ind(std::move(ili))
-{
-  int num_per = control.getNumPer();
-  JournalRecordPair pa(journal);
-  pa << "Calculating IRFs against control for " << (int) irf_list_ind.size() << " shocks and for "
-     << num_per << " periods" << endrec;
-  const TwoDMatrix &vcov = mod.getVcov();
-  for (int ishock : irf_list_ind)
-    {
-      double stderror = sqrt(vcov.get(ishock, ishock));
-      irf_res.push_back(new SimResultsIRF(control, model.numeq(), num_per,
-                                          ishock, stderror));
-      irf_res.push_back(new SimResultsIRF(control, model.numeq(), num_per,
-                                          ishock, -stderror));
-    }
-
-  for (unsigned int ii = 0; ii < irf_list_ind.size(); ii++)
-    {
-      irf_res[2*ii]->simulate(dr, journal);
-      irf_res[2*ii+1]->simulate(dr, journal);
-    }
-}
-
-IRFResults::~IRFResults()
-{
-  for (auto & irf_re : irf_res)
-    delete irf_re;
-}
-
-void
-IRFResults::writeMat(mat_t *fd, const char *prefix) const
-{
-  for (unsigned int i = 0; i < irf_list_ind.size(); i++)
-    {
-      char tmp[100];
-      int ishock = irf_list_ind[i];
-      const char *shockname = model.getExogNames().getName(ishock);
-      sprintf(tmp, "%s_irfp_%s", prefix, shockname);
-      irf_res[2*i]->writeMat(fd, tmp);
-      sprintf(tmp, "%s_irfm_%s", prefix, shockname);
-      irf_res[2*i+1]->writeMat(fd, tmp);
-    }
-}
-
-void
-SimulationWorker::operator()(std::mutex &mut)
-{
-  auto *esr = new ExplicitShockRealization(sr, np);
-  TwoDMatrix *m = dr.simulate(em, np, st, *esr);
-  {
-    std::unique_lock<std::mutex> lk{mut};
-    res.addDataSet(m, esr, st);
-  }
-}
-
-/* Here we create a new instance of |ExplicitShockRealization| of the
-   corresponding control, add the impulse, and simulate. */
-
-void
-SimulationIRFWorker::operator()(std::mutex &mut)
-{
-  auto *esr
-    = new ExplicitShockRealization(res.control.getShocks(idata));
-  esr->addToShock(ishock, 0, imp);
-  TwoDMatrix *m = dr.simulate(em, np, res.control.getStart(idata), *esr);
-  m->add(-1.0, res.control.getData(idata));
-  {
-    std::unique_lock<std::mutex> lk{mut};
-    res.addDataSet(m, esr, res.control.getStart(idata));
-  }
-}
-
-void
-RTSimulationWorker::operator()(std::mutex &mut)
-{
-  NormalConj nc(res.nc.getDim());
-  const PartitionY &ypart = dr.getYPart();
-  int nu = dr.nexog();
-  const Vector &ysteady = dr.getSteady();
-
-  // initialize vectors and subvectors for simulation
-  Vector dyu(ypart.nys()+nu);
-  ConstVector ystart_pred(ystart, ypart.nstat, ypart.nys());
-  ConstVector ysteady_pred(ysteady, ypart.nstat, ypart.nys());
-  Vector dy(dyu, 0, ypart.nys());
-  Vector u(dyu, ypart.nys(), nu);
-  Vector y(nc.getDim());
-  ConstVector ypred(y, ypart.nstat, ypart.nys());
-
-  // simulate the first real-time period
-  int ip = 0;
-  dy = ystart_pred;
-  dy.add(-1.0, ysteady_pred);
-  sr.get(ip, u);
-  dr.eval(em, y, dyu);
-  if (ip >= res.num_burn)
-    nc.update(y);
-
-  // simulate other real-time periods@>=
-  while (y.isFinite() && ip < res.num_burn + res.num_per)
-    {
-      ip++;
-      dy = ypred;
-      sr.get(ip, u);
-      dr.eval(em, y, dyu);
-      if (ip >= res.num_burn)
-        nc.update(y);
-    }
-  {
-    std::unique_lock<std::mutex> lk{mut};
-    res.nc.update(nc);
-    if (res.num_per-ip > 0)
-      {
-        res.incomplete_simulations++;
-        res.thrown_periods += res.num_per-ip;
-      }
-  }
-}
-
-/* This calculates factorization $FF^T=V$ in the Cholesky way. It does
-   not work for semidefinite matrices. */
-
-void
-RandomShockRealization::choleskyFactor(const ConstTwoDMatrix &v)
-{
-  factor = v;
-  lapack_int rows = factor.nrows(), lda = factor.getLD();
-  for (int i = 0; i < rows; i++)
-    for (int j = i+1; j < rows; j++)
-      factor.get(i, j) = 0.0;
-  lapack_int info;
-
-  dpotrf("L", &rows, factor.base(), &lda, &info);
-  KORD_RAISE_IF(info != 0,
-                "Info!=0 in RandomShockRealization::choleskyFactor");
-}
-
-/* This calculates $FF^T=V$ factorization by symmetric Schur
-   decomposition. It works for semidifinite matrices. */
-
-void
-RandomShockRealization::schurFactor(const ConstTwoDMatrix &v)
-{
-  SymSchurDecomp ssd(v);
-  ssd.getFactor(factor);
-}
-
-void
-RandomShockRealization::get(int n, Vector &out)
-{
-  KORD_RAISE_IF(out.length() != numShocks(),
-                "Wrong length of out vector in RandomShockRealization::get");
-  Vector d(out.length());
-  for (int i = 0; i < d.length(); i++)
-    {
-      d[i] = mtwister.normal();
-    }
-  out.zeros();
-  factor.multaVec(out, ConstVector(d));
-}
-
-ExplicitShockRealization::ExplicitShockRealization(ShockRealization &sr,
-                                                   int num_per)
-  : shocks(sr.numShocks(), num_per)
-{
-  for (int j = 0; j < num_per; j++)
-    {
-      Vector jcol{shocks.getCol(j)};
-      sr.get(j, jcol);
-    }
-}
-
-void
-ExplicitShockRealization::get(int n, Vector &out)
-{
-  KORD_RAISE_IF(out.length() != numShocks(),
-                "Wrong length of out vector in ExplicitShockRealization::get");
-  int i = n % shocks.ncols();
-  ConstVector icol{shocks.getCol(i)};
-  out = icol;
-}
-
-void
-ExplicitShockRealization::addToShock(int ishock, int iper, double val)
-{
-  KORD_RAISE_IF(ishock < 0 || ishock > numShocks(),
-                "Wrong index of shock in ExplicitShockRealization::addToShock");
-  int j = iper % shocks.ncols();
-  shocks.get(ishock, j) += val;
-}
-
-void
-GenShockRealization::get(int n, Vector &out)
-{
-  KORD_RAISE_IF(out.length() != numShocks(),
-                "Wrong length of out vector in GenShockRealization::get");
-  ExplicitShockRealization::get(n, out);
-  Vector r(numShocks());
-  RandomShockRealization::get(n, r);
-  for (int j = 0; j < numShocks(); j++)
-    if (!std::isfinite(out[j]))
-      out[j] = r[j];
-}

dynare++/kord/decision_rule.hh

-// Copyright 2004, Ondra Kamenik
-
-// Decision rule and simulation
-
-/* The main purpose of this file is a decision rule representation which
-   can run a simulation. So we define an interface for classes providing
-   realizations of random shocks, and define the class
-   |DecisionRule|. The latter basically takes tensor container of
-   derivatives of policy rules, and adds them up with respect to
-   $\sigma$. The class allows to specify the $\sigma$ different from $1$.
-
-   In addition, we provide classes for running simulations and storing
-   the results, calculating some statistics and generating IRF. The class
-   |DRFixPoint| allows for calculation of the fix point of a given
-   decision rule. */
-
-#ifndef DECISION_RULE_H
-#define DECISION_RULE_H
-
-#include <matio.h>
-
-#include "kord_exception.hh"
-#include "korder.hh"
-#include "normal_conjugate.hh"
-#include "mersenne_twister.hh"
-
-/* This is a general interface to a shock realizations. The interface
-   has only one method returning the shock realizations at the given
-   time. This method is not constant, since it may change a state of the
-   object. */
-
-class ShockRealization
-{
-public:
-  virtual ~ShockRealization()
-  = default;
-  virtual void get(int n, Vector &out) = 0;
-  virtual int numShocks() const = 0;
-};
-
-/* This class is an abstract interface to decision rule. Its main
-   purpose is to define a common interface for simulation of a decision
-   rule. We need only a simulate, evaluate, cetralized clone and output
-   method. The |simulate| method simulates the rule for a given
-   realization of the shocks. |eval| is a primitive evaluation (it takes
-   a vector of state variables (predetermined, both and shocks) and
-   returns the next period variables. Both input and output are in
-   deviations from the rule's steady. |evaluate| method makes only one
-   step of simulation (in terms of absolute values, not
-   deviations). |centralizedClone| returns a new copy of the decision
-   rule, which is centralized about provided fix-point. And finally
-   |writeMat| writes the decision rule to the MAT file. */
-
-class DecisionRule
-{
-public:
-  enum emethod { horner, trad };
-  virtual ~DecisionRule()
-  = default;
-  virtual TwoDMatrix *simulate(emethod em, int np, const ConstVector &ystart,
-                               ShockRealization &sr) const = 0;
-  virtual void eval(emethod em, Vector &out, const ConstVector &v) const = 0;
-  virtual void evaluate(emethod em, Vector &out, const ConstVector &ys,
-                        const ConstVector &u) const = 0;
-  virtual void writeMat(mat_t *fd, const char *prefix) const = 0;
-  virtual DecisionRule *centralizedClone(const Vector &fixpoint) const = 0;
-  virtual const Vector&getSteady() const = 0;
-  virtual int nexog() const = 0;
-  virtual const PartitionY&getYPart() const = 0;
-};
-
-/* The main purpose of this class is to implement |DecisionRule|
-   interface, which is a simulation. To be able to do this we have to
-   know the partitioning of state vector $y$ since we will need to pick
-   only predetermined part $y^*$. Also, we need to know the steady state.
-
-   The decision rule will take the form: $$y_t-\bar
-   y=\sum_{i=0}^n\left[g_{(yu)^i}\right]_{\alpha_1\ldots\alpha_i}\prod_{m=1}^i
-   \left[\matrix{y^*_{t-1}-\bar y^*\cr u_t}\right]^{\alpha_m},$$ where
-   the tensors $\left[g_{(yu)^i}\right]$ are tensors of the constructed
-   container, and $\bar y$ is the steady state.
-
-   If we know the fix point of the rule (conditional zero shocks)
-   $\tilde y$, the rule can be transformed to so called ``centralized''
-   form. This is very similar to the form above but the zero dimensional
-   tensor is zero:
-   $$y_t-\tilde y=\sum_{i=1}^n
-   \left[\tilde g_{(yu)^i}\right]_{\alpha_1\ldots\alpha_i}\prod_{m=1}^i
-   \left[\matrix{y^*_{t-1}-\tilde y^*\cr u_t}\right]^{\alpha_m}.$$
-   We provide a method and a constructor to transform a rule to the centralized form.
-
-   The class is templated, the template argument is either |KOrder::fold|
-   or |KOrder::unfold|. So, there are two implementations of |DecisionRule|
-   interface. */
-
-template <int t>
-class DecisionRuleImpl : public ctraits<t>::Tpol, public DecisionRule
-{
-protected:
-  using _Tparent = typename ctraits<t>::Tpol;
-  const Vector ysteady;
-  const PartitionY ypart;
-  const int nu;
-public:
-  DecisionRuleImpl(const _Tparent &pol, const PartitionY &yp, int nuu,
-                   const ConstVector &ys)
-    : ctraits<t>::Tpol(pol), ysteady(ys), ypart(yp), nu(nuu)
-  {
-  }
-  DecisionRuleImpl(_Tparent &pol, const PartitionY &yp, int nuu,
-                   const ConstVector &ys)
-    : ctraits<t>::Tpol(0, yp.ny(), pol), ysteady(ys), ypart(yp),
-    nu(nuu)
-  {
-  }
-  DecisionRuleImpl(const _Tg &g, const PartitionY &yp, int nuu,
-                   const ConstVector &ys, double sigma)
-    : ctraits<t>::Tpol(yp.ny(), yp.nys()+nuu), ysteady(ys), ypart(yp), nu(nuu)
-  {
-    fillTensors(g, sigma);
-  }
-  DecisionRuleImpl(const DecisionRuleImpl<t> &dr, const ConstVector &fixpoint)
-    : ctraits<t>::Tpol(dr.ypart.ny(), dr.ypart.nys()+dr.nu),
-    ysteady(fixpoint), ypart(dr.ypart), nu(dr.nu)
-  {
-    centralize(dr);
-  }
-  const Vector &
-  getSteady() const override
-  {
-    return ysteady;
-  }
-  TwoDMatrix *simulate(emethod em, int np, const ConstVector &ystart,
-                       ShockRealization &sr) const override;
-  void evaluate(emethod em, Vector &out, const ConstVector &ys,
-                const ConstVector &u) const override;
-  DecisionRule *centralizedClone(const Vector &fixpoint) const override;
-  void writeMat(mat_t *fd, const char *prefix) const override;
-
-  int
-  nexog() const override
-  {
-    return nu;
-  }
-  const PartitionY &
-  getYPart() const override
-  {
-    return ypart;
-  }
-protected:
-  void fillTensors(const _Tg &g, double sigma);
-  void centralize(const DecisionRuleImpl &dr);
-  void eval(emethod em, Vector &out, const ConstVector &v) const override;
-};
-
-/* Here we have to fill the tensor polynomial. This involves two
-   separated actions. First is to evaluate the approximation at a given
-   $\sigma$, the second is to compile the tensors $[g_{{(yu)}^{i+j}}]$ from
-   $[g_{y^iu^j}]$. The first action is done here, the second is done by
-   method |addSubTensor| of a full symmetry tensor.
-
-   The way how the evaluation is done is described here:
-
-   The $q-$order approximation to the solution can be written as:
-
-   $$
-   \eqalign{
-   y_t-\bar y &= \sum_{l=1}^q{1\over l!}\left[\sum_{i+j+k=l}
-   \left(\matrix{l\cr i,j,k}\right)\left[g_{y^iu^j\sigma^k}\right]
-   _{\alpha_1\ldots\alpha_j\beta_1\ldots\beta_j}
-   \prod_{m=1}^i[y^*_{t-1}-\bar y^*]^{\alpha_m}
-   \prod_{n=1}^j[u_t]^{\beta_m}\sigma^k\right]\cr
-   &= \sum_{l=1}^q\left[\sum_{i+j\leq l}\left(\matrix{i+j\cr i}\right)
-   \left[\sum_{k=0}^{l-i-j}{1\over l!}
-   \left(\matrix{l\cr k}\right)\left[g_{y^iu^j\sigma^k}\right]\sigma^k\right]
-   \prod_{m=1}^i[y^*_{t-1}-\bar y^*]^{\alpha_m}
-   \prod_{n=1}^j[u_t]^{\beta_m}\sigma^k\right]
-   }
-   $$
-
-   This means that for each $i+j+k=l$ we have to add
-   $${1\over l!}\left(\matrix{l\cr
-   k}\right)\left[g_{y^iu^j\sigma^k}\right]\cdot\sigma^k=
-   {1\over (i+j)!k!}\left[g_{y^iu^j\sigma^k}\right]\cdot\sigma^k$$ to
-   $g_{(yu)^{i+j}}$. In addition, note that the multiplier
-   $\left(\matrix{i+j\cr i}\right)$ is applied when the fully symmetric
-   tensor $[g_{(yu)^{i+j}}]$ is evaluated.
-
-   So we go through $i+j=d=0\ldots q$ and in each loop we form the fully
-   symmetric tensor $[g_{(yu)^l}]$ and insert it to the container. */
-
-template <int t>
-void
-DecisionRuleImpl<t>::fillTensors(const _Tg &g, double sigma)
-{
-  IntSequence tns(2);
-  tns[0] = ypart.nys(); tns[1] = nu;
-  int dfact = 1;
-  for (int d = 0; d <= g.getMaxDim(); d++, dfact *= d)
-    {
-      auto *g_yud = new _Ttensym(ypart.ny(), ypart.nys()+nu, d);
-      g_yud->zeros();
-
-      // fill tensor of |g_yud| of dimension |d|
-      /* Here we have to fill the tensor $\left[g_{(yu)^d}\right]$. So we go
-         through all pairs $(i,j)$ giving $i+j=d$, and through all $k$ from
-         zero up to maximal dimension minus $d$. In this way we go through all
-         symmetries of $g_{y^iu^j\sigma^k}$ which will be added to $g_{(yu)^d}$.
-
-         Note that at the beginning, |dfact| is a factorial of |d|. We
-         calculate |kfact| is equal to $k!$. As indicated in
-         |@<|DecisionRuleImpl::fillTensors| code@>|, the added tensor is thus
-         multiplied with ${1\over d!k!}\sigma^k$. */
-
-      for (int i = 0; i <= d; i++)
-        {
-          int j = d-i;
-          int kfact = 1;
-          _Ttensor tmp(ypart.ny(),
-                       TensorDimens(Symmetry(i, j), tns));
-          tmp.zeros();
-          for (int k = 0; k+d <= g.getMaxDim(); k++, kfact *= k)
-            {
-              Symmetry sym(i, j, 0, k);
-              if (g.check(sym))
-                {
-                  double mult = pow(sigma, k)/dfact/kfact;
-                  tmp.add(mult, *(g.get(sym)));
-                }
-            }
-          g_yud->addSubTensor(tmp);
-        }
-
-      this->insert(g_yud);
-    }
-}
-
-/* The centralization is straightforward. We suppose here that the
-   object's steady state is the fix point $\tilde y$. It is clear that
-   the new derivatives $\left[\tilde g_{(yu)^i}\right]$ will be equal to
-   the derivatives of the original decision rule |dr| at the new steady
-   state $\tilde y$. So, the new derivatives are obtained by derivating the
-   given decision rule $dr$ and evaluating its polynomial at
-   $$dstate=\left[\matrix{\tilde y^*-\bar y^*\cr 0}\right],$$
-   where $\bar y$ is the steady state of the original rule |dr|. */
-
-template <int t>
-void
-DecisionRuleImpl<t>::centralize(const DecisionRuleImpl &dr)
-{
-  Vector dstate(ypart.nys() + nu);
-  dstate.zeros();
-  Vector dstate_star(dstate, 0, ypart.nys());
-  ConstVector newsteady_star(ysteady, ypart.nstat, ypart.nys());
-  ConstVector oldsteady_star(dr.ysteady, ypart.nstat, ypart.nys());
-  dstate_star.add(1.0, newsteady_star);
-  dstate_star.add(-1.0, oldsteady_star);
-
-  _Tpol pol(dr);
-  int dfac = 1;
-  for (int d = 1; d <= dr.getMaxDim(); d++, dfac *= d)
-    {
-      pol.derivative(d-1);
-      _Ttensym *der = pol.evalPartially(d, dstate);
-      der->mult(1.0/dfac);
-      this->insert(der);
-    }
-}
-
-/* Here we evaluate repeatedly the polynomial storing results in the
-   created matrix. For exogenous shocks, we use |ShockRealization|
-   class, for predetermined variables, we use |ystart| as the first
-   state. The |ystart| vector is required to be all state variables
-   |ypart.ny()|, although only the predetermined part of |ystart| is
-   used.
-
-   We simulate in terms of $\Delta y$, this is, at the beginning the
-   |ysteady| is canceled from |ystart|, we simulate, and at the end
-   |ysteady| is added to all columns of the result. */
-
-template <int t>
-TwoDMatrix *
-DecisionRuleImpl<t>::simulate(emethod em, int np, const ConstVector &ystart,
-                              ShockRealization &sr) const
-{
-  KORD_RAISE_IF(ysteady.length() != ystart.length(),
-                "Start and steady lengths differ in DecisionRuleImpl::simulate");
-  auto *res = new TwoDMatrix(ypart.ny(), np);
-
-  // initialize vectors and subvectors for simulation
-  /* Here allocate the stack vector $(\Delta y^*, u)$, define the
-     subvectors |dy|, and |u|, then we pickup predetermined parts of
-     |ystart| and |ysteady|. */
-  Vector dyu(ypart.nys()+nu);
-  ConstVector ystart_pred(ystart, ypart.nstat, ypart.nys());
-  ConstVector ysteady_pred(ysteady, ypart.nstat, ypart.nys());
-  Vector dy(dyu, 0, ypart.nys());
-  Vector u(dyu, ypart.nys(), nu);
-
-  // perform the first step of simulation
-  /* We cancel |ysteady| from |ystart|, get realization to |u|, and
-     evaluate the polynomial. */
-  dy = ystart_pred;
-  dy.add(-1.0, ysteady_pred);
-  sr.get(0, u);
-  Vector out{res->getCol(0)};
-  eval(em, out, dyu);
-
-  // perform all other steps of simulations
-  /* Also clear. If the result at some period is not finite, we pad the
-     rest of the matrix with zeros. */
-  int i = 1;
-  while (i < np)
-    {
-      ConstVector ym{res->getCol(i-1)};
-      ConstVector dym(ym, ypart.nstat, ypart.nys());
-      dy = dym;
-      sr.get(i, u);
-      Vector out{res->getCol(i)};
-      eval(em, out, dyu);
-      if (!out.isFinite())
-        {
-          if (i+1 < np)
-            {
-              TwoDMatrix rest(*res, i+1, np-i-1);
-              rest.zeros();
-            }
-          break;
-        }
-      i++;
-    }
-
-  // add the steady state to columns of |res|
-  /* Even clearer. We add the steady state to the numbers computed above
-     and leave the padded columns to zero. */
-  for (int j = 0; j < i; j++)
-    {
-      Vector col{res->getCol(j)};
-      col.add(1.0, ysteady);
-    }
-
-  return res;
-}
-
-/* This is one period evaluation of the decision rule. The simulation
-   is a sequence of repeated one period evaluations with a difference,
-   that the steady state (fix point) is cancelled and added once. Hence
-   we have two special methods. */
-
-template <int t>
-void
-DecisionRuleImpl<t>::evaluate(emethod em, Vector &out, const ConstVector &ys,
-                              const ConstVector &u) const
-{
-  KORD_RAISE_IF(ys.length() != ypart.nys() || u.length() != nu,
-                "Wrong dimensions of input vectors in DecisionRuleImpl::evaluate");
-  KORD_RAISE_IF(out.length() != ypart.ny(),
-                "Wrong dimension of output vector in DecisionRuleImpl::evaluate");
-  ConstVector ysteady_pred(ysteady, ypart.nstat, ypart.nys());
-  Vector ys_u(ypart.nys()+nu);
-  Vector ys_u1(ys_u, 0, ypart.nys());
-  ys_u1 = ys;
-  ys_u1.add(-1.0, ysteady_pred);
-  Vector ys_u2(ys_u, ypart.nys(), nu);
-  ys_u2 = u;
-  eval(em, out, ys_u);
-  out.add(1.0, ysteady);
-}
-
-/* This is easy. We just return the newly created copy using the
-   centralized constructor. */
-
-template <int t>
-DecisionRule *
-DecisionRuleImpl<t>::centralizedClone(const Vector &fixpoint) const
-{
-  return new DecisionRuleImpl<t>(*this, fixpoint);
-}
-
-/* Here we only encapsulate two implementations to one, deciding
-   according to the parameter. */
-
-template <int t>
-void
-DecisionRuleImpl<t>::eval(emethod em, Vector &out, const ConstVector &v) const
-{
-  if (em == DecisionRule::horner)
-    _Tparent::evalHorner(out, v);
-  else
-    _Tparent::evalTrad(out, v);
-}
-
-/* Write the decision rule and steady state to the MAT file. */
-
-template <int t>
-void
-DecisionRuleImpl<t>::writeMat(mat_t *fd, const char *prefix) const
-{
-  ctraits<t>::Tpol::writeMat(fd, prefix);
-  TwoDMatrix dum(ysteady.length(), 1);
-  dum.getData() = ysteady;
-  char tmp[100];
-  sprintf(tmp, "%s_ss", prefix);
-  ConstTwoDMatrix(dum).writeMat(fd, tmp);
-}
-
-/* This is exactly the same as |DecisionRuleImpl<KOrder::fold>|. The
-   only difference is that we have a conversion from
-   |UnfoldDecisionRule|, which is exactly
-   |DecisionRuleImpl<KOrder::unfold>|. */
-
-class UnfoldDecisionRule;
-class FoldDecisionRule : public DecisionRuleImpl<KOrder::fold>
-{
-  friend class UnfoldDecisionRule;
-public:
-  FoldDecisionRule(const ctraits<KOrder::fold>::Tpol &pol, const PartitionY &yp, int nuu,
-                   const ConstVector &ys)
-    : DecisionRuleImpl<KOrder::fold>(pol, yp, nuu, ys)
-  {
-  }
-  FoldDecisionRule(ctraits<KOrder::fold>::Tpol &pol, const PartitionY &yp, int nuu,
-                   const ConstVector &ys)
-    : DecisionRuleImpl<KOrder::fold>(pol, yp, nuu, ys)
-  {
-  }
-  FoldDecisionRule(const ctraits<KOrder::fold>::Tg &g, const PartitionY &yp, int nuu,
-                   const ConstVector &ys, double sigma)
-    : DecisionRuleImpl<KOrder::fold>(g, yp, nuu, ys, sigma)
-  {
-  }
-  FoldDecisionRule(const DecisionRuleImpl<KOrder::fold> &dr, const ConstVector &fixpoint)
-    : DecisionRuleImpl<KOrder::fold>(dr, fixpoint)
-  {
-  }
-  FoldDecisionRule(const UnfoldDecisionRule &udr);
-};
-
-/* This is exactly the same as |DecisionRuleImpl<KOrder::unfold>|, but
-   with a conversion from |FoldDecisionRule|, which is exactly
-   |DecisionRuleImpl<KOrder::fold>|. */
-
-class UnfoldDecisionRule : public DecisionRuleImpl<KOrder::unfold>
-{
-  friend class FoldDecisionRule;
-public:
-  UnfoldDecisionRule(const ctraits<KOrder::unfold>::Tpol &pol, const PartitionY &yp, int nuu,
-                     const ConstVector &ys)
-    : DecisionRuleImpl<KOrder::unfold>(pol, yp, nuu, ys)
-  {
-  }
-  UnfoldDecisionRule(ctraits<KOrder::unfold>::Tpol &pol, const PartitionY &yp, int nuu,
-                     const ConstVector &ys)
-    : DecisionRuleImpl<KOrder::unfold>(pol, yp, nuu, ys)
-  {
-  }
-  UnfoldDecisionRule(const ctraits<KOrder::unfold>::Tg &g, const PartitionY &yp, int nuu,
-                     const ConstVector &ys, double sigma)
-    : DecisionRuleImpl<KOrder::unfold>(g, yp, nuu, ys, sigma)
-  {
-  }
-  UnfoldDecisionRule(const DecisionRuleImpl<KOrder::unfold> &dr, const ConstVector &fixpoint)
-    : DecisionRuleImpl<KOrder::unfold>(dr, fixpoint)
-  {
-  }
-  UnfoldDecisionRule(const FoldDecisionRule &udr);
-};
-
-/* This class serves for calculation of the fix point of the decision
-   rule given that the shocks are zero. The class is very similar to the
-   |DecisionRuleImpl|. Besides the calculation of the fix point, the only
-   difference between |DRFixPoint| and |DecisionRuleImpl| is that the
-   derivatives wrt. shocks are ignored (since shocks are zero during the
-   calculations). That is why have a different |fillTensor| method.
-
-   The solution algorithm is Newton and is described in
-   |@<|DRFixPoint::solveNewton| code@>|. It solves $F(y)=0$, where
-   $F=g(y,0)-y$. The function $F$ is given by its derivatives |bigf|. The
-   Jacobian of the solved system is given by derivatives stored in
-   |bigfder|. */
-
-template <int t>
-class DRFixPoint : public ctraits<t>::Tpol
-{
-  using _Tparent = typename ctraits<t>::Tpol;
-  static int max_iter;
-  static int max_newton_iter;
-  static int newton_pause;
-  static double tol;
-  const Vector ysteady;
-  const PartitionY ypart;
-  _Tparent *bigf;
-  _Tparent *bigfder;
-public:
-  using emethod = typename DecisionRule::emethod;
-  DRFixPoint(const _Tg &g, const PartitionY &yp,
-             const Vector &ys, double sigma);
-  
-  ~DRFixPoint() override;
-
-  bool calcFixPoint(emethod em, Vector &out);
-
-  int
-  getNumIter() const
-  {
-    return iter;
-  }
-  int
-  getNewtonLastIter() const
-  {
-    return newton_iter_last;
-  }
-  int
-  getNewtonTotalIter() const
-  {
-    return newton_iter_total;
-  }
-protected:
-  void fillTensors(const _Tg &g, double sigma);
-  bool solveNewton(Vector &y);
-private:
-  int iter;
-  int newton_iter_last;
-  int newton_iter_total;
-};
-
-/* Here we have to setup the function $F=g(y,0)-y$ and ${\partial
-   F\over\partial y}$. The former is taken from the given derivatives of
-   $g$ where a unit matrix is subtracted from the first derivative
-   (|Symmetry(1)|). Then the derivative of the $F$ polynomial is
-   calculated. */
-
-template <int t>
-DRFixPoint<t>::DRFixPoint(const _Tg &g, const PartitionY &yp,
-                          const Vector &ys, double sigma)
-  : ctraits<t>::Tpol(yp.ny(), yp.nys()),
-  ysteady(ys), ypart(yp), bigf(nullptr), bigfder(nullptr)
-{
-  fillTensors(g, sigma);
-  _Tparent yspol(ypart.nstat, ypart.nys(), *this);
-  bigf = new _Tparent((const _Tparent &) yspol);
-  _Ttensym *frst = bigf->get(Symmetry(1));
-  for (int i = 0; i < ypart.nys(); i++)
-    frst->get(i, i) = frst->get(i, i) - 1;
-  bigfder = new _Tparent(*bigf, 0);
-}
-
-template <int t>
-DRFixPoint<t>::~DRFixPoint()
-{
-  if (bigf)
-    delete bigf;
-  if (bigfder)
-    delete bigfder;
-}
-
-/* Here we fill the tensors for the |DRFixPoint| class. We ignore the
-   derivatives $g_{y^iu^j\sigma^k}$ for which $j>0$. So we go through all
-   dimensions |d|, and all |k| such that |d+k| is between the maximum
-   dimension and |d|, and add ${\sigma^k\over d!k!}g_{y^d\sigma^k}$ to
-   the tensor $g_{(y)^d}$. */
-
-template <int t>
-void
-DRFixPoint<t>::fillTensors(const _Tg &g, double sigma)
-{
-  int dfact = 1;
-  for (int d = 0; d <= g.getMaxDim(); d++, dfact *= d)
-    {
-      auto *g_yd = new _Ttensym(ypart.ny(), ypart.nys(), d);
-      g_yd->zeros();
-      int kfact = 1;
-      for (int k = 0; d+k <= g.getMaxDim(); k++, kfact *= k)
-        {
-          if (g.check(Symmetry(d, 0, 0, k)))
-            {
-              const _Ttensor *ten = g.get(Symmetry(d, 0, 0, k));
-              double mult = pow(sigma, k)/dfact/kfact;
-              g_yd->add(mult, *ten);
-            }
-        }
-      this->insert(g_yd);
-    }
-}
-
-/* This tries to solve polynomial equation $F(y)=0$, where $F$
-   polynomial is |bigf| and its derivative is in |bigfder|. It returns
-   true if the Newton converged. The method takes the given vector as
-   initial guess, and rewrites it with a solution. The method guarantees
-   to return the vector, which has smaller norm of the residual. That is
-   why the input/output vector |y| is always changed.
-
-   The method proceeds with a Newton step, if the Newton step improves
-   the residual error. So we track residual errors in |flastnorm| and
-   |fnorm| (former and current). In addition, at each step we search for
-   an underrelaxation parameter |urelax|, which improves the residual. If
-   |urelax| is less that |urelax_threshold|, we stop searching and stop
-   the Newton. */
-
-template <int t>
-bool
-DRFixPoint<t>::solveNewton(Vector &y)
-{
-  const double urelax_threshold = 1.e-5;
-  Vector sol((const Vector &)y);
-  Vector delta(y.length());
-  newton_iter_last = 0;
-  bool delta_finite = true;
-  double flastnorm = 0.0;
-  double fnorm = 0.0;
-  bool converged = false;
-  double urelax = 1.0;
-
-  do
-    {
-      _Ttensym *jacob = bigfder->evalPartially(1, sol);
-      bigf->evalHorner(delta, sol);
-      if (newton_iter_last == 0)
-        flastnorm = delta.getNorm();
-      delta_finite = delta.isFinite();
-      if (delta_finite)
-        {
-          ConstTwoDMatrix(*jacob).multInvLeft(delta);
-
-          // find |urelax| improving residual
-          /* Here we find the |urelax|. We cycle as long as the new residual size
-             |fnorm| is greater than last residual size |flastnorm|. If the urelax
-             is less than |urelax_threshold| we give up. The |urelax| is damped by
-             the ratio of |flastnorm| and |fnorm|. It the ratio is close to one, we
-             damp by one half. */
-          bool urelax_found = false;
-          urelax = 1.0;
-          while (!urelax_found && urelax > urelax_threshold)
-            {
-              Vector soltmp((const Vector &)sol);
-              soltmp.add(-urelax, delta);
-              Vector f(sol.length());
-              bigf->evalHorner(f, soltmp);
-              fnorm = f.getNorm();
-              if (fnorm <= flastnorm)
-                urelax_found = true;
-              else
-                urelax *= std::min(0.5, flastnorm/fnorm);
-            }
-
-          sol.add(-urelax, delta);
-          delta_finite = delta.isFinite();
-        }
-      delete jacob;
-      newton_iter_last++;
-      converged = delta_finite && fnorm < tol;
-      flastnorm = fnorm;
-    }
-  while (!converged && newton_iter_last < max_newton_iter
-         &&urelax > urelax_threshold);
-
-  newton_iter_total += newton_iter_last;
-  if (!converged)
-    newton_iter_last = 0;
-  y = (const Vector &) sol;
-  return converged;
-}
-
-/* This method solves the fix point of the no-shocks rule
-   $y_{t+1}=f(y_t)$. It combines dull steps with Newton attempts. The
-   dull steps correspond to evaluations setting $y_{t+1}=f(y_t)$. For
-   reasonable models the dull steps converge to the fix-point but very
-   slowly. That is why we make Newton attempt from time to time. The
-   frequency of the Newton attempts is given by |newton_pause|. We
-   perform the calculations in deviations from the steady state. So, at
-   the end, we have to add the steady state.
-
-   The method also sets the members |iter|, |newton_iter_last| and
-   |newton_iter_total|. These numbers can be examined later.
-
-   The |out| vector is not touched if the algorithm has not convered. */
-
-template <int t>
-bool
-DRFixPoint<t>::calcFixPoint(emethod em, Vector &out)
-{
-  KORD_RAISE_IF(out.length() != ypart.ny(),
-                "Wrong length of out in DRFixPoint::calcFixPoint");
-
-  Vector delta(ypart.nys());
-  Vector ystar(ypart.nys());
-  ystar.zeros();
-
-  iter = 0;
-  newton_iter_last = 0;
-  newton_iter_total = 0;
-  bool converged = false;
-  do
-    {
-      if ((iter/newton_pause)*newton_pause == iter)
-        converged = solveNewton(ystar);
-      if (!converged)
-        {
-          bigf->evalHorner(delta, ystar);
-          KORD_RAISE_IF_X(!delta.isFinite(),
-                          "NaN or Inf asserted in DRFixPoint::calcFixPoint",
-                          KORD_FP_NOT_FINITE);
-          ystar.add(1.0, delta);
-          converged = delta.getNorm() < tol;
-        }
-      iter++;
-    }
-  while (iter < max_iter && !converged);
-
-  if (converged)
-    {
-      _Tparent::evalHorner(out, ystar);
-      out.add(1.0, ysteady);
-    }
-
-  return converged;
-}
-
-/* This is a basically a number of matrices of the same dimensions,
-   which can be obtained as simulation results from a given decision rule
-   and shock realizations. We also store the realizations of shocks and the
-   starting point of each simulation. */
-
-class ExplicitShockRealization;
-class SimResults
-{
-protected:
-  int num_y;
-  int num_per;
-  int num_burn;
-  std::vector<TwoDMatrix *> data;
-  std::vector<ExplicitShockRealization *> shocks;
-  std::vector<ConstVector> start;
-public:
-  SimResults(int ny, int nper, int nburn = 0)
-    : num_y(ny), num_per(nper), num_burn(nburn)
-  {
-  }
-  virtual
-  ~SimResults();
-  void simulate(int num_sim, const DecisionRule &dr, const Vector &start,
-                const TwoDMatrix &vcov, Journal &journal);
-  void simulate(int num_sim, const DecisionRule &dr, const Vector &start,
-                const TwoDMatrix &vcov);
-  int
-  getNumPer() const
-  {
-    return num_per;
-  }
-  int
-  getNumBurn() const
-  {
-    return num_burn;
-  }
-  int
-  getNumSets() const
-  {
-    return (int) data.size();
-  }
-  const TwoDMatrix &
-  getData(int i) const
-  {
-    return *(data[i]);
-  }
-  const ExplicitShockRealization &
-  getShocks(int i) const
-  {
-    return *(shocks[i]);
-  }
-  const ConstVector &
-  getStart(int i) const
-  {
-    return start[i];
-  }
-
-  bool addDataSet(TwoDMatrix *d, ExplicitShockRealization *sr, const ConstVector &st);
-  void writeMat(const char *base, const char *lname) const;
-  void writeMat(mat_t *fd, const char *lname) const;
-};
-
-/* This does the same as |SimResults| plus it calculates means and
-   covariances of the simulated data. */
-
-class SimResultsStats : public SimResults
-{
-protected:
-  Vector mean;
-  TwoDMatrix vcov;
-public:
-  SimResultsStats(int ny, int nper, int nburn = 0)
-    : SimResults(ny, nper, nburn), mean(ny), vcov(ny, ny)
-  {
-  }
-  void simulate(int num_sim, const DecisionRule &dr, const Vector &start,
-                const TwoDMatrix &vcov, Journal &journal);
-  void writeMat(mat_t *fd, const char *lname) const;
-protected:
-  void calcMean();
-  void calcVcov();
-};
-
-/* This does the similar thing as |SimResultsStats| but the statistics are
-   not calculated over all periods but only within each period. Then we
-   do not calculate covariances with periods but only variances. */
-
-class SimResultsDynamicStats : public SimResults
-{
-protected:
-  TwoDMatrix mean;
-  TwoDMatrix variance;
-public:
-  SimResultsDynamicStats(int ny, int nper, int nburn = 0)
-    : SimResults(ny, nper, nburn), mean(ny, nper), variance(ny, nper)
-  {
-  }
-  void simulate(int num_sim, const DecisionRule &dr, const Vector &start,
-                const TwoDMatrix &vcov, Journal &journal);
-  void writeMat(mat_t *fd, const char *lname) const;
-protected:
-  void calcMean();
-  void calcVariance();
-};
-
-/* This goes through control simulation results, and for each control
-   it adds a given impulse to a given shock and runs a simulation. The
-   control simulation is then cancelled and the result is stored. After
-   that these results are averaged with variances calculated.
-
-   The means and the variances are then written to the MAT-4 file. */
-
-class SimulationIRFWorker;
-class SimResultsIRF : public SimResults
-{
-  friend class SimulationIRFWorker;
-protected:
-  const SimResults &control;
-  int ishock;
-  double imp;
-  TwoDMatrix means;
-  TwoDMatrix variances;
-public:
-  SimResultsIRF(const SimResults &cntl, int ny, int nper, int i, double impulse)
-    : SimResults(ny, nper, 0), control(cntl),
-      ishock(i), imp(impulse),
-      means(ny, nper), variances(ny, nper)
-  {
-  }
-  void simulate(const DecisionRule &dr, Journal &journal);
-  void simulate(const DecisionRule &dr);
-  void writeMat(mat_t *fd, const char *lname) const;
-protected:
-  void calcMeans();
-  void calcVariances();
-};
-
-/* This simulates and gathers all statistics from the real time
-   simulations. In the |simulate| method, it runs |RTSimulationWorker|s
-   which accummulate information from their own estimates. The estimation
-   is done by means of |NormalConj| class, which is a conjugate family of
-   densities for normal distibutions. */
-
-class RTSimulationWorker;
-class RTSimResultsStats
-{
-  friend class RTSimulationWorker;
-protected:
-  Vector mean;
-  TwoDMatrix vcov;
-  int num_per;
-  int num_burn;
-  NormalConj nc;
-  int incomplete_simulations;
-  int thrown_periods;
-public:
-  RTSimResultsStats(int ny, int nper, int nburn = 0)
-    : mean(ny), vcov(ny, ny),
-      num_per(nper), num_burn(nburn), nc(ny),
-      incomplete_simulations(0), thrown_periods(0)
-  {
-  }
-  void simulate(int num_sim, const DecisionRule &dr, const Vector &start,
-                const TwoDMatrix &vcov, Journal &journal);
-  void simulate(int num_sim, const DecisionRule &dr, const Vector &start,
-                const TwoDMatrix &vcov);
-  void writeMat(mat_t *fd, const char *lname);
-};
-
-/* For each shock, this simulates plus and minus impulse. The class
-   maintains a vector of simulation results, each gets a particular shock
-   and sign (positive/negative). The results of type |SimResultsIRF| are
-   stored in a vector so that even ones are positive, odd ones are
-   negative.
-
-   The constructor takes a reference to the control simulations, which
-   must be finished before the constructor is called. The control
-   simulations are passed to all |SimResultsIRF|s.
-
-   The constructor also takes the vector of indices of exogenous
-   variables (|ili|) for which the IRFs are generated. The list is kept
-   (as |irf_list_ind|) for other methods. */
-
-class DynamicModel;
-class IRFResults
-{
-  std::vector<SimResultsIRF *> irf_res;
-  const DynamicModel &model;
-  std::vector<int> irf_list_ind;
-public:
-  IRFResults(const DynamicModel &mod, const DecisionRule &dr,
-             const SimResults &control, std::vector<int> ili,
-             Journal &journal);
-  ~IRFResults();
-  void writeMat(mat_t *fd, const char *prefix) const;
-};
-
-/* This worker simulates the given decision rule and inserts the result
-   to |SimResults|. */
-
-class SimulationWorker : public sthread::detach_thread
-{
-protected:
-  SimResults &res;
-  const DecisionRule &dr;
-  DecisionRule::emethod em;
-  int np;
-  const Vector &st;
-  ShockRealization &sr;
-public:
-  SimulationWorker(SimResults &sim_res,
-                   const DecisionRule &dec_rule,
-                   DecisionRule::emethod emet, int num_per,
-                   const Vector &start, ShockRealization &shock_r)
-    : res(sim_res), dr(dec_rule), em(emet), np(num_per), st(start), sr(shock_r)
-  {
-  }
-  void operator()(std::mutex &mut) override;
-};
-
-/* This worker simulates a given impulse |imp| to a given shock
-   |ishock| based on a given control simulation with index |idata|. The
-   control simulations are contained in |SimResultsIRF| which is passed
-   to the constructor. */
-
-class SimulationIRFWorker : public sthread::detach_thread
-{
-  SimResultsIRF &res;
-  const DecisionRule &dr;
-  DecisionRule::emethod em;
-  int np;
-  int idata;
-  int ishock;
-  double imp;
-public:
-  SimulationIRFWorker(SimResultsIRF &sim_res,
-                      const DecisionRule &dec_rule,
-                      DecisionRule::emethod emet, int num_per,
-                      int id, int ishck, double impulse)
-    : res(sim_res), dr(dec_rule), em(emet), np(num_per),
-      idata(id), ishock(ishck), imp(impulse)
-  {
-  }
-  void operator()(std::mutex &mut) override;
-};
-
-/* This class does the real time simulation job for
-   |RTSimResultsStats|. It simulates the model period by period. It
-   accummulates the information in the |RTSimResultsStats::nc|. If NaN or
-   Inf is observed, it ends the simulation and adds to the
-   |thrown_periods| of |RTSimResultsStats|. */
-
-class RTSimulationWorker : public sthread::detach_thread
-{
-protected:
-  RTSimResultsStats &res;
-  const DecisionRule &dr;
-  DecisionRule::emethod em;
-  int np;
-  const Vector &ystart;
-  ShockRealization &sr;
-public:
-  RTSimulationWorker(RTSimResultsStats &sim_res,
-                     const DecisionRule &dec_rule,
-                     DecisionRule::emethod emet, int num_per,
-                     const Vector &start, ShockRealization &shock_r)
-    : res(sim_res), dr(dec_rule), em(emet), np(num_per), ystart(start), sr(shock_r)
-  {
-  }
-  void operator()(std::mutex &mut) override;
-};
-
-/* This class generates draws from Gaussian distribution with zero mean
-   and the given variance-covariance matrix. It stores the factor of vcov
-   $V$ matrix, yielding $FF^T = V$. */
-
-class RandomShockRealization : virtual public ShockRealization
-{
-protected:
-  MersenneTwister mtwister;
-  TwoDMatrix factor;
-public:
-  RandomShockRealization(const ConstTwoDMatrix &v, unsigned int iseed)
-    : mtwister(iseed), factor(v.nrows(), v.nrows())
-  {
-    schurFactor(v);
-  }
-  RandomShockRealization(const RandomShockRealization &sr)
-    : mtwister(sr.mtwister), factor(sr.factor)
-  {
-  }
-  ~RandomShockRealization()
-  override = default;
-  void get(int n, Vector &out) override;
-  int
-  numShocks() const override
-  {
-    return factor.nrows();
-  }
-protected:
-  void choleskyFactor(const ConstTwoDMatrix &v);
-  void schurFactor(const ConstTwoDMatrix &v);
-};
-
-/* This is just a matrix of finite numbers. It can be constructed from
-   any |ShockRealization| with a given number of periods. */
-
-class ExplicitShockRealization : virtual public ShockRealization
-{
-  TwoDMatrix shocks;
-public:
-  ExplicitShockRealization(const ConstTwoDMatrix &sh)
-    : shocks(sh)
-  {
-  }
-  ExplicitShockRealization(const ExplicitShockRealization &sr)
-    : shocks(sr.shocks)
-  {
-  }
-  ExplicitShockRealization(ShockRealization &sr, int num_per);
-  void get(int n, Vector &out) override;
-  int
-  numShocks() const override
-  {
-    return shocks.nrows();
-  }
-  const TwoDMatrix &
-  getShocks()
-  {
-    return shocks;
-  }
-  void addToShock(int ishock, int iper, double val);
-  void
-  print() const
-  {
-    shocks.print();
-  }
-};
-
-/* This represents a user given shock realization. The first matrix of
-   the constructor is a covariance matrix of shocks, the second matrix is
-   a rectangular matrix, where columns correspond to periods, rows to
-   shocks. If an element of the matrix is {\tt NaN}, or {\tt Inf}, or
-   {\tt -Inf}, then the random shock is taken instead of that element.
-
-   In this way it is a generalization of both |RandomShockRealization|
-   and |ExplicitShockRealization|. */
-
-class GenShockRealization : public RandomShockRealization, public ExplicitShockRealization
-{
-public:
-  GenShockRealization(const ConstTwoDMatrix &v, const ConstTwoDMatrix &sh, int seed)
-    : RandomShockRealization(v, seed), ExplicitShockRealization(sh)
-  {
-    KORD_RAISE_IF(sh.nrows() != v.nrows() || v.nrows() != v.ncols(),
-                  "Wrong dimension of input matrix in GenShockRealization constructor");
-  }
-  void get(int n, Vector &out) override;
-  int
-  numShocks() const override
-  {
-    return RandomShockRealization::numShocks();
-  }
-};
-
-#endif

dynare++/kord/dynamic_model.cc

-// Copyright 2005, Ondra Kamenik
-
-#include "dynamic_model.hh"
-
-void
-NameList::print() const
-{
-  for (int i = 0; i < getNum(); i++)
-    printf("%s\n", getName(i));
-}
-
-void
-NameList::writeMat(mat_t *fd, const char *vname) const
-{
-  int maxlen = 0;
-  for (int i = 0; i < getNum(); i++)
-    if (maxlen < (int) strlen(getName(i)))
-      maxlen = (int) strlen(getName(i));
-
-  if (maxlen == 0)
-    return;
-
-  auto *m = new char[getNum()*maxlen];
-
-  for (int i = 0; i < getNum(); i++)
-    for (int j = 0; j < maxlen; j++)
-      if (j < (int) strlen(getName(i)))
-        m[j*getNum()+i] = getName(i)[j];
-      else
-        m[j*getNum()+i] = ' ';
-
-#if MATIO_MAJOR_VERSION > 1 || (MATIO_MAJOR_VERSION == 1 && MATIO_MINOR_VERSION >= 5)
-  size_t dims[2];
-  const matio_compression compression = MAT_COMPRESSION_NONE;
-#else
-  int dims[2];
-  const int compression = COMPRESSION_NONE;
-#endif
-  dims[0] = getNum();
-  dims[1] = maxlen;
-
-  matvar_t *v = Mat_VarCreate(vname, MAT_C_CHAR, MAT_T_UINT8, 2, dims, m, 0);
-
-  Mat_VarWrite(fd, v, compression);
-
-  Mat_VarFree(v);
-  delete[] m;
-}
-
-void
-NameList::writeMatIndices(mat_t *fd, const char *prefix) const
-{
-  char tmp[100];
-  TwoDMatrix aux(1, 1);
-  for (int i = 0; i < getNum(); i++)
-    {
-      sprintf(tmp, "%s_i_%s", prefix, getName(i));
-      aux.get(0, 0) = i+1;
-      aux.writeMat(fd, tmp);
-    }
-}

dynare++/kord/dynamic_model.hh

-// Copyright 2005, Ondra Kamenik
-
-// Dynamic model abstraction
-
-/* This file only defines a generic interface to an SDGE model. The model
-   takes the form:
-   $$E_t\left[f(g^{**}(g^*(y,u_t),u_{t+1}),g(y,u),y,u_t)\right]=0$$
-   The interface is defined via pure virtual class |DynamicModel|. */
-
-#ifndef DYNAMIC_MODEL_H
-#define DYNAMIC_MODEL_H
-
-#include "t_container.hh"
-#include "sparse_tensor.hh"
-
-#include "Vector.hh"
-
-/* The class is a virtual pure class which provides an access to names
-   of the variables. */
-
-class NameList
-{
-public:
-  virtual ~NameList()
-  = default;
-  virtual int getNum() const = 0;
-  virtual const char *getName(int i) const = 0;
-  void print() const;
-  void writeMat(mat_t *fd, const char *vname) const;
-  void writeMatIndices(mat_t *fd, const char *prefix) const;
-};
-
-/* This is the interface to an information on a generic SDGE
-   model. It is sufficient for calculations of policy rule Taylor
-   approximations at some (not necessarily deterministic) steady state.
-
-   We need to know a partitioning of endogenous variables $y$. We suppose
-   that $y$ is partitioned as
-   $$y=\left[\matrix{\hbox{static}\cr\hbox{pred}\cr\hbox{both}\cr\hbox{forward}}\right]$$
-   of which we define
-   $$y^*=\left[\matrix{\hbox{pred}\cr\hbox{both}}\right]\quad
-   y^{**}=\left[\matrix{\hbox{both}\cr\hbox{forward}}\right]$$
-   where ``static'' are meant those variables, which appear only at time
-   $t$; ``pred'' are meant those variables, which appear only at $t$ and
-   $t-1$; ``both'' are meant those variables, which appear at least at
-   $t-1$ and $t+1$; and ``forward'' are meant those variables, which
-   appear only at $t$ and $t+1$. This partitioning is given by methods
-   |nstat()|, |npred()|, |nboth()|, and |nforw()|. The number of
-   equations |numeq()| must be the same as a number of endogenous
-   variables.
-
-   In order to complete description, we need to know a number of
-   exogenous variables, which is a size of $u$, hence |nexog()| method.
-
-   The model contains an information about names of variables, the
-   variance-covariance matrix of the shocks, the derivatives of equations
-   of $f$ at some steady state, and the steady state. These can be
-   retrieved by the corresponding methods.
-
-   The derivatives of the system are calculated with respect to stacked
-   variables, the stack looks as:
-   $$\left[\matrix{y^{**}_{t+1}\cr y_t\cr y^*_{t-1}\cr u_t}\right].$$
-
-   There are only three operations. The first
-   |solveDeterministicSteady()| solves the deterministic steady steate
-   which can be retrieved by |getSteady()| later. The method
-   |evaluateSystem| calculates $f(y^{**},y,y^*,u)$, where $y$ and $u$ are
-   passed, or $f(y^{**}_{t+1}, y_t, y^*_{t-1}, u)$, where $y^{**}_{t+1}$,
-   $y_t$, $y^*_{t-1}$, $u$ are passed. Finally, the method
-   |calcDerivativesAtSteady()| calculates derivatives of $f$ at the
-   current steady state, and zero shocks. The derivatives can be
-   retrieved with |getModelDerivatives()|. All the derivatives are done
-   up to a given order in the model, which can be retrieved by |order()|.
-
-   The model initialization is done in a constructor of the implementing
-   class. The constructor usually calls a parser, which parses a given
-   file (usually a text file), and retrieves all necessary information
-   about the model, inluding variables, partitioning, variance-covariance
-   matrix, information helpful for calculation of the deterministic
-   steady state, and so on. */
-
-class DynamicModel
-{
-public:
-  virtual DynamicModel *clone() const = 0;
-  virtual ~DynamicModel()
-  = default;
-
-  virtual int nstat() const = 0;
-  virtual int nboth() const = 0;
-  virtual int npred() const = 0;
-  virtual int nforw() const = 0;
-  virtual int nexog() const = 0;
-  virtual int order() const = 0;
-  int
-  numeq() const
-  {
-    return nstat()+nboth()+npred()+nforw();
-  }
-
-  virtual const NameList&getAllEndoNames() const = 0;
-  virtual const NameList&getStateNames() const = 0;
-  virtual const NameList&getExogNames() const = 0;
-  virtual const TwoDMatrix&getVcov() const = 0;
-  virtual const TensorContainer<FSSparseTensor>&getModelDerivatives() const = 0;
-  virtual const Vector&getSteady() const = 0;
-  virtual Vector&getSteady() = 0;
-
-  virtual void solveDeterministicSteady() = 0;
-  virtual void evaluateSystem(Vector &out, const ConstVector &yy, const Vector &xx) = 0;
-  virtual void evaluateSystem(Vector &out, const ConstVector &yym, const ConstVector &yy,
-                              const ConstVector &yyp, const Vector &xx) = 0;
-  virtual void calcDerivativesAtSteady() = 0;
-};
-
-#endif

dynare++/kord/faa_di_bruno.cc

-// Copyright 2005, Ondra Kamenik
-
-#include "faa_di_bruno.hh"
-#include "fine_container.hh"
-
-#include <cmath>
-
-double FaaDiBruno::magic_mult = 1.5;
-
-// |FaaDiBruno::calculate| folded sparse code
-/* We take an opportunity to refine the stack container to avoid
-   allocation of more memory than available. */
-void
-FaaDiBruno::calculate(const StackContainer<FGSTensor> &cont,
-                      const TensorContainer<FSSparseTensor> &f,
-                      FGSTensor &out)
-{
-  out.zeros();
-  for (int l = 1; l <= out.dimen(); l++)
-    {
-      int mem_mb, p_size_mb;
-      int max = estimRefinment(out.getDims(), out.nrows(), l, mem_mb, p_size_mb);
-      FoldedFineContainer fine_cont(cont, max);
-      fine_cont.multAndAdd(l, f, out);
-      JournalRecord recc(journal);
-      recc << "dim=" << l << " avmem=" << mem_mb << " tmpmem=" << p_size_mb << " max=" << max
-           << " stacks=" << cont.numStacks() << "->" << fine_cont.numStacks() << endrec;
-    }
-}
-
-// |FaaDiBruno::calculate| folded dense code
-/* Here we just simply evaluate |multAndAdd| for the dense
-   container. There is no opportunity for tuning. */
-void
-FaaDiBruno::calculate(const FoldedStackContainer &cont, const FGSContainer &g,
-                      FGSTensor &out)
-{
-  out.zeros();
-  for (int l = 1; l <= out.dimen(); l++)
-    {
-      long int mem = SystemResources::availableMemory();
-      cont.multAndAdd(l, g, out);
-      JournalRecord rec(journal);
-      int mem_mb = mem/1024/1024;
-      rec << "dim=" << l << " avmem=" << mem_mb << endrec;
-    }
-}
-
-// |FaaDiBruno::calculate| unfolded sparse code
-/* This is the same as |@<|FaaDiBruno::calculate| folded sparse
-   code@>|. The only difference is that we construct unfolded fine
-   container. */
-void
-FaaDiBruno::calculate(const StackContainer<UGSTensor> &cont,
-                      const TensorContainer<FSSparseTensor> &f,
-                      UGSTensor &out)
-{
-  out.zeros();
-  for (int l = 1; l <= out.dimen(); l++)
-    {
-      int mem_mb, p_size_mb;
-      int max = estimRefinment(out.getDims(), out.nrows(), l, mem_mb, p_size_mb);
-      UnfoldedFineContainer fine_cont(cont, max);
-      fine_cont.multAndAdd(l, f, out);
-      JournalRecord recc(journal);
-      recc << "dim=" << l << " avmem=" << mem_mb << " tmpmem=" << p_size_mb << " max=" << max
-           << " stacks=" << cont.numStacks() << "->" << fine_cont.numStacks() << endrec;
-    }
-}
-
-// |FaaDiBruno::calculate| unfolded dense code
-/* Again, no tuning opportunity here. */
-void
-FaaDiBruno::calculate(const UnfoldedStackContainer &cont, const UGSContainer &g,
-                      UGSTensor &out)
-{
-  out.zeros();
-  for (int l = 1; l <= out.dimen(); l++)
-    {
-      long int mem = SystemResources::availableMemory();
-      cont.multAndAdd(l, g, out);
-      JournalRecord rec(journal);
-      int mem_mb = mem/1024/1024;
-      rec << "dim=" << l << " avmem=" << mem_mb << endrec;
-    }
-}
-
-/* This function returns a number of maximum rows used for refinement of
-   the stacked container. We want to set the maximum so that the expected
-   memory consumption for the number of paralel threads would be less
-   than available memory. On the other hand we do not want to be too
-   pesimistic since a very fine refinement can be very slow.
-
-   Besides memory needed for a dense unfolded slice of a tensor from
-   |f|, each thread needs |magic_mult*per_size| bytes of memory. In the
-   worst case, |magic_mult| will be equal to 2, this means memory
-   |per_size| for target temporary (permuted symmetry) tensor plus one
-   copy for intermediate result. However, this shows to be too
-   pesimistic, so we set |magic_mult| to 1.5. The memory for permuted
-   symmetry temporary tensor |per_size| is estimated as a weigthed
-   average of unfolded memory of the |out| tensor and unfolded memory of
-   a symetric tensor with the largest coordinate size. Some experiments
-   showed that the best combination of the two is to take 100\% if the
-   latter, so we set |lambda| to zero.
-
-   The |max| number of rows in the refined |cont| must be such that each
-   slice fits to remaining memory. Number of columns of the slice are
-   never greater $max^l$. (This is not true, since stacks corresponing to
-   unit/zero matrices cannot be further refined). We get en equation:
-
-   $$nthreads\cdot max^l\cdot 8\cdot r = mem -
-   magic\_mult\cdot nthreads\cdot per\_size\cdot 8\cdot r,$$
-   where |mem| is available memory in bytes, |nthreads| is a number of
-   threads, $r$ is a number of rows, and $8$ is |sizeof(double)|.
-
-   If the right hand side is less than zero, we set |max| to 10, just to
-   let it do something. */
-
-int
-FaaDiBruno::estimRefinment(const TensorDimens &tdims, int nr, int l,
-                           int &avmem_mb, int &tmpmem_mb)
-{
-  int nthreads = sthread::detach_thread_group::max_parallel_threads;
-  long int per_size1 = tdims.calcUnfoldMaxOffset();
-  auto per_size2 = (long int) pow((double) tdims.getNVS().getMax(), l);
-  double lambda = 0.0;
-  long int per_size = sizeof(double)*nr
-    *(long int) (lambda*per_size1+(1-lambda)*per_size2);
-  long int mem = SystemResources::availableMemory();
-  int max = 0;
-  double num_cols = ((double) (mem-magic_mult*nthreads*per_size))
-    /nthreads/sizeof(double)/nr;
-  if (num_cols > 0)
-    {
-      double maxd = pow(num_cols, ((double) 1)/l);
-      max = (int) floor(maxd);
-    }
-  if (max == 0)
-    {
-      max = 10;
-      JournalRecord rec(journal);
-      rec << "dim=" << l << " run out of memory, imposing max=" << max;
-      if (nthreads > 1)
-        rec << " (decrease number of threads)";
-      rec << endrec;
-    }
-  avmem_mb = mem/1024/1024;
-  tmpmem_mb = (nthreads*per_size)/1024/1024;
-  return max;
-}