diff --git a/dynare++/integ/cc/precalc_quadrature.hh b/dynare++/integ/cc/precalc_quadrature.hh
index 231e928d255f7b32d76e976cf6d734a5d7689792..f997e0f178cc5422b980e3d1c409daff2fd0ca15 100644
--- a/dynare++/integ/cc/precalc_quadrature.hh
+++ b/dynare++/integ/cc/precalc_quadrature.hh
@@ -18,25 +18,25 @@
* along with Dynare. If not, see <http://www.gnu.org/licenses/>.
*/
-// The file contains one dimensional quadrature points and weights for
-// a few quadratures. The format of data is clear. There is a class
-// OneDPrecalcQuadrature which implements an interface OneDQuadrature
-// using the data of this format.
+// The file contains one dimensional quadrature points and weights for a few
+// quadratures. The format of data is clear. There is a class
+// OneDPrecalcQuadrature which implements an interface OneDQuadrature using the
+// data of this format.
// Gauss-Hermite quadrature; prefix gh
-// number of levels
+// Number of levels
static const int gh_num_levels = 26;
-// number of points in each level
+// Number of points in each level
static const int gh_num_points[] = {
1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20,
30, 32, 40, 50, 60, 64
};
-// weights, starting with the first level
+// Weights, starting with the first level
static const double gh_weights[] = {
- // weights 1 = sqrt(pi)
+ // weights 1 = √π
1.77245385090551588191942755656782537698745727539062,
// weights 2
0.886226925452758013649083741671e+00,
@@ -550,7 +550,7 @@ static const double gh_weights[] = {
0.553570653584e-48
};
-// points, starting with the first level
+// Points, starting with the first level
static const double gh_points[] = {
// points 1
0.0,
@@ -1068,16 +1068,16 @@ static const double gh_points[] = {
// Gauss-Legendre quadrature; prefix gl
-// number of levels
+// Number of levels
static const int gl_num_levels = 22;
-// number of points in each level
+// Number of points in each level
static const int gl_num_points[] = {
1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20,
32, 64
};
-// weights, starting with the first level
+// Weights, starting with the first level
static const double gl_weights[] = {
// weight 1
2.0e+00,
@@ -1409,7 +1409,7 @@ static const double gl_weights[] = {
0.178328072169643294729607914497e-02
};
-// points, starting with the first level
+// Points, starting with the first level
static const double gl_points[] = {
// points 1
0.0e+00,
diff --git a/dynare++/integ/cc/product.cc b/dynare++/integ/cc/product.cc
index 711abb32c2680a316b350f9a4b281f159c4fdf25..d3515455087cf07895c16458863b1fea30fa1e96 100644
--- a/dynare++/integ/cc/product.cc
+++ b/dynare++/integ/cc/product.cc
@@ -24,7 +24,7 @@
#include <iostream>
#include <iomanip>
-/* This constructs a product iterator corresponding to index $(j0,0\ldots,0)$. */
+/* This constructs a product iterator corresponding to index (j0,0,…,0). */
prodpit::prodpit(const ProductQuadrature &q, int j0, int l)
: prodq(q), level(l), npoints(q.uquad.numPoints(l)),
@@ -51,7 +51,6 @@ prodpit::operator==(const prodpit &ppit) const
prodpit &
prodpit::operator++()
{
- // todo: throw if |prodq==NULL| or |jseq==NULL| or |sig==NULL| or |end_flag==true|
int i = prodq.dimen()-1;
jseq[i]++;
while (i >= 0 && jseq[i] == npoints)
@@ -77,8 +76,6 @@ prodpit::operator++()
void
prodpit::setPointAndWeight()
{
- // todo: raise if |prodq==NULL| or |jseq==NULL| or |sig==NULL| or
- // |p==NULL| or |end_flag==true|
w = 1.0;
for (int i = 0; i < prodq.dimen(); i++)
{
@@ -107,26 +104,26 @@ prodpit::print() const
ProductQuadrature::ProductQuadrature(int d, const OneDQuadrature &uq)
: QuadratureImpl<prodpit>(d), uquad(uq)
{
- // todo: check |d>=1|
+ // TODO: check d≥1
}
-/* This calls |prodpit| constructor to return an iterator which points
- approximatelly at |ti|-th portion out of |tn| portions. First we find
+/* This calls prodpit constructor to return an iterator which points
+ approximatelly at ‘ti’-th portion out of ‘tn’ portions. First we find
out how many points are in the level, and then construct an interator
- $(j0,0,\ldots,0)$ where $j0=$|ti*npoints/tn|. */
+ (j0,0,…,0) where j0=ti·npoints/tn. */
prodpit
ProductQuadrature::begin(int ti, int tn, int l) const
{
- // todo: raise is |l<dimen()|
- // todo: check |l<=uquad.numLevels()|
+ // TODO: raise if l<dimen()
+ // TODO: check l ≤ uquad.numLevels()
int npoints = uquad.numPoints(l);
return prodpit(*this, ti*npoints/tn, l);
}
-/* This just starts at the first level and goes to a higher level as
- long as a number of evaluations (which is $n_k^d$ for $k$ being the
- level) is less than the given number of evaluations. */
+/* This just starts at the first level and goes to a higher level as long as a
+ number of evaluations (which is nₖᵈ for k being the level) is less than the
+ given number of evaluations. */
void
ProductQuadrature::designLevelForEvals(int max_evals, int &lev, int &evals) const
@@ -143,5 +140,4 @@ ProductQuadrature::designLevelForEvals(int max_evals, int &lev, int &evals) cons
while (lev < uquad.numLevels()-2 && evals < max_evals);
lev--;
evals = last_evals;
-
}
diff --git a/dynare++/integ/cc/product.hh b/dynare++/integ/cc/product.hh
index 512580aac1a872da2146b9945032d7031d7230ae..894cc7add400a153e66ebc9c8bc5bca594f3345d 100644
--- a/dynare++/integ/cc/product.hh
+++ b/dynare++/integ/cc/product.hh
@@ -20,16 +20,20 @@
// 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$$
+/* This file defines a product multidimensional quadrature. If Qₖ$ denotes the
+ one dimensional quadrature, then the product quadrature Q of k level and
+ dimension d takes the form:
- Here we define the product quadrature iterator |prodpit| and plug it
- into |QuadratureImpl| to obtains |ProductQuadrature|. */
+ nₖ nₖ
+ Qf = ∑ … ∑ w_i₁·…·w_{i_d} f(x_i₁,…,x_{i_d})
+ i₁=1 i_d=1
+
+ which can be written in terms of the one dimensional quadrature Qₖ as
+
+ Qf=(Qₖ⊗…⊗Qₖ)f
+
+ Here we define the product quadrature iterator prodpit and plug it into
+ QuadratureImpl to obtains ProductQuadrature. */
#ifndef PRODUCT_H
#define PRODUCT_H
@@ -38,18 +42,18 @@
#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.
+/* 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 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|. */
+ The constructor prodpit(const ProductQuadrature& q, int j0, int l)
+ constructs an iterator pointing to (j0,0,…,0), which is used by begin()
+ dictated by QuadratureImpl. */
class ProductQuadrature;
@@ -97,13 +101,12 @@ 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. */
+/* 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>
{
diff --git a/dynare++/integ/cc/quadrature.hh b/dynare++/integ/cc/quadrature.hh
index dcbc3d41ea27ec078b962c5c7f1f0479e8a1c29e..7c6cb8e2c4221b028ae4c5f67f181e8bb04a2be7 100644
--- a/dynare++/integ/cc/quadrature.hh
+++ b/dynare++/integ/cc/quadrature.hh
@@ -21,29 +21,29 @@
// 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. */
+ OneDQuadrature, and a parent for all multi-dimensional quadratures. This
+ parent class Quadrature presents a general concept of quadrature, this is
+
+ N
+ ∫ f(x)dx ≃ ∑ wᵢxᵢ
+ ⁱ⁼¹
+
+ 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ᵢ and wᵢ for i=1,…,N. All the iterators must be able to go through only a
+ portion of the set i=1,…,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
@@ -56,10 +56,10 @@
#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. */
+/* 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
{
@@ -71,18 +71,17 @@ public:
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|.
+/* 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. */
+ 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
{
@@ -104,11 +103,11 @@ public:
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.
+/* 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. */
+ See QuadratureImpl class declaration for details. */
template <typename _Tpit>
class QuadratureImpl;
@@ -129,14 +128,13 @@ public:
{
}
- /* 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 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. */
+ This method just everything up as it is coming. This might be imply large
+ numerical errors, perhaps in future something smarter should be implemented. */
void
operator()(std::mutex &mut) override
@@ -147,8 +145,8 @@ public:
tmpall.zeros();
Vector tmp(outvec.length());
- // note that since beg came from begin, it has empty signal
- // and first evaluation gets no signal
+ /* 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);
@@ -162,13 +160,13 @@ public:
}
};
-/* 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.
+/* 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. */
+ 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
@@ -183,8 +181,8 @@ public:
void
integrate(VectorFunctionSet &fs, int level, Vector &out) const override
{
- // todo: out.length()==func.outdim()
- // todo: dim == func.indim()
+ // TODO: out.length()==func.outdim()
+ // TODO: dim == func.indim()
out.zeros();
sthread::detach_thread_group gr;
for (int ti = 0; ti < fs.getNum(); ti++)
@@ -207,7 +205,7 @@ public:
std::ofstream fd{fname, std::ios::out | std::ios::trunc};
if (fd.fail())
{
- // todo: raise
+ // TODO: raise
std::cerr << "Cannot open file " << fname << " for writing." << std::endl;
exit(EXIT_FAILURE);
}
@@ -238,16 +236,16 @@ 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.
+/* This is only an interface to a precalculated data in file
+ precalc_quadrature.hh 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. */
+ The implementing subclasses just fill the necessary data from the file, the
+ rest is calculated here. */
class OneDPrecalcQuadrature : public OneDQuadrature
{
@@ -289,17 +287,21 @@ 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$.
+/* Just precalculated Gauss-Hermite quadrature. This quadrature integrates
+ exactly integrals
- 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.
+ +∞
+ ∫ xᵏe^{−x²}dx
+ −∞
- Check {\tt precalc\_quadrature.dat} for available levels. */
+ 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 precalc_quadrature.hh for available levels. */
class GaussHermite : public OneDPrecalcQuadrature
{
@@ -307,11 +309,16 @@ public:
GaussHermite();
};
-/* Just precalculated Gauss--Legendre quadrature. This quadrature integrates exactly integrals
- $$\int_0^1x^k{\rm d}x$$
- for level $k$.
+/* Just precalculated Gauss-Legendre quadrature. This quadrature integrates
+ exactly integrals
+
+ ₁
+ ∫ xᵏdx
+ ⁰
+
+ for level k.
- Check {\tt precalc\_quadrature.dat} for available levels. */
+ Check precalc_quadrature.hh for available levels. */
class GaussLegendre : public OneDPrecalcQuadrature
{
diff --git a/dynare++/integ/cc/quasi_mcarlo.cc b/dynare++/integ/cc/quasi_mcarlo.cc
index 2da9d99450cc6fe952f39be1ab8004af0e7d5b57..fef6bfdf1199df695addcee4af6eda66f9060a39 100644
--- a/dynare++/integ/cc/quasi_mcarlo.cc
+++ b/dynare++/integ/cc/quasi_mcarlo.cc
@@ -25,9 +25,9 @@
#include <iomanip>
#include <array>
-/* 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. */
+/* 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),
@@ -44,23 +44,26 @@ RadicalInverse::RadicalInverse(int n, int b, int mxn)
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}
- $$ */
+/* This evaluates the radical inverse. If there was no permutation, we have to
+ calculate:
+
+ c₀ c₁ cⱼ
+ ── + ── + … + ────
+ b b² bʲ⁺¹
+
+ which is evaluated as:
+
+ ⎛ ⎛⎛cⱼ 1 cⱼ₋₁⎞ 1 cⱼ₋₂⎞ ⎞ 1 c₀
+ ⎢…⎢⎢──·─ + ────⎥·─ + ────⎥…⎥·─ + ──
+ ⎝ ⎝⎝ b b b ⎠ b b ⎠ ⎠ b b
+
+
+ Now with permutation π, we have:
+
+ ⎛ ⎛⎛π(cⱼ) 1 π(cⱼ₋₁)⎞ 1 π(cⱼ₋₂)⎞ ⎞ 1 π(c₀)
+ ⎢…⎢⎢─────·─ + ───────⎥·─ + ───────⎥…⎥·─ + ─────
+ ⎝ ⎝⎝ b b b ⎠ b b ⎠ ⎠ b b
+*/
double
RadicalInverse::eval(const PermutationScheme &p) const
@@ -79,7 +82,7 @@ RadicalInverse::eval(const PermutationScheme &p) const
void
RadicalInverse::increase()
{
- // todo: raise if |num+1 > maxn|
+ // TODO: raise if num+1 > maxn
num++;
int i = 0;
coeff[i]++;
@@ -100,8 +103,8 @@ RadicalInverse::print() const
coeff.print();
}
-/* Here we have the first 170 primes. This means that we are not able
- to integrate dimensions greater than 170. */
+/* 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,
@@ -123,21 +126,21 @@ std::array<int, 170> HaltonSequence::primes = {
947, 953, 967, 971, 977, 983, 991, 997, 1009, 1013
};
-/* This takes first |dim| primes and constructs |dim| radical inverses
- and calls |eval|. */
+/* 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|
+ // 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|. */
+/* This calls RadicalInverse::increase() for all radical inverses and calls
+ eval(). */
void
HaltonSequence::increase()
@@ -149,7 +152,7 @@ HaltonSequence::increase()
eval();
}
-/* This sets point |pt| to radical inverse evaluations in each dimension. */
+/* This sets point ‘pt’ to radical inverse evaluations in each dimension. */
void
HaltonSequence::eval()
{
@@ -187,7 +190,6 @@ qmcpit::operator==(const qmcpit &qpit) const
qmcpit &
qmcpit::operator++()
{
- // todo: raise if |halton == null || qmcq == NULL|
halton.increase();
return *this;
}
@@ -198,14 +200,12 @@ 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
{
diff --git a/dynare++/integ/cc/quasi_mcarlo.hh b/dynare++/integ/cc/quasi_mcarlo.hh
index c423073c5b7fac731a1d70470366a14df57169d0..8df3e63fa550ad13b6d9b60c344f3777fa924d78 100644
--- a/dynare++/integ/cc/quasi_mcarlo.hh
+++ b/dynare++/integ/cc/quasi_mcarlo.hh
@@ -22,23 +22,29 @@
/* 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|. */
+ QMCarloCubeQuadrature and integrates:
+
+ ∫ f(x)dx
+ [0,1]ⁿ
+
+ The quadrature for a function of normally distributed parameters is named
+ QMCarloNormalQuadrature and integrates:
+
+ 1
+ ──────── ∫ f(x)e^{−½xᵀx}dx
+ √{(2π)ⁿ} [−∞,+∞]ⁿ
+
+ 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
@@ -50,9 +56,9 @@
#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$. */
+/* 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,…,base−1. */
class PermutationScheme
{
@@ -62,15 +68,14 @@ public:
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|.
+/* This class represents an integer number ‘num’ as c₀+c₁b+c₂b²+…+cⱼbʲ, where b
+ is ‘base’ and c₀,…,cⱼ 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. */
+ 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
{
@@ -88,11 +93,11 @@ public:
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. */
+/* This is a vector of RadicalInverses, 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
{
@@ -124,10 +129,9 @@ 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. */
+/* 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
{
@@ -158,10 +162,10 @@ public:
}
};
-/* 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. */
+/* 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
{
@@ -199,10 +203,9 @@ public:
}
};
-/* 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|. */
+/* 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
{
diff --git a/dynare++/integ/cc/smolyak.cc b/dynare++/integ/cc/smolyak.cc
index 067cff5969bbdb2798dcb17b2a9ce137bffce433..0a41f0645c177996b1a17fb402e1c290749fc850 100644
--- a/dynare++/integ/cc/smolyak.cc
+++ b/dynare++/integ/cc/smolyak.cc
@@ -24,9 +24,9 @@
#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. */
+/* 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),
@@ -44,14 +44,13 @@ 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. */
+/* 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]++;
@@ -73,14 +72,13 @@ smolpit::operator++()
return *this;
}
-/* Here we set the point coordinates according to |jseq| and
- |isummand|. Also the weight is set here. */
+/* 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()|
+ // todo: raise if isummand ≥ smolq.numSummands()
int l = smolq.level;
int d = smolq.dimen();
int sumk = (smolq.levels[isummand]).sum();
@@ -114,21 +112,19 @@ smolpit::print() const
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$. */
+/* 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≤|k|≤l+d−1 and all kᵢ are positive integers. This is
+ equivalent to going through all k such that l−d≤|k|≤l−1 and all kᵢ are
+ non-negative integers. This is equivalent to going through d+1 dimensional
+ sequences (k,x) such that |(k,x)|=l−1 and x=0,…,d−1. The resulting sequence
+ of positive integers is obtained by adding 1 to all kᵢ. */
SmolyakQuadrature::SmolyakQuadrature(int d, int l, const OneDQuadrature &uq)
: QuadratureImpl<smolpit>(d), level(l), uquad(uq)
{
- // todo: check |l>1|, |l>=d|
- // todo: check |l>=uquad.miLevel()|, |l<=uquad.maxLevel()|
+ // TODO: check l>1, l≥d
+ // TODO: check l≥uquad.miLevel(), l≤uquad.maxLevel()
int cum = 0;
for (const auto &si : SymmetrySet(l-1, d+1))
{
@@ -147,10 +143,10 @@ SmolyakQuadrature::SmolyakQuadrature(int d, int l, const OneDQuadrature &uq)
}
}
-/* 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. */
+/* 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
@@ -161,14 +157,14 @@ SmolyakQuadrature::numEvals(int l) const
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|. */
+/* 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|
+ // TODO: raise is level≠l
if (ti == tn)
return smolpit(*this, numSummands());
@@ -180,9 +176,9 @@ SmolyakQuadrature::begin(int ti, int tn, int l) const
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. */
+/* 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
@@ -203,8 +199,8 @@ SmolyakQuadrature::calcNumEvaluations(int lev) const
return cum;
}
-/* This returns a maximum level such that the number of evaluations is
- less than the given number. */
+/* 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
diff --git a/dynare++/integ/cc/smolyak.hh b/dynare++/integ/cc/smolyak.hh
index dc88523612c130fda81d2a5aff9d62d7b9ed3777..a51a98bb75c7b869ca9cc3076a55bf0717008da1 100644
--- a/dynare++/integ/cc/smolyak.hh
+++ b/dynare++/integ/cc/smolyak.hh
@@ -20,17 +20,19 @@
// 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|. */
+/* This file defines Smolyak (sparse grid) multidimensional quadrature for
+ non-nested underlying one dimensional quadrature. Let Q¹ₗ denote the one
+ dimensional quadrature of l level. Let nₗ denote a number of points in the l
+ level. Than the Smolyak quadrature can be defined as
+
+ ⎛ d−1 ⎞
+ Qᵈf = ∑ (−1)^{l+d−|k|−1} ⎢ ⎥ (Q¹_k₁⊗…⊗Q¹_{k_d})f
+ l≤|k|≤l+d−1 ⎝|k|−l⎠
+
+ where d is the dimension, k is d-dimensional sequence of integers, and |k|
+ denotes the sum of the sequence.
+
+ Here we define smolpit as Smolyak iterator and SmolyakQuadrature. */
#ifndef SMOLYAK_H
#define SMOLYAK_H
@@ -41,18 +43,16 @@
#include "quadrature.hh"
#include "pascal_triangle.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|.
+/* 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. */
+ 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;
@@ -97,23 +97,30 @@ 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|. */
+/* Here we define the class SmolyakQuadrature. It maintains an array of
+ summands of the Smolyak quadrature formula:
+
+ ⎛ d−1 ⎞
+ ∑ (−1)^{l+d−|k|−1} ⎢ ⎥ (Q¹_k₁⊗…⊗Q¹_{k_d})f
+ l≤|k|≤l+d−1 ⎝|k|−l⎠
+
+ 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:
+
+ d
+ ∑ ∏ n_kᵢ
+ ᵏ ⁱ⁼¹
+
+ 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>
{
diff --git a/dynare++/integ/cc/vector_function.cc b/dynare++/integ/cc/vector_function.cc
index e92912a8a541e7256e849bb5666a26fa37bfa557..bca1b1afabc8f1213b94a354272df434b293eb58 100644
--- a/dynare++/integ/cc/vector_function.cc
+++ b/dynare++/integ/cc/vector_function.cc
@@ -25,16 +25,16 @@
#include <cmath>
#include <algorithm>
-/* Just an easy constructor of sequence of booleans defaulting to
- change everywhere. */
+/* 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). */
+/* This sets ‘false’ (no change) before a given parameter, and ‘true’ (change)
+ after the given parameter (including). */
void
ParameterSignal::signalAfter(int l)
@@ -56,8 +56,8 @@ VectorFunctionSet::VectorFunctionSet(const VectorFunction &f, int n)
}
}
-/* This constructs a function set with shallow copy in the first and
- hard copies in others. */
+/* This constructs a function set with shallow copy in the first and hard
+ copies in others. */
VectorFunctionSet::VectorFunctionSet(VectorFunction &f, int n)
: funcs(n)
{
@@ -71,25 +71,23 @@ VectorFunctionSet::VectorFunctionSet(VectorFunction &f, int n)
}
-/* 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$. */
+/* Here we construct the object from the given function f and given
+ variance-covariance matrix Σ=vcov. The matrix A is calculated as lower
+ triangular and yields Σ=AAᵀ. */
GaussConverterFunction::GaussConverterFunction(VectorFunction &f, const GeneralMatrix &vcov)
: VectorFunction(f), func(&f), A(vcov.nrows(), vcov.nrows()),
multiplier(calcMultiplier())
{
- // todo: raise if |A.nrows() != indim()|
+ // TODO: raise if A.nrows() ≠ 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.nrows(), vcov.nrows()),
multiplier(calcMultiplier())
{
- // todo: raise if |A.nrows() != indim()|
+ // TODO: raise if A.nrows() ≠ indim()
calcCholeskyFactor(vcov);
}
@@ -100,11 +98,12 @@ GaussConverterFunction::GaussConverterFunction(const GaussConverterFunction &f)
}
/* 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|. */
+
+ g(y) = 1/√(πⁿ) f(√2·Ay).
+
+ Since the matrix A is lower triangular, the change signal for the function f
+ will look like (0,…,0,1,…,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)
@@ -125,7 +124,7 @@ GaussConverterFunction::eval(const Vector &point, const ParameterSignal &sig, Ve
out.mult(multiplier);
}
-/* This returns $1\over\sqrt{\pi^n}$. */
+/* This returns 1/√(πⁿ). */
double
GaussConverterFunction::calcMultiplier() const
@@ -145,5 +144,5 @@ GaussConverterFunction::calcCholeskyFactor(const GeneralMatrix &vcov)
lapack_int info;
dpotrf("L", &rows, A.base(), &lda, &info);
- // todo: raise if |info!=1|
+ // TODO: raise if info≠1
}
diff --git a/dynare++/integ/cc/vector_function.hh b/dynare++/integ/cc/vector_function.hh
index 0ab061ecc085f793695c8ad84714bdccbf6c00e3..6a0f7dec6ce920de0d876a04743298cc9dc604f4 100644
--- a/dynare++/integ/cc/vector_function.hh
+++ b/dynare++/integ/cc/vector_function.hh
@@ -20,16 +20,15 @@
// 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.
+/* 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|. */
+ From the signalling reason, and from other reasons, the function evaluation
+ is not const. */
#ifndef VECTOR_FUNCTION_H
#define VECTOR_FUNCTION_H
@@ -40,12 +39,11 @@
#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.
+/* 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. */
+ ‘true’ means that a parameter was changed. */
class ParameterSignal
{
@@ -68,13 +66,12 @@ public:
}
};
-/* 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.
+/* 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. */
+ 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
{
@@ -102,13 +99,13 @@ public:
}
};
-/* 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.
+/* 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). */
+ 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
{
@@ -133,27 +130,34 @@ public:
}
};
-/* 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
+/* This class wraps another VectorFunction to allow integration of a function
+ through normally distributed inputs. Namely, if one wants to integrate
+
+ 1
+ ─────────── ∫ f(x)e^{−½xᵀ|Σ|⁻¹x}dx
+ √{(2π)ⁿ|Σ|}
+
+ then if we write Σ=AAᵀ and x=√2·Ay, we get integral
+
+ 1 1
+ ─────────── ∫ f(√2·Ay)e^{−½yᵀy} √(2ⁿ)|A|dy = ───── ∫ f(√2·Ay)e^{−½yᵀy}dy
+ √{(2π)ⁿ|Σ|} √(πⁿ)
+
+ which means that a given function f we have to wrap to yield a function
+
+ g(y) = 1/√(πⁿ) f(√2·Ay).
+
+ 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²}, 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$. */
+ the function and in that case the class takes ownership of f. */
class GaussConverterFunction : public VectorFunction
{
diff --git a/dynare++/integ/src/quadrature-points.cc b/dynare++/integ/src/quadrature-points.cc
index fda5360b9a512f610b055d5d5d2f1de6135c2afa..0e62fef04b6a77c25641155cd712d227915246c5 100644
--- a/dynare++/integ/src/quadrature-points.cc
+++ b/dynare++/integ/src/quadrature-points.cc
@@ -54,7 +54,7 @@ QuadParams::QuadParams(int argc, char **argv)
{
if (argc == 1)
{
- // print the help and exit
+ // 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);
}
@@ -131,7 +131,7 @@ main(int argc, char **argv)
{
QuadParams params(argc, argv);
- // open output file for writing
+ // Open output file for writing
std::ofstream fout{params.outname, std::ios::out | std::ios::trunc};
if (fout.fail())
{
@@ -146,26 +146,26 @@ main(int argc, char **argv)
buffer << f.rdbuf();
std::string contents{buffer.str()};
- // parse the vcov matrix
+ // Parse the vcov matrix
ogp::MatrixParser mp;
mp.parse(contents);
if (mp.nrows() != mp.ncols())
throw ogu::Exception(__FILE__, __LINE__,
"VCOV matrix not square");
- // and put to the GeneralMatrix
+ // 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
+ // Calculate the factor A of vcov, so that A·Aᵀ=VCOV
GeneralMatrix A(vcov.nrows(), vcov.nrows());
SymSchurDecomp ssd(vcov);
ssd.getFactor(A);
- // construct Gauss-Hermite quadrature
+ // Construct Gauss-Hermite quadrature
GaussHermite ghq;
- // construct Smolyak quadrature
+ // Construct Smolyak quadrature
int level = params.max_level;
SmolyakQuadrature sq(vcov.nrows(), level, ghq);
@@ -173,11 +173,11 @@ main(int argc, char **argv)
<< "Maximum level: " << level << std::endl
<< "Total number of nodes: " << sq.numEvals(level) << std::endl;
- // put the points to the vector
+ // 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_cast<const Vector &>(qit.point())));
- // sort and uniq
+ // 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());
@@ -185,7 +185,7 @@ main(int argc, char **argv)
std::cout << "Duplicit nodes removed: " << static_cast<unsigned long>(sq.numEvals(level)-points.size())
<< std::endl;
- // calculate weights and mass
+ // Calculate weights and mass
double mass = 0.0;
std::vector<double> weights;
for (auto & point : points)
@@ -194,7 +194,7 @@ main(int argc, char **argv)
mass += weights.back();
}
- // calculate discarded mass
+ // Calculate discarded mass
double discard_mass = 0.0;
for (double weight : weights)
if (weight/mass < params.discard_weight)
@@ -202,7 +202,7 @@ main(int argc, char **argv)
std::cout << "Total mass discarded: " << std::fixed << discard_mass/mass << std::endl;
- // dump the results
+ // Dump the results
int npoints = 0;
double upscale_weight = 1/(mass-discard_mass);
Vector x(vcov.nrows());
@@ -210,11 +210,11 @@ main(int argc, char **argv)
for (int i = 0; i < static_cast<int>(weights.size()); i++)
if (weights[i]/mass >= params.discard_weight)
{
- // print the upscaled weight
+ // Print the upscaled weight
fout << std::setw(20) << upscale_weight*weights[i];
- // multiply point with the factor A and sqrt(2)
+ // Multiply point with the factor A and √2
A.multVec(0.0, x, std::sqrt(2.), *(points[i]));
- // print the coordinates
+ // Print the coordinates
for (int j = 0; j < x.length(); j++)
fout << ' ' << std::setw(20) << x[j];
fout << std::endl;
diff --git a/dynare++/integ/testing/tests.cc b/dynare++/integ/testing/tests.cc
index 664967e49489d93a884dab9a7fa7318cd6eddf85..030fe5b1a7ad6a4ade4c706a2b7500ef9d555a58 100644
--- a/dynare++/integ/testing/tests.cc
+++ b/dynare++/integ/testing/tests.cc
@@ -40,8 +40,8 @@
#include <memory>
#include <cstdlib>
-// evaluates unfolded (Dx)^k power, where x is a vector, D is a
-// Cholesky factor (lower triangular)
+/* Evaluates unfolded (Dx)ᵏ power, where x is a vector, D is a Cholesky factor
+ (lower triangular) */
class MomentFunction : public VectorFunction
{
GeneralMatrix D;
@@ -109,9 +109,8 @@ TensorPower::eval(const Vector &point, const ParameterSignal &sig, Vector &out)
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
+/* Evaluates (1+1/d)ᵈ(x₁·…·x_d)^(1/d), its integral over [0,1]ᵈ
+ is 1.0, and its variation grows exponentially */
class Function1 : public VectorFunction
{
int dim;
@@ -148,8 +147,8 @@ Function1::eval(const Vector &point, const ParameterSignal &sig, Vector &out)
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
+// Evaluates Function1 but with transformation xᵢ=0.5(yᵢ+1)
+// This makes the new function integrate over [−1,1]ᵈ to 1.0
class Function1Trans : public Function1
{
public:
@@ -177,9 +176,8 @@ Function1Trans::eval(const Vector &point, const ParameterSignal &sig, Vector &ou
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
+/* 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;
@@ -208,7 +206,7 @@ class TestRunnable
public:
const std::string name;
int dim; // dimension of the solved problem
- int nvar; // number of variable of the solved problem
+ int nvar; // number of variables of the solved problem
TestRunnable(std::string name_arg, int d, int nv)
: name{move(name_arg)}, dim(d), nvar(nv)
{
@@ -252,15 +250,15 @@ TestRunnable::test() const
bool
TestRunnable::smolyak_normal_moments(const GeneralMatrix &m, int imom, int level)
{
- // first make m*m' and then Cholesky factor
+ // First make m·mᵀ and then Cholesky factor
GeneralMatrix msq(m * transpose(m));
- // make vector function
+ // Make vector function
int dim = m.nrows();
TensorPower tp(dim, imom);
GaussConverterFunction func(tp, msq);
- // smolyak quadrature
+ // Smolyak quadrature
Vector smol_out(UFSTensor::calcMaxOffset(dim, imom));
{
WallTimer tim("\tSmolyak quadrature time: ");
@@ -270,7 +268,7 @@ TestRunnable::smolyak_normal_moments(const GeneralMatrix &m, int imom, int level
std::cout << "\tNumber of Smolyak evaluations: " << quad.numEvals(level) << std::endl;
}
- // check against theoretical moments
+ // 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;
@@ -280,15 +278,15 @@ TestRunnable::smolyak_normal_moments(const GeneralMatrix &m, int imom, int level
bool
TestRunnable::product_normal_moments(const GeneralMatrix &m, int imom, int level)
{
- // first make m*m' and then Cholesky factor
+ // First make m·mᵀ and then Cholesky factor
GeneralMatrix msq(m * transpose(m));
- // make vector function
+ // Make vector function
int dim = m.nrows();
TensorPower tp(dim, imom);
GaussConverterFunction func(tp, msq);
- // product quadrature
+ // Product quadrature
Vector prod_out(UFSTensor::calcMaxOffset(dim, imom));
{
WallTimer tim("\tProduct quadrature time: ");
@@ -298,7 +296,7 @@ TestRunnable::product_normal_moments(const GeneralMatrix &m, int imom, int level
std::cout << "\tNumber of product evaluations: " << quad.numEvals(level) << std::endl;
}
- // check against theoretical moments
+ // 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;
@@ -350,7 +348,6 @@ TestRunnable::qmc_cube(const VectorFunction &func, double res, double tol, int l
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;
@@ -360,7 +357,6 @@ TestRunnable::qmc_cube(const VectorFunction &func, double res, double tol, int l
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;
@@ -370,7 +366,6 @@ TestRunnable::qmc_cube(const VectorFunction &func, double res, double tol, int l
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;
@@ -466,7 +461,7 @@ public:
}
};
-// note that here we pass 1,1 to tls since smolyak has its own PascalTriangle
+// Note that here we pass 1,1 to tls since smolyak has its own PascalTriangle
class F1GaussLegendre : public TestRunnable
{
public:
@@ -505,7 +500,7 @@ int
main()
{
std::vector<std::unique_ptr<TestRunnable>> all_tests;
- // fill in vector of 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>());
@@ -513,7 +508,7 @@ main()
all_tests.push_back(std::make_unique<F1GaussLegendre>());
all_tests.push_back(std::make_unique<F1QuasiMCarlo>());
- // find maximum dimension and maximum nvar
+ // Find maximum dimension and maximum nvar
int dmax = 0;
int nvmax = 0;
for (const auto &test : all_tests)
@@ -523,7 +518,7 @@ main()
}
TLStatic::init(dmax, nvmax); // initialize library
- // launch the tests
+ // Launch the tests
int success = 0;
for (const auto &test : all_tests)
{
diff --git a/dynare++/src/dynare3.cc b/dynare++/src/dynare3.cc
index 226570eaabdb9fc8a2d43c381e50398b3b67c383..3b220cfbcf93d5963e1e00b95288ad610e5ed008 100644
--- a/dynare++/src/dynare3.cc
+++ b/dynare++/src/dynare3.cc
@@ -179,7 +179,7 @@ Dynare::solveDeterministicSteady(Vector &steady)
"Could not obtain convergence in non-linear solver");
}
-// evaluate system at given y_t=y_{t+1}=y_{t-1}, and given shocks x_t
+// Evaluate system at given yₜ=yₜ₊₁=yₜ₋₁, and given shocks xₜ
void
Dynare::evaluateSystem(Vector &out, const ConstVector &yy, const Vector &xx)
{
@@ -188,9 +188,8 @@ Dynare::evaluateSystem(Vector &out, const ConstVector &yy, const Vector &xx)
evaluateSystem(out, yym, yy, yyp, xx);
}
-// evaluate system at given y^*_{t-1}, y_t, y^{**}_{t+1} and at
-// exogenous x_t, all three vectors yym, yy, and yyp have the
-// respective lengths of y^*_{t-1}, y_t, y^{**}_{t+1}
+/* Evaluate system at given y*ₜ₋₁, yₜ, y**ₜ₊₁ and at exogenous xₜ, all three
+ vectors yym, yy, and yyp have the respective lengths of y*ₜ₋₁, yₜ, y**ₜ₊₁ */
void
Dynare::evaluateSystem(Vector &out, const ConstVector &yym, const ConstVector &yy,
const ConstVector &yyp, const Vector &xx)
@@ -317,8 +316,8 @@ DynareEvalLoader::DynareEvalLoader(const ogp::FineAtoms &a, Vector &out)
throw DynareException(__FILE__, __LINE__, "Wrong length of out vector in DynareEvalLoader constructor");
}
-/** This clears the container of model derivatives and initializes it
- * inserting empty sparse tensors up to the given order. */
+/* This clears the container of model derivatives and initializes it inserting
+ empty sparse tensors up to the given order. */
DynareDerEvalLoader::DynareDerEvalLoader(const ogp::FineAtoms &a,
TensorContainer<FSSparseTensor> &mod_ders,
int order)
diff --git a/dynare++/src/dynare3.hh b/dynare++/src/dynare3.hh
index e2aede483f0cad116949ba7d57f9278a072df2ad..4990a91f33e7258815e5f23442ba5107dd4d92b8 100644
--- a/dynare++/src/dynare3.hh
+++ b/dynare++/src/dynare3.hh
@@ -51,9 +51,9 @@ public:
{
return names[i];
}
- /** This for each string of the input vector calculates its index
- * in the names. And returns the resulting vector of indices. If
- * the name cannot be found, then an exception is raised. */
+ /* This for each string of the input vector calculates its index in the
+ names. And returns the resulting vector of indices. If the name cannot be
+ found, then an exception is raised. */
std::vector<int> selectIndices(const std::vector<std::string> &ns) const;
};
@@ -112,8 +112,8 @@ class Dynare : public DynamicModel
std::unique_ptr<ogp::FormulaDerEvaluator> fde;
const double ss_tol;
public:
- /** Parses the given model file and uses the given order to
- * override order from the model file (if it is != -1). */
+ /* Parses the given model file and uses the given order to
+ override order from the model file (if it is ≠ −1). */
Dynare(const std::string &modname, int ord, double sstol, Journal &jr);
/** Parses the given equations with explicitly given names. */
Dynare(const std::vector<std::string> &endo,
@@ -121,7 +121,7 @@ public:
const std::vector<std::string> &par,
const std::string &equations, int ord,
double sstol, Journal &jr);
- /** Makes a deep copy of the object. */
+ /* Makes a deep copy of the object. */
Dynare(const Dynare &dyn);
Dynare(Dynare &&) = default;
std::unique_ptr<DynamicModel>
diff --git a/dynare++/src/dynare_atoms.hh b/dynare++/src/dynare_atoms.hh
index f389fe2d289fd3a607d5b4e650cc3702b87b285b..c9639be182db47c3355d6a9fa00c81e355955df1 100644
--- a/dynare++/src/dynare_atoms.hh
+++ b/dynare++/src/dynare_atoms.hh
@@ -37,9 +37,9 @@ namespace ogdyn
using std::vector;
using std::string;
- /** A definition of a type mapping a string to an integer. Used as
- * a substitution map, saying what names are substituted for what
- * expressions represented by tree indices. */
+ /* A definition of a type mapping a string to an integer. Used as a
+ substitution map, saying what names are substituted for what expressions
+ represented by tree indices. */
using Tsubstmap = map<string, int>;
class DynareStaticAtoms : public ogp::StaticAtoms
@@ -51,14 +51,13 @@ namespace ogdyn
}
DynareStaticAtoms(const DynareStaticAtoms &a) = default;
~DynareStaticAtoms() override = default;
- /** This registers a unique varname identifier. It throws an
- * exception if the variable name is duplicate. It checks the
- * uniqueness and then it calls StaticAtoms::register_name. */
+ /* This registers a unique varname identifier. It throws an exception if
+ the variable name is duplicate. It checks the uniqueness and then it
+ calls StaticAtoms::register_name. */
void register_name(string name) override;
protected:
- /** This returns a tree index of the given variable, and if
- * the variable has not been registered, it throws an
- * exception. */
+ /* This returns a tree index of the given variable, and if the variable has
+ not been registered, it throws an exception. */
int check_variable(const string &name) const override;
};
@@ -68,52 +67,50 @@ namespace ogdyn
enum class atype {endovar, exovar, param};
protected:
using Tatypemap = map<string, atype>;
- /** The map assigining a type to each name. */
+ /* The map assigining a type to each name. */
Tatypemap atom_type;
public:
DynareDynamicAtoms()
: ogp::SAtoms()
{
}
- /** This parses a variable of the forms: varname(+3),
- * varname(3), varname, varname(-3), varname(0), varname(+0),
- * varname(-0). */
+ /* This parses a variable of the forms: varname(+3), varname(3), varname,
+ varname(-3), varname(0), varname(+0), varname(-0). */
void parse_variable(const string &in, std::string &out, int &ll) const override;
- /** Registers unique name of endogenous variable. */
+ /* Registers unique name of endogenous variable. */
void register_uniq_endo(string name) override;
- /** Registers unique name of exogenous variable. */
+ /* Registers unique name of exogenous variable. */
void register_uniq_exo(string name) override;
- /** Registers unique name of parameter. */
+ /* Registers unique name of parameter. */
void register_uniq_param(string name) override;
- /** Return true if the name is a given type. */
+ /* Return true if the name is a given type. */
bool is_type(const string &name, atype tp) const;
- /** Debug print. */
+ /* Debug print. */
void print() const override;
- /** Implement NularyStringConvertor::convert. */
+ /* Implement NularyStringConvertor::convert. */
std::string convert(int t) const override;
};
- /** This class represents the atom values for dynare, where
- * exogenous variables can occur only at time t, and endogenous at
- * times t-1, t, and t+1. */
+ /* This class represents the atom values for dynare, where exogenous
+ variables can occur only at time t, and endogenous at times t−1, t, and
+ t+1. */
class DynareAtomValues : public ogp::AtomValues
{
protected:
- /** Reference to the atoms (we suppose that they are only at
- * t-1,t,t+1. */
+ /* Reference to the atoms (we suppose that they are only at t−1,t,t+1. */
const ogp::FineAtoms &atoms;
- /** De facto reference to the values of parameters. */
+ /* De facto reference to the values of parameters. */
const ConstVector paramvals;
- /** De facto reference to the values of endogenous at time t-1. Only
- * predetermined and both part. */
+ /* De facto reference to the values of endogenous at time t−1. Only
+ predetermined and both part. */
const ConstVector yym;
- /** De facto reference to the values of endogenous at time t. Ordering
- * given by the atoms. */
+ /* De facto reference to the values of endogenous at time t. Ordering given
+ by the atoms. */
const ConstVector yy;
- /** De facto reference to the values of endogenous at time t+1. Only
- * both and forward looking part. */
+ /* De facto reference to the values of endogenous at time t+1. Only both
+ and forward looking part. */
const ConstVector yyp;
- /** De facto reference to the values of exogenous at time t. */
+ /* De facto reference to the values of exogenous at time t. */
const ConstVector xx;
public:
DynareAtomValues(const ogp::FineAtoms &a, const Vector &pvals, const Vector &ym,
@@ -129,20 +126,20 @@ namespace ogdyn
void setValues(ogp::EvalTree &et) const override;
};
- /** This class represents the atom values at the steady state. It
- * makes only appropriate subvector yym and yyp of the y vector,
- * makes a vector of zero exogenous variables and uses
- * DynareAtomValues with more general interface. */
+ /* This class represents the atom values at the steady state. It makes only
+ appropriate subvector yym and yyp of the y vector, makes a vector of zero
+ exogenous variables and uses DynareAtomValues with more general
+ interface. */
class DynareSteadyAtomValues : public ogp::AtomValues
{
protected:
- /** Subvector of yy. */
+ /* Subvector of yy. */
const ConstVector yym;
- /** Subvector of yy. */
+ /* Subvector of yy. */
const ConstVector yyp;
- /** Vector of zeros for exogenous variables. */
+ /* Vector of zeros for exogenous variables. */
Vector xx;
- /** Atom values using this yym, yyp and xx. */
+ /* Atom values using this yym, yyp and xx. */
DynareAtomValues av;
public:
DynareSteadyAtomValues(const ogp::FineAtoms &a, const Vector &pvals, const Vector &y)
@@ -163,21 +160,19 @@ namespace ogdyn
class DynareStaticSteadyAtomValues : public ogp::AtomValues
{
protected:
- /** Reference to static atoms over which the tree, where the
- * values go, is defined. */
+ /* Reference to static atoms over which the tree, where the values go, is
+ defined. */
const ogp::StaticFineAtoms &atoms_static;
- /** Reference to dynamic atoms for which the class gets input
- * data. */
+ /* Reference to dynamic atoms for which the class gets input data. */
const ogp::FineAtoms &atoms;
- /** De facto reference to input data, this is a vector of
- * endogenous variables in internal ordering of the dynamic
- * atoms. */
+ /* De facto reference to input data, this is a vector of endogenous
+ variables in internal ordering of the dynamic atoms. */
ConstVector yy;
- /** De facto reference to input parameters corresponding to
- * ordering defined by the dynamic atoms. */
+ /* De facto reference to input parameters corresponding to ordering defined
+ by the dynamic atoms. */
ConstVector paramvals;
public:
- /** Construct the object. */
+ /* Construct the object. */
DynareStaticSteadyAtomValues(const ogp::FineAtoms &a, const ogp::StaticFineAtoms &sa,
const Vector &pvals, const Vector &yyy)
: atoms_static(sa),
@@ -186,18 +181,17 @@ namespace ogdyn
paramvals(pvals)
{
}
- /** Set the values to the tree defined over the static atoms. */
+ /* Set the values to the tree defined over the static atoms. */
void setValues(ogp::EvalTree &et) const override;
};
- /** This class takes a vector of endogenous variables and a
- * substitution map. It supposes that variables at the right hand
- * sides of the substitutions are set in the endogenous vector. It
- * evaluates the substitutions and if the variables corresponding
- * to left hand sides are not set in the endogenous vector it sets
- * them to calculated values. If a variable is already set, it
- * does not override its value. It has no methods, everything is
- * done in the constructor. */
+ /* This class takes a vector of endogenous variables and a substitution map.
+ It supposes that variables at the right hand sides of the substitutions
+ are set in the endogenous vector. It evaluates the substitutions and if
+ the variables corresponding to left hand sides are not set in the
+ endogenous vector it sets them to calculated values. If a variable is
+ already set, it does not override its value. It has no methods, everything
+ is done in the constructor. */
class DynareSteadySubstitutions : public ogp::FormulaEvalLoader
{
protected:
@@ -213,11 +207,10 @@ namespace ogdyn
vector<int> right_hand_sides;
};
- /** This class is a static version of DynareSteadySustitutions. It
- * works for static atoms and static tree and substitution map
- * over the static tree. It also needs dynamic version of the
- * atoms, since it defines ordering of the vectors pvals, and
- * yy. */
+ /* This class is a static version of DynareSteadySustitutions. It works for
+ static atoms and static tree and substitution map over the static tree. It
+ also needs dynamic version of the atoms, since it defines ordering of the
+ vectors pvals, and yy. */
class DynareStaticSteadySubstitutions : public ogp::FormulaEvalLoader
{
protected:
diff --git a/dynare++/src/dynare_model.cc b/dynare++/src/dynare_model.cc
index fbf20247cd8ec02c530bc7affb91c53173cf6f81..d7d8abb6a2ce82c05457944e248fbdac45d898cb 100644
--- a/dynare++/src/dynare_model.cc
+++ b/dynare++/src/dynare_model.cc
@@ -192,8 +192,8 @@ DynareModel::check_model() const
throw DynareException(__FILE__, __LINE__, "Model has " + std::to_string(eqs.nformulas())
+ " equations for " + std::to_string(atoms.ny()) + " endogenous variables");
- // check whether all nulary terms of all formulas in eqs are
- // either constant or assigned to a name
+ /* check whether all nulary terms of all formulas in eqs are either constant
+ or assigned to a name */
for (int i = 0; i < eqs.nformulas(); i++)
{
int ft = eqs.formula(i);
@@ -295,12 +295,11 @@ DynareModel::get_nonlinear_subterms(int t) const
void
DynareModel::substitute_atom_for_term(const string &name, int ll, int t)
{
- // if the term t is itself a named atom (parameter, exo, endo),
- // then we have to unassign it first
+ /* if the term t is itself a named atom (parameter, exo, endo), then we have
+ to unassign it first */
if (atoms.is_named_atom(t))
atoms.unassign_variable(atoms.name(t), atoms.lead(t), t);
- // assign allocated tree index
- // for the term now to name(ll)
+ /* assign allocated tree index for the term now to name(ll) */
atoms.assign_variable(name, ll, t);
// make operation t nulary in operation tree
eqs.nularify(t);
@@ -311,8 +310,8 @@ DynareModel::final_job()
{
if (t_plobjective != -1 && t_pldiscount != -1)
{
- // at this moment include all equations and all variables; in
- // future we will exclude purely exogenous processes; todo:
+ /* at this moment include all equations and all variables; in future we
+ will exclude purely exogenous processes; TODO */
PlannerBuilder::Tvarset vset;
for (int i = 0; i < atoms.ny(); i++)
vset.insert(atoms.get_endovars()[i]);
@@ -320,8 +319,7 @@ DynareModel::final_job()
for (int i = 0; i < eqs.nformulas(); i++)
eset.push_back(i);
- // construct the planner builder, this adds a lot of stuff to
- // the model
+ // construct the planner builder, this adds a lot of stuff to the model
pbuilder = std::make_unique<PlannerBuilder>(*this, vset, eset);
}
@@ -334,8 +332,8 @@ DynareModel::final_job()
old_atoms = std::make_unique<DynareDynamicAtoms>(atoms);
// construct empty substitutions from old_atoms to atoms
atom_substs = std::make_unique<ogp::AtomSubstitutions>(*old_atoms, atoms);
- // do the actual substitution, it will also call
- // parsing_finished for atoms which creates internal orderings
+ /* do the actual substitution, it will also call parsing_finished for atoms
+ which creates internal orderings */
atoms.substituteAllLagsAndExo1Leads(eqs, *atom_substs);
}
@@ -505,8 +503,7 @@ DynareParser::add_name(string name, int flag)
void
DynareParser::error(string mes)
{
- // throwing zero offset since this exception will be caugth at
- // constructor
+ // throwing zero offset since this exception will be caugth at constructor
throw ogp::ParserException(std::move(mes), 0);
}
@@ -523,12 +520,11 @@ DynareParser::print() const
<< "initval position: " << initval_beg << ' ' << initval_end << '\n';
}
-/** A global symbol for passing info to the DynareParser from
- * parser. */
+/* A global symbol for passing info to the DynareParser from parser. */
DynareParser *dynare_parser;
-/** The declarations of functions defined in dynglob_ll.cc and
- * dynglob_tab.cc generated from dynglob.lex and dynglob.y */
+/* The declarations of functions defined in dynglob_ll.cc and dynglob_tab.cc
+ generated from dynglob.lex and dynglob.y */
void *dynglob__scan_string(const char *);
void dynglob__destroy_buffer(void *);
void dynglob_parse();
@@ -601,13 +597,13 @@ DynareParser::calc_init()
(*init_vals)[i] = aae.get_value(atoms.get_endovars()[outer]);
}
- // if the planner's FOCs have been added, then add estimate of
- // Lagrange multipliers to the vector
+ /* if the planner's FOCs have been added, then add estimate of Lagrange
+ multipliers to the vector */
if (pbuilder)
MultInitSS mis(*pbuilder, *param_vals, *init_vals);
- // if forward substitution builder has been created, we have to
- // its substitutions and evaluate them
+ /* if forward substitution builder has been created, we have to its
+ substitutions and evaluate them */
if (fbuilder)
ogdyn::DynareSteadySubstitutions dss(atoms, eqs.getTree(),
fbuilder->get_aux_map(), *param_vals, *init_vals);
@@ -664,8 +660,8 @@ NLSelector::operator()(int t) const
return false;
else
return true;
- // DIVIDE is linear if the denominator is constant, or if
- // the nominator is zero
+ /* DIVIDE is linear if the denominator is constant, or if the nominator
+ is zero */
if (op.getCode() == ogp::code_t::DIVIDE
&& (op.getOp1() == ogp::OperationTree::zero
|| model.is_constant_term(op.getOp2())))
diff --git a/dynare++/src/dynare_model.hh b/dynare++/src/dynare_model.hh
index ea977df4e4be7cbe97bbf21c18acae1a4bfc2de8..18c428df009a69f66ef5c2b0fbe414bae79a93f0 100644
--- a/dynare++/src/dynare_model.hh
+++ b/dynare++/src/dynare_model.hh
@@ -43,10 +43,9 @@ namespace ogdyn
using std::unordered_set;
using std::map;
- /** This represents an interval in a string by the pair of
- * positions (including the first, excluding the second). A
- * position is given by the line and the column within the line
- * (both starting from 1). */
+ /* This represents an interval in a string by the pair of positions
+ (including the first, excluding the second). A position is given by the
+ line and the column within the line (both starting from 1). */
struct PosInterval
{
int fl;
@@ -59,7 +58,7 @@ namespace ogdyn
{
}
PosInterval &operator=(const PosInterval &pi) = default;
- /** Debug print. */
+ /* Debug print. */
void
print() const
{
@@ -67,12 +66,12 @@ namespace ogdyn
}
};
- /** This class is basically a GeneralMatrix but is created from
- * parsed matrix data. */
+ /* This class is basically a GeneralMatrix but is created from parsed matrix
+ data. */
class ParsedMatrix : public TwoDMatrix
{
public:
- /** Construct the object from the parsed data of ogp::MatrixParser. */
+ /* Construct the object from the parsed data of ogp::MatrixParser. */
ParsedMatrix(const ogp::MatrixParser &mp);
};
@@ -83,8 +82,8 @@ namespace ogdyn
class MultInitSS;
class ModelSSWriter;
- /** A subclass is responsible for creating param_vals, init_vals,
- * and vcov_mat. */
+ /* A subclass is responsible for creating param_vals, init_vals, and
+ vcov_mat. */
class DynareModel
{
friend class PlannerBuilder;
@@ -92,46 +91,43 @@ namespace ogdyn
friend class MultInitSS;
friend class ModelSSWriter;
protected:
- /** All atoms for whole model. */
+ /* All atoms for whole model. */
DynareDynamicAtoms atoms;
- /** Parsed model equations. */
+ /* Parsed model equations. */
ogp::FormulaParser eqs;
- /** Order of approximation. */
+ /* Order of approximation. */
int order{-1};
- /** A vector of parameters values created by a subclass. It
- * is stored with natural ordering (outer) of the parameters
- * given by atoms. */
+ /* A vector of parameters values created by a subclass. It is stored with
+ natural ordering (outer) of the parameters given by atoms. */
std::unique_ptr<Vector> param_vals;
- /** A vector of initial values created by a subclass. It is
- * stored with internal ordering given by atoms. */
+ /* A vector of initial values created by a subclass. It is stored with
+ internal ordering given by atoms. */
std::unique_ptr<Vector> init_vals;
- /** A matrix for vcov. It is created by a subclass. */
+ /* A matrix for vcov. It is created by a subclass. */
std::unique_ptr<TwoDMatrix> vcov_mat;
- /** Tree index of the planner objective. If there was no
- * planner objective keyword, the value is set to -1. */
+ /* Tree index of the planner objective. If there was no planner objective
+ keyword, the value is set to −1. */
int t_plobjective{-1};
- /** Tree index of the planner discount. If there was no
- * planner discount keyword, the value is set to -1. */
+ /* Tree index of the planner discount. If there was no planner discount
+ keyword, the value is set to −1. */
int t_pldiscount{-1};
- /** Pointer to PlannerBuilder, which is created only if the
- * planner's FOC are added to the model. */
+ /* Pointer to PlannerBuilder, which is created only if the planner's FOC
+ are added to the model. */
std::unique_ptr<PlannerBuilder> pbuilder;
- /** Pointer to an object which builds auxiliary variables and
- * equations to rewrite a model containing multiple leads to
- * an equivalent model having only +1 leads. */
+ /* Pointer to an object which builds auxiliary variables and equations to
+ rewrite a model containing multiple leads to an equivalent model having
+ only +1 leads. */
std::unique_ptr<ForwSubstBuilder> fbuilder;
- /** Pointer to AtomSubstitutions which are created when the
- * atoms are being substituted because of multiple lags
- * etc. It uses also an old copy of atoms, which is
- * created. */
+ /* Pointer to AtomSubstitutions which are created when the atoms are being
+ substituted because of multiple lags etc. It uses also an old copy of
+ atoms, which is created. */
std::unique_ptr<ogp::AtomSubstitutions> atom_substs;
- /** Pointer to a copy of original atoms before substitutions
- * took place. */
+ /* Pointer to a copy of original atoms before substitutions took place. */
std::unique_ptr<ogp::SAtoms> old_atoms;
public:
- /** Initializes the object to an empty state. */
+ /* Initializes the object to an empty state. */
DynareModel();
- /** Construct a new deep copy. */
+ /* Construct a new deep copy. */
DynareModel(const DynareModel &dm);
virtual ~DynareModel() = default;
virtual std::unique_ptr<DynareModel> clone() const = 0;
@@ -150,7 +146,7 @@ namespace ogdyn
{
return order;
}
- /** Return the vector of parameter values. */
+ /* Return the vector of parameter values. */
const Vector &
getParams() const
{
@@ -161,7 +157,7 @@ namespace ogdyn
{
return *param_vals;
}
- /** Return the vector of initial values of endo variables. */
+ /* Return the vector of initial values of endo variables. */
const Vector &
getInit() const
{
@@ -172,7 +168,7 @@ namespace ogdyn
{
return *init_vals;
}
- /** Return the vcov matrix. */
+ /* Return the vcov matrix. */
const TwoDMatrix &
getVcov() const
{
@@ -183,94 +179,85 @@ namespace ogdyn
{
return *vcov_mat;
}
- /** Return planner info. */
+ /* Return planner info. */
const PlannerInfo *get_planner_info() const;
- /** Return forward substitutions info. */
+ /* Return forward substitutions info. */
const ForwSubstInfo *get_forw_subst_info() const;
- /** Return substitutions info. */
+ /* Return substitutions info. */
const ogp::SubstInfo *get_subst_info() const;
- /** This sets initial values given in outer ordering. */
+ /* This sets initial values given in outer ordering. */
void setInitOuter(const Vector &x);
- /** This returns true if the given term is a function of
- * hardwired constants, numerical constants and parameters. */
+ /* This returns true if the given term is a function of hardwired
+ constants, numerical constants and parameters. */
bool is_constant_term(int t) const;
- /** Debug print. */
+ /* Debug print. */
void print() const;
- /** Dump the model to the output stream. This includes
- * variable declarations, parameter values, model code,
- * initval, vcov and order. */
+ /* Dump the model to the output stream. This includes variable
+ declarations, parameter values, model code, initval, vcov and order. */
void dump_model(std::ostream &os) const;
protected:
- /** Adds a name of endogenous, exogenous or a parameter. The
- * sort is governed by the flag. See dynglob.y for values of
- * the flag. This is used by a subclass when declaring the
- * names. */
+ /* Adds a name of endogenous, exogenous or a parameter. The sort is
+ governed by the flag. See dynglob.yy for values of the flag. This is
+ used by a subclass when declaring the names. */
void add_name(std::string name, int flag);
- /** This checks the model consistency. Thus includes: number
- * of endo variables and number of equations, min and max lag
- * of endogenous variables and occurrrences of exogenous
- * variables. It throws an exception, if there is a problem. */
+ /* This checks the model consistency. Thus includes: number of endo
+ variables and number of equations, min and max lag of endogenous
+ variables and occurrrences of exogenous variables. It throws an
+ exception, if there is a problem. */
void check_model() const;
- /** This shifts the given variable identified by the tree
- * index in time. So if the given tree index represents a(+3)
- * and the tshift is -4, the method returns tree index of the
- * a(-1). If a(-1) doesn't exist, it is added to the tree. If
- * it exists, its tree index is returned. If the tree index
- * doesn't correspond to an endogenous nor exogenous variable,
- * an exception is thrown. */
+ /* This shifts the given variable identified by the tree index in time. So
+ if the given tree index represents a(+3) and the tshift is −4, the
+ method returns tree index of the a(-1). If a(-1) doesn't exist, it is
+ added to the tree. If it exists, its tree index is returned. If the tree
+ index doesn't correspond to an endogenous nor exogenous variable, an
+ exception is thrown. */
int variable_shift(int t, int tshift);
- /** For the given set of atoms identified by tree indices and
- * given time shift, this method returns a map mapping each
- * variable in the given set to its time shifted variable. The
- * map is passed through the reference and is cleared in the
- * beginning. */
+ /* For the given set of atoms identified by tree indices and given time
+ shift, this method returns a map mapping each variable in the given set
+ to its time shifted variable. The map is passed through the reference
+ and is cleared in the beginning. */
void variable_shift_map(const unordered_set<int> &a_set, int tshift,
map<int, int> &s_map);
- /** This returns maximum lead and minimum lag of an endogenous
- * or exogenous variable in the given term. If there are no
- * endo or exo variables, than it returns the least integer as
- * max lead and the greatest integer as min lag. */
+ /* This returns maximum lead and minimum lag of an endogenous or exogenous
+ variable in the given term. If there are no endo or exo variables, than
+ it returns the least integer as max lead and the greatest integer as min
+ lag. */
void termspan(int t, int &mlead, int &mlag) const;
- /** This function returns a set of non-linear subterms of the
- * given term, these are terms whose linear combination
- * constitutes the given term. */
+ /* This function returns a set of non-linear subterms of the given term,
+ these are terms whose linear combination constitutes the given term. */
unordered_set<int> get_nonlinear_subterms(int t) const;
- /** This method assigns already used tree index of some term
- * to the not-yet used atom name with the given lead/lag. In
- * this way, all occurrences of term t are substituted with
- * the atom name(ll). The method handles also rewriting
- * operation tree including derivatives of the term t. */
+ /* This method assigns already used tree index of some term to the not-yet
+ used atom name with the given lead/lag. In this way, all occurrences of
+ term t are substituted with the atom name(ll). The method handles also
+ rewriting operation tree including derivatives of the term t. */
void substitute_atom_for_term(const string &name, int ll, int t);
- /** This performs a final job after the model is parsed. It
- * creates the PlannerBuilder object if the planner's FOC are
- * needed, then it creates ForwSubstBuilder handling multiple
- * leads and finally it creates the substitution object saving
- * old atoms and performs the substitutions. */
+ /* This performs a final job after the model is parsed. It creates the
+ PlannerBuilder object if the planner's FOC are needed, then it creates
+ ForwSubstBuilder handling multiple leads and finally it creates the
+ substitution object saving old atoms and performs the substitutions. */
void final_job();
};
- /** This class constructs DynareModel from dynare++ model file. It
- * parses variable declarations, model equations, parameter
- * assignments, initval assignments, vcov matrix and order of
- * approximation. */
+ /* This class constructs DynareModel from dynare++ model file. It parses
+ variable declarations, model equations, parameter assignments, initval
+ assignments, vcov matrix and order of approximation. */
class DynareParser : public DynareModel
{
protected:
- /** Static atoms for parameter assignments. */
+ /* Static atoms for parameter assignments. */
DynareStaticAtoms pa_atoms;
- /** Assignments for the parameters. */
+ /* Assignments for the parameters. */
ogp::AtomAssignings paramset;
- /** Static atoms for initval assignments. */
+ /* Static atoms for initval assignments. */
DynareStaticAtoms ia_atoms;
- /** Assignments for the initval. */
+ /* Assignments for the initval. */
ogp::AtomAssignings initval;
- /** Matrix parser for vcov. */
+ /* Matrix parser for vcov. */
ogp::MatrixParser vcov;
public:
- /** This, in fact, creates DynareModel from the given string
- * of the given length corresponding to the Dynare++ model
- * file. If the given ord is not -1, then it overrides setting
- * in the model file. */
+ /* This, in fact, creates DynareModel from the given string of the given
+ length corresponding to the Dynare++ model file. If the given ord is not
+ −1, then it overrides setting in the model file. */
DynareParser(const string &str, int ord);
DynareParser(const DynareParser &dp);
std::unique_ptr<DynareModel>
@@ -278,85 +265,80 @@ namespace ogdyn
{
return std::make_unique<DynareParser>(*this);
}
- /** Adds a name of endogenous, exogenous or a parameter. This
- * addss the name to the parent class DynareModel and also
- * registers the name to either paramset, or initval. */
+ /* Adds a name of endogenous, exogenous or a parameter. This addss the name
+ to the parent class DynareModel and also registers the name to either
+ paramset, or initval. */
void add_name(string name, int flag);
- /** Sets position of the model section. Called from
- * dynglob.y. */
+ /* Sets position of the model section. Called from dynglob.yy. */
void
set_model_pos(int off1, int off2)
{
model_beg = off1;
model_end = off2;
}
- /** Sets position of the section setting parameters. Called
- * from dynglob.y. */
+ /* Sets position of the section setting parameters. Called from
+ dynglob.yy. */
void
set_paramset_pos(int off1, int off2)
{
paramset_beg = off1;
paramset_end = off2;
}
- /** Sets position of the initval section. Called from
- * dynglob.y. */
+ /* Sets position of the initval section. Called from dynglob.yy. */
void
set_initval_pos(int off1, int off2)
{
initval_beg = off1;
initval_end = off2;
}
- /** Sets position of the vcov section. Called from
- * dynglob.y. */
+ /* Sets position of the vcov section. Called from dynglob.yy. */
void
set_vcov_pos(int off1, int off2)
{
vcov_beg = off1;
vcov_end = off2;
}
- /** Parser the given string as integer and set to as the
- * order. */
+ /* Parser the given string as integer and set to as the order. */
void
set_order_pos(int off1, int off2)
{
order_beg = off1;
order_end = off2;
}
- /** Sets position of the planner_objective section. Called
- * from dynglob.y. */
+ /* Sets position of the planner_objective section. Called from
+ dynglob.yy. */
void
set_pl_objective_pos(int off1, int off2)
{
plobjective_beg = off1;
plobjective_end = off2;
}
- /** Sets position of the planner_discount section. Called from
- * dynglob.y. */
+ /* Sets position of the planner_discount section. Called from
+ dynglob.yy. */
void
set_pl_discount_pos(int off1, int off2)
{
pldiscount_beg = off1;
pldiscount_end = off2;
}
- /** Processes a syntax error from bison. */
+ /* Processes a syntax error from bison. */
void error(string mes);
- /** Debug print. */
+ /* Debug print. */
void print() const;
protected:
void parse_glob(const string &stream);
int parse_order(const string &stream);
int parse_pldiscount(const string &stream);
- /** Evaluate paramset assignings and set param_vals. */
+ /* Evaluate paramset assignings and set param_vals. */
void calc_params();
- /** Evaluate initval assignings and set init_vals. */
+ /* Evaluate initval assignings and set init_vals. */
void calc_init();
- /** Do the final job. This includes building the planner
- * problem (if any) and substituting for multiple lags, and
- * one period leads of exogenous variables, and calculating
- * initial guess of lagrange multipliers in the social planner
- * problem. Precondtion: everything parsed and calculated
- * parameters, postcondition: calculated initvals vector and
- * parsing_finished for expanded vectors. */
+ /* Do the final job. This includes building the planner problem (if any)
+ and substituting for multiple lags, and one period leads of exogenous
+ variables, and calculating initial guess of lagrange multipliers in the
+ social planner problem. Precondtion: everything parsed and calculated
+ parameters, postcondition: calculated initvals vector and
+ parsing_finished for expanded vectors. */
void final_job();
private:
int model_beg, model_end;
@@ -368,12 +350,11 @@ namespace ogdyn
int pldiscount_beg, pldiscount_end;
};
- /** Semiparsed model. The equations are given by a string,
- * everything other by C/C++ objects. The initial values are set
- * manually after the creation of this object. This implies that
- * no automatic substitutions cannot be done here, which in turn
- * implies that we cannot do here a social planner nor substitutions
- * of multiple lags. */
+ /* Semiparsed model. The equations are given by a string, everything other by
+ C++ objects. The initial values are set manually after the creation of
+ this object. This implies that no automatic substitutions cannot be done
+ here, which in turn implies that we cannot do here a social planner nor
+ substitutions of multiple lags. */
class DynareSPModel : public DynareModel
{
public:
@@ -390,10 +371,9 @@ namespace ogdyn
}
};
- /** This class implements a selector of operations which correspond
- * to non-linear functions. This inherits from ogp::opselector and
- * is used to calculate non-linear subterms in
- * DynareModel::get_nonlinear_subterms(). */
+ /* This class implements a selector of operations which correspond to
+ non-linear functions. This inherits from ogp::opselector and is used to
+ calculate non-linear subterms in DynareModel::get_nonlinear_subterms(). */
class NLSelector : public ogp::opselector
{
private:
@@ -405,10 +385,9 @@ namespace ogdyn
bool operator()(int t) const override;
};
- /** This class writes a mathematical code evaluating the system of
- * equations and the first derivatives at zero shocks and at the
- * given (static) state. Static means that lags and leads are
- * ignored. */
+ /* This class writes a mathematical code evaluating the system of equations
+ and the first derivatives at zero shocks and at the given (static) state.
+ Static means that lags and leads are ignored. */
class ModelSSWriter : public ogp::DefaultOperationFormatter
{
protected:
@@ -419,15 +398,14 @@ namespace ogdyn
model(m)
{
}
- /** This writes the evaluation of the system. It calls pure
- * virtual methods for writing a preamble, then assignment of
- * atoms, and then assignment for resulting object. These are
- * language dependent and are implemented in the subclass. */
+ /* This writes the evaluation of the system. It calls pure virtual methods
+ for writing a preamble, then assignment of atoms, and then assignment
+ for resulting object. These are language dependent and are implemented
+ in the subclass. */
void write_der0(std::ostream &os);
- /** This writes the evaluation of the first order derivative of
- the system. It calls pure virtual methods for writing a
- preamble, assignment, and assignemnt of the resulting
- objects. */
+ /* This writes the evaluation of the first order derivative of the system.
+ It calls pure virtual methods for writing a preamble, assignment, and
+ assignemnt of the resulting objects. */
void write_der1(std::ostream &os);
protected:
virtual void write_der0_preamble(std::ostream &os) const = 0;
@@ -440,7 +418,7 @@ namespace ogdyn
class MatlabSSWriter : public ModelSSWriter
{
protected:
- /** Identifier used in function names. */
+ /* Identifier used in function names. */
std::string id;
public:
MatlabSSWriter(const DynareModel &dm, std::string id_arg);
@@ -448,27 +426,26 @@ namespace ogdyn
// from ModelSSWriter
void write_der0_preamble(std::ostream &os) const override;
void write_der1_preamble(std::ostream &os) const override;
- /** This writes atom assignments. We have four kinds of atoms
- * set here: endogenous vars coming from one parameter,
- * parameter values given by the second parameter, constants,
- * and the OperationTree::num_constants hardwired constants in
- * ogp::OperationTree. */
+ /* This writes atom assignments. We have four kinds of atoms set here:
+ endogenous vars coming from one parameter, parameter values given by the
+ second parameter, constants, and the OperationTree::num_constants
+ hardwired constants in ogp::OperationTree. */
void write_atom_assignment(std::ostream &os) const override;
void write_der0_assignment(std::ostream &os) const override;
void write_der1_assignment(std::ostream &os) const override;
- /** This prints t10 for t=10. */
+ /* This prints t10 for t=10. */
void format_term(int t, std::ostream &os) const override;
- /** This prints a10 for t=10. The atoms a10 are supposed to be
- * set by write_atom_assignments(). */
+ /* This prints a10 for t=10. The atoms a10 are supposed to be set by
+ write_atom_assignments(). */
void format_nulary(int t, std::ostream &os) const override;
private:
void write_common1_preamble(std::ostream &os) const;
void write_common2_preamble(std::ostream &os) const;
};
- /** This class implements OperationFormatter for debugging
- * purposes. It renders atoms in a more friendly way than the
- * ogp::DefaulOperationFormatter. */
+ /* This class implements OperationFormatter for debugging purposes. It
+ renders atoms in a more friendly way than the
+ ogp::DefaulOperationFormatter. */
class DebugOperationFormatter : public ogp::DefaultOperationFormatter
{
protected:
diff --git a/dynare++/src/dynare_params.hh b/dynare++/src/dynare_params.hh
index ca8318cbb4edda08a0ccf6b2936de6ca846d23c2..79e8a1ea0b7d0f99bb0c81bcbc38fdf4874a5c1d 100644
--- a/dynare++/src/dynare_params.hh
+++ b/dynare++/src/dynare_params.hh
@@ -20,8 +20,8 @@
/*
along shocks: m mult max_evals
- ellipse: m mult max_evals (10*m) (0.5*mult)
- simul: m max_evals (10*m)
+ ellipse: m mult max_evals (10·m) (0.5·mult)
+ simul: m max_evals (10·m)
--check-scale 2.0 --check-evals 1000 --check-num 10 --check PES
*/
@@ -45,7 +45,7 @@ struct DynareParams
std::string prefix;
int seed;
int order;
- /** Tolerance used for steady state calcs. */
+ /* Tolerance used for steady state calcs. */
double ss_tol;
bool check_along_path;
bool check_along_shocks;
@@ -53,9 +53,9 @@ struct DynareParams
int check_evals;
int check_num;
double check_scale;
- /** Flag for doing IRFs even if the irf_list is empty. */
+ /* Flag for doing IRFs even if the irf_list is empty. */
bool do_irfs_all;
- /** List of shocks for which IRF will be calculated. */
+ /* List of shocks for which IRF will be calculated. */
std::vector<std::string> irf_list;
bool do_centralize;
double qz_criterium;
@@ -95,8 +95,8 @@ private:
check_evals, check_scale, check_num, noirfs, irfs,
help, version, centralize, no_centralize, qz_criterium};
void processCheckFlags(const std::string &flags);
- /** This gathers strings from argv[optind] and on not starting
- * with '-' to the irf_list. It stops one item before the end,
- * since this is the model file. */
+ /* This gathers strings from argv[optind] and on not starting with '-' to the
+ irf_list. It stops one item before the end, since this is the model
+ file. */
void processIRFList(int argc, char **argv);
};
diff --git a/dynare++/src/forw_subst_builder.cc b/dynare++/src/forw_subst_builder.cc
index e062ec652263d98d5e77c905a51a1a6ba072a230..e4e4f4cf25677db75dac1a98aed9c67bb1a1a043 100644
--- a/dynare++/src/forw_subst_builder.cc
+++ b/dynare++/src/forw_subst_builder.cc
@@ -51,8 +51,8 @@ ForwSubstBuilder::ForwSubstBuilder(DynareModel &m)
// unassign all variables with lead greater than 1
unassign_gt_1_leads();
- // forget the derivatives in the tree because some variables could
- // have been unassigned
+ /* Forget the derivatives in the tree because some variables could have been
+ unassigned */
model.eqs.getTree().forget_derivative_maps();
info.num_new_terms += model.getParser().getTree().get_num_op();
@@ -97,11 +97,9 @@ ForwSubstBuilder::substitute_for_term(int t, int i, int j)
auxt = model.eqs.add_nulary(name);
// add AUXLD_*_*_{ll+1} = AUXLD_*_*_{ll}(+1)
model.eqs.add_formula(model.eqs.add_binary(ogp::code_t::MINUS, auxt, lastauxt_lead));
- // add substitution to the map; todo: this
- // works well because in the context where
- // aux_map is used the timing doesn't matter,
- // however, it is misleading, needs to be
- // changed
+ /* add substitution to the map; TODO: this works well because in the
+ context where aux_map is used the timing doesn't matter, however,
+ it is misleading, needs to be changed */
aux_map.emplace(name, lagt);
}
diff --git a/dynare++/src/forw_subst_builder.hh b/dynare++/src/forw_subst_builder.hh
index 4083a9b6a9e2b1750882393ef8bab32c56cc632d..49f635b3cb66421d523d9a801b35b1a9f961ecc5 100644
--- a/dynare++/src/forw_subst_builder.hh
+++ b/dynare++/src/forw_subst_builder.hh
@@ -27,8 +27,8 @@
namespace ogdyn
{
- /** This struct encapsulates information about the process of
- * forward substitutions. */
+ /* This struct encapsulates information about the process of forward
+ substitutions. */
struct ForwSubstInfo
{
int num_affected_equations{0};
@@ -43,57 +43,53 @@ namespace ogdyn
{
using Ttermauxmap = map<int, string>;
protected:
- /** Reference to the model, to which we will add equations and
- * change some equations. */
+ /* Reference to the model, to which we will add equations and change some
+ equations. */
DynareModel &model;
- /** A map mapping new auxiliary variables to the terms in the
- * tree in the DynareModel. */
+ /* A map mapping new auxiliary variables to the terms in the tree in the
+ DynareModel. */
Tsubstmap aux_map;
- /** Information about the substitutions. */
+ /* Information about the substitutions. */
ForwSubstInfo info;
public:
- /** Do all the jobs needed. This scans all equations in the
- * model, and for equations containing forward looking
- * variables greater than 1 lead, it makes corresponding
- * substitutions. Basically, it breaks each equation to its
- * non-linear components and creates substitutions for these
- * components, not for whole equation. This is because the
- * expectation operator can go through the linear part of the
- * function. This will save us many occurrences of other
- * variables involved in the equation. */
+ /* Do all the jobs needed. This scans all equations in the model, and for
+ equations containing forward looking variables greater than 1 lead, it
+ makes corresponding substitutions. Basically, it breaks each equation to
+ its non-linear components and creates substitutions for these
+ components, not for whole equation. This is because the expectation
+ operator can go through the linear part of the function. This will save
+ us many occurrences of other variables involved in the equation. */
ForwSubstBuilder(DynareModel &m);
ForwSubstBuilder(const ForwSubstBuilder &b) = delete;
- /** Copy constructor with a new instance of the model. */
+ /* Copy constructor with a new instance of the model. */
ForwSubstBuilder(const ForwSubstBuilder &b, DynareModel &m);
- /** Return the auxiliary variable mapping. */
+ /* Return the auxiliary variable mapping. */
const Tsubstmap &
get_aux_map() const
{
return aux_map;
}
- /** Return the information. */
+ /* Return the information. */
const ForwSubstInfo &
get_info() const
{
return info;
}
private:
- /** This method takes a nonlinear term t, and if it has leads
- * of greater than 1, then it substitutes the term for the new
- * variable (or string of variables). Note that the
- * substitution is done by DynamicAtoms::assign_variable. This
- * means that the substitution is made for all other
- * ocurrences of t in the model. So there is no need of
- * tracking already substituted terms. The other two
- * parameters are just for identification of the new auxiliary
- * variables. When called from the constructor, i is an
- * equation number, j is an order of the non-linear term in
- * the equation. */
+ /* This method takes a nonlinear term t, and if it has leads of greater
+ than 1, then it substitutes the term for the new variable (or string of
+ variables). Note that the substitution is done by
+ DynamicAtoms::assign_variable. This means that the substitution is made
+ for all other ocurrences of t in the model. So there is no need of
+ tracking already substituted terms. The other two parameters are just
+ for identification of the new auxiliary variables. When called from the
+ constructor, i is an equation number, j is an order of the non-linear
+ term in the equation. */
void substitute_for_term(int t, int i, int j);
- /** This is called just at the end of the job. It unassigns
- * all nulary terms with a lead greater than 1. */
+ /* This is called just at the end of the job. It unassigns all nulary terms
+ with a lead greater than 1. */
void unassign_gt_1_leads();
- /** This unassigns all leads greater than 1 of the given name. */
+ /* This unassigns all leads greater than 1 of the given name. */
void unassign_gt_1_leads(const string &name);
};
};
diff --git a/dynare++/src/nlsolve.cc b/dynare++/src/nlsolve.cc
index a841522e5b3725aea455cd49001785e22e8e7648..566a9ad1c82d691969429a029464a00e7b260300 100644
--- a/dynare++/src/nlsolve.cc
+++ b/dynare++/src/nlsolve.cc
@@ -113,9 +113,9 @@ GoldenSectionSearch::init_bracket(OneDFunction &f, double x1, double &x2, double
// now we know that f1, f2, and fb are finite
if (std::isfinite(fbsym))
{
- // we have four numbers f1, fb, f2, fbsym, we test for the
- // following combinations to find the bracket:
- // [f1,f2,fbsym], [f1,fb,fbsym] and [f1,fb,fbsym]
+ /* we have four numbers f1, fb, f2, fbsym, we test for the following
+ combinations to find the bracket: [f1,f2,fbsym], [f1,fb,fbsym] and
+ [f1,fb,fbsym] */
if (f1 > f2 && f2 < fbsym)
{
bracket_found = true;
@@ -135,8 +135,8 @@ GoldenSectionSearch::init_bracket(OneDFunction &f, double x1, double &x2, double
// choose the smallest value in case we end
if (f1 > fbsym)
{
- // the smallest value is on the other end, we do
- // not want to continue
+ /* the smallest value is on the other end, we do not want to
+ continue */
b = bsym;
return false;
}
@@ -148,9 +148,9 @@ GoldenSectionSearch::init_bracket(OneDFunction &f, double x1, double &x2, double
}
else
{
- // we have only three numbers, we test for the bracket,
- // and if not found, we set b as potential result and
- // shorten x2 as potential init value for next cycle
+ /* we have only three numbers, we test for the bracket, and if not
+ found, we set b as potential result and shorten x2 as potential
+ init value for next cycle */
if (f1 > fb && fb < f2)
bracket_found = true;
else
@@ -172,8 +172,8 @@ GoldenSectionSearch::init_bracket(OneDFunction &f, double x1, double &x2, double
return bracket_found;
}
-/** This moves x2 toward to x1 until the function at x2 is finite and
- * b as a golden section between x1 and x2 yields also finite f. */
+/* This moves x2 toward to x1 until the function at x2 is finite and b as a
+ golden section between x1 and x2 yields also finite f. */
bool
GoldenSectionSearch::search_for_finite(OneDFunction &f, double x1, double &x2, double &b)
{
diff --git a/dynare++/src/nlsolve.hh b/dynare++/src/nlsolve.hh
index b7aeb9448130568db590c98d987bd754e4840cd3..b45ace078b7dd62a5f305c456678db6317e3b768 100644
--- a/dynare++/src/nlsolve.hh
+++ b/dynare++/src/nlsolve.hh
@@ -38,25 +38,22 @@ namespace ogu
class GoldenSectionSearch
{
protected:
- /** This should not be greater than DBL_EPSILON^(1/2). */
constexpr static double tol = 1.e-4;
- /** This is equal to the golden section ratio. */
+ /* This is equal to the golden section ratio. */
static double golden;
public:
static double search(OneDFunction &f, double x1, double x2);
protected:
- /** This initializes a bracket by moving x2 and b (as a golden
- * section of x1,x2) so that f(x1)>f(b) && f(b)<f(x2). The point
- * x1 is not moved, since it is considered as reliable and f(x1)
- * is supposed to be finite. If initialization of a bracket
- * succeeded, then [x1, b, x2] is the bracket and true is
- * returned. Otherwise, b is the minimum found and false is
- * returned. */
+ /* This initializes a bracket by moving x2 and b (as a golden section of
+ x1,x2) so that f(x1)>f(b) && f(b)<f(x2). The point x1 is not moved,
+ since it is considered as reliable and f(x1) is supposed to be finite.
+ If initialization of a bracket succeeded, then [x1,b,x2] is the bracket
+ and true is returned. Otherwise, b is the minimum found and false is
+ returned. */
static bool init_bracket(OneDFunction &f, double x1, double &x2, double &b);
- /** This supposes that f(x1) is finite and it moves x2 toward x1
- * until x2 and b (as a golden section of x1,x2) are finite. If
- * succeeded, the routine returns true and x2, and b. Otherwise,
- * it returns false. */
+ /* This supposes that f(x1) is finite and it moves x2 toward x1 until x2
+ and b (as a golden section of x1,x2) are finite. If succeeded, the
+ routine returns true and x2, and b. Otherwise, it returns false. */
static bool search_for_finite(OneDFunction &f, double x1, double &x2, double &b);
};
@@ -67,9 +64,9 @@ namespace ogu
virtual ~VectorFunction() = default;
virtual int inDim() const = 0;
virtual int outDim() const = 0;
- /** Check dimensions of eval parameters. */
+ /* Check dimensions of eval parameters. */
void check_for_eval(const ConstVector &in, Vector &out) const;
- /** Evaluate the vector function. */
+ /* Evaluate the vector function. */
virtual void eval(const ConstVector &in, Vector &out) = 0;
};
@@ -103,13 +100,12 @@ namespace ogu
xnewton.zeros(); xcauchy.zeros(); x.zeros();
}
~NLSolver() override = default;
- /** Returns true if the problem has converged. xx as input is the
- * starting value, as output it is a solution. */
+ /* Returns true if the problem has converged. xx as input is the starting
+ value, as output it is a solution. */
bool solve(Vector &xx, int &iter);
- /** To implement OneDFunction interface. It returns
- * func(xx)^T*func(xx), where
- * xx=x+lambda*xcauchy+(1-lambda)*xnewton. It is non-const only
- * because it calls func, x, xnewton, xcauchy is not changed. */
+ /* To implement OneDFunction interface. It returns func(xx)ᵀ·func(xx),
+ where xx=x+lambda·xcauchy+(1−lambda)·xnewton. It is non-const only
+ because it calls func, x, xnewton, xcauchy is not changed. */
double eval(double lambda) override;
};
};
diff --git a/dynare++/src/planner_builder.cc b/dynare++/src/planner_builder.cc
index 9f0d998d522074afdc59c2d56267fc86f5ee0520..90a15516ba2015fa02f000b7b31550408aa4f177 100644
--- a/dynare++/src/planner_builder.cc
+++ b/dynare++/src/planner_builder.cc
@@ -155,19 +155,18 @@ PlannerBuilder::shift_derivatives_of_f()
diff_f(yi, fi, ll-minlag)
= model.eqs.add_substitution(diff_f(yi, fi, ll-minlag), subst);
}
- // now do it for lags, these are put as leads under
- // expectations after time t, so we have to introduce
- // auxiliary variables at time t, and make leads of them here
+ /* now do it for lags, these are put as leads under expectations after
+ time t, so we have to introduce auxiliary variables at time t, and
+ make leads of them here */
for (int ll = minlag; ll < 0; ll++)
{
int ft = diff_f(yi, fi, ll-minlag);
if (ft != ogp::OperationTree::zero)
{
- // if the ft term has a lead, than we need to
- // introduce an auxiliary variable z_t, define it
- // as E_t[ft] and put z_{t-ll} to the
- // equation. Otherwise, we just put leaded ft to
- // the equation directly.
+ /* if the ft term has a lead, than we need to introduce an
+ auxiliary variable zₜ, define it as 𝔼ₜ[ft] and put z_{t-ll}
+ to the equation. Otherwise, we just put leaded ft to the
+ equation directly. */
int ft_maxlead, ft_minlag;
model.termspan(ft, ft_maxlead, ft_minlag);
if (ft_maxlead > 0)
@@ -187,11 +186,10 @@ PlannerBuilder::shift_derivatives_of_f()
}
else
{
- // no auxiliary variable is needed and the
- // term ft can be leaded in place
+ /* no auxiliary variable is needed and the term ft can be
+ leaded in place */
model.variable_shift_map(model.eqs.nulary_of_term(ft), -ll, subst);
- diff_f(yi, fi, ll-minlag)
- = model.eqs.add_substitution(ft, subst);
+ diff_f(yi, fi, ll-minlag) = model.eqs.add_substitution(ft, subst);
}
}
}
@@ -291,7 +289,7 @@ PlannerBuilder::lagrange_mult_f()
void
PlannerBuilder::form_equations()
{
- // add planner's FOCs
+ // add planner’s FOCs
for (int yi = 0; yi < diff_f.dim1(); yi++)
{
int eq = ogp::OperationTree::zero;
@@ -344,10 +342,9 @@ MultInitSS::MultInitSS(const PlannerBuilder &pb, const Vector &pvals, Vector &yy
builder.static_tree,
builder.static_aux_map, pvals, yy);
- // gather all the terms from builder.diff_b_static and
- // builder.diff_f_static to the vector, the ordering is important,
- // since the index of this vector will have to be decoded to the
- // position in b and F.
+ /* gather all the terms from builder.diff_b_static and builder.diff_f_static
+ to the vector, the ordering is important, since the index of this vector
+ will have to be decoded to the position in b and F. */
vector<int> terms;
for (int yi = 0; yi < builder.diff_b_static.nrows(); yi++)
for (int l = 0; l < builder.diff_b_static.ncols(); l++)
@@ -357,8 +354,8 @@ MultInitSS::MultInitSS(const PlannerBuilder &pb, const Vector &pvals, Vector &yy
for (int l = 0; l < builder.diff_f_static.dim3(); l++)
terms.push_back(builder.diff_f_static(yi, fi, l));
- // evaluate the terms, it will call a series of load(i,res), which
- // sum the results through lags/leads to b and F
+ /* evaluate the terms, it will call a series of load(i,res), which sum the
+ results through lags/leads to b and F */
DynareStaticSteadyAtomValues dssav(builder.model.atoms, builder.static_atoms, pvals, yy);
ogp::FormulaCustomEvaluator fe(builder.static_tree, terms);
fe.eval(dssav, *this);
@@ -378,8 +375,7 @@ MultInitSS::MultInitSS(const PlannerBuilder &pb, const Vector &pvals, Vector &yy
if (!std::isfinite(yy[iy]))
yy[iy] = lambda[fi];
- // go through all substitutions of the multiplier and set them
- // as well
+ /* go through all substitutions of the multiplier and set them as well */
if (builder.model.atom_substs)
{
const ogp::AtomSubstitutions::Toldnamemap &old2new
@@ -404,8 +400,8 @@ MultInitSS::MultInitSS(const PlannerBuilder &pb, const Vector &pvals, Vector &yy
void
MultInitSS::load(int i, double res)
{
- // we can afford it, since the evaluator sets res to exact zero if
- // the term is zero
+ /* we can afford it, since the evaluator sets res to exact zero if the term
+ is zero */
if (res == 0)
return;
// decode i and add to either b or F
diff --git a/dynare++/tl/cc/t_container.cc b/dynare++/tl/cc/t_container.cc
index 4e846e618176ac37f7e4077845d103bd7fe0bd66..c40920153d95639b57a7026444d88fab0bebf72b 100644
--- a/dynare++/tl/cc/t_container.cc
+++ b/dynare++/tl/cc/t_container.cc
@@ -23,7 +23,7 @@
#include "ps_tensor.hh"
#include "pyramid_prod.hh"
-// |UGSContainer| conversion from |FGSContainer|
+// UGSContainer conversion from FGSContainer
UGSContainer::UGSContainer(const FGSContainer &c)
: TensorContainer<UGSTensor>(c.num())
{
@@ -31,17 +31,16 @@ UGSContainer::UGSContainer(const FGSContainer &c)
insert(std::make_unique<UGSTensor>(*(it.second)));
}
-/* We set |l| to dimension of |t|, this is a tensor which multiplies
- tensors from the container from the left. Also we set |k| to a
- dimension of the resulting tensor. We go through all equivalences on
- |k| element set and pickup only those which have $l$ classes.
+/* We set ‘l’ to dimension of ‘t’, this is a tensor which multiplies tensors
+ from the container from the left. Also we set ‘k’ to a dimension of the
+ resulting tensor. We go through all equivalences on ‘k’ element set and
+ pickup only those which have ‘l’ classes.
- In each loop, we fetch all necessary tensors for the product to the
- vector |ts|. Then we form Kronecker product |KronProdAll| and feed it
- with tensors from |ts|. Then we form unfolded permuted symmetry tensor
- |UPSTensor| as matrix product of |t| and Kronecker product |kp|. Then
- we add the permuted data to |out|. This is done by |UPSTensor| method
- |addTo|. */
+ In each loop, we fetch all necessary tensors for the product to the vector
+ ‘ts’. Then we form Kronecker product KronProdAll and feed it with tensors
+ from ‘ts’. Then we form unfolded permuted symmetry tensor UPSTensor as
+ matrix product of ‘t’ and Kronecker product ‘kp’. Then we add the permuted
+ data to ‘out’. This is done by UPSTensor method addTo(). */
void
UGSContainer::multAndAdd(const UGSTensor &t, UGSTensor &out) const
@@ -53,8 +52,7 @@ UGSContainer::multAndAdd(const UGSTensor &t, UGSTensor &out) const
for (const auto &it : eset)
if (it.numClasses() == l)
{
- std::vector<const UGSTensor *> ts
- = fetchTensors(out.getSym(), it);
+ std::vector<const UGSTensor *> ts = fetchTensors(out.getSym(), it);
KronProdAllOptim kp(l);
for (int i = 0; i < l; i++)
kp.setMat(i, *(ts[i]));
@@ -64,7 +62,7 @@ UGSContainer::multAndAdd(const UGSTensor &t, UGSTensor &out) const
}
}
-// |FGSContainer| conversion from |UGSContainer|
+// FGSContainer conversion from UGSContainer
FGSContainer::FGSContainer(const UGSContainer &c)
: TensorContainer<FGSTensor>(c.num())
{
@@ -72,9 +70,9 @@ FGSContainer::FGSContainer(const UGSContainer &c)
insert(std::make_unique<FGSTensor>(*(it.second)));
}
-// |FGSContainer::multAndAdd| folded code
+// FGSContainer::multAndAdd() folded code
/* Here we perform one step of the Faà Di Bruno operation. We call the
- |multAndAdd| for unfolded tensor. */
+ multAndAdd() for unfolded tensor. */
void
FGSContainer::multAndAdd(const FGSTensor &t, FGSTensor &out) const
{
@@ -82,10 +80,9 @@ FGSContainer::multAndAdd(const FGSTensor &t, FGSTensor &out) const
multAndAdd(ut, out);
}
-// |FGSContainer::multAndAdd| unfolded code
-/* This is the same as |@<|UGSContainer::multAndAdd| code@>|
- but we do not construct |UPSTensor| from the Kronecker
- product, but |FPSTensor|. */
+// FGSContainer::multAndAdd() unfolded code
+/* This is the same as UGSContainer::multAndAdd() but we do not construct
+ UPSTensor from the Kronecker product, but FPSTensor. */
void
FGSContainer::multAndAdd(const UGSTensor &t, FGSTensor &out) const
{
@@ -107,10 +104,9 @@ FGSContainer::multAndAdd(const UGSTensor &t, FGSTensor &out) const
}
}
-/* This fills a given vector with integer sequences corresponding to
- first |num| indices from interval |start| (including) to |end|
- (excluding). If there are not |num| of such indices, the shorter vector
- is returned. */
+/* This fills a given vector with integer sequences corresponding to first
+ ‘num’ indices from interval ‘start’ (including) to ‘end’ (excluding). If
+ there are not ‘num’ of such indices, the shorter vector is returned. */
Tensor::index
FGSContainer::getIndices(int num, std::vector<IntSequence> &out,
const Tensor::index &start,
diff --git a/dynare++/tl/cc/tensor.cc b/dynare++/tl/cc/tensor.cc
index 2d1c4fcc59dcb72d26092d7822040359c4cfb711..875876d7d835f4ed1d962c43fe2a4424da41bfcf 100644
--- a/dynare++/tl/cc/tensor.cc
+++ b/dynare++/tl/cc/tensor.cc
@@ -23,10 +23,10 @@
#include "tl_static.hh"
#include "pascal_triangle.hh"
-/* Here we increment a given sequence within full symmetry given by
- |nv|, which is number of variables in each dimension. The underlying
- tensor is unfolded, so we increase the rightmost by one, and if it is
- |nv| we zero it and increase the next one to the left. */
+/* Here we increment a given sequence within full symmetry given by ‘nv’, which
+ is number of variables in each dimension. The underlying tensor is unfolded,
+ so we increase the rightmost by one, and if it is ‘nv’ we zero it and
+ increase the next one to the left. */
void
UTensor::increment(IntSequence &v, int nv)
@@ -42,7 +42,7 @@ UTensor::increment(IntSequence &v, int nv)
}
}
-/* This is dual to |UTensor::increment(IntSequence& v, int nv)|. */
+/* This is dual to UTensor::increment(IntSequence& v, int nv). */
void
UTensor::decrement(IntSequence &v, int nv)
@@ -58,11 +58,10 @@ UTensor::decrement(IntSequence &v, int nv)
}
}
-/* Here we increment index for general symmetry for unfolded
- storage. The sequence |nvmx| assigns for each coordinate a number of
- variables. Since the storage is unfolded, we do not need information
- about what variables are symmetric, everything necessary is given by
- |nvmx|. */
+/* Here we increment index for general symmetry for unfolded storage. The
+ sequence ‘nvmx’ assigns for each coordinate a number of variables. Since the
+ storage is unfolded, we do not need information about what variables are
+ symmetric, everything necessary is given by ‘nvmx’. */
void
UTensor::increment(IntSequence &v, const IntSequence &nvmx)
@@ -78,8 +77,8 @@ UTensor::increment(IntSequence &v, const IntSequence &nvmx)
}
}
-/* This is a dual code to |UTensor::increment(IntSequence& v, const
- IntSequence& nvmx)|. */
+/* This is a dual code to UTensor::increment(IntSequence& v, const IntSequence&
+ nvmx). */
void
UTensor::decrement(IntSequence &v, const IntSequence &nvmx)
@@ -95,8 +94,8 @@ UTensor::decrement(IntSequence &v, const IntSequence &nvmx)
}
}
-/* Here we return an offset for a given coordinates of unfolded full
- symmetry tensor. This is easy. */
+/* Here we return an offset for a given coordinates of unfolded full symmetry
+ tensor. This is easy. */
int
UTensor::getOffset(const IntSequence &v, int nv)
@@ -124,12 +123,11 @@ UTensor::getOffset(const IntSequence &v, const IntSequence &nvmx)
return res;
}
-/* Decrementing of coordinates of folded index is not that easy. Note
- that if a trailing part of coordinates is $(b, a, a, a)$ (for
- instance) with $b<a$, then a preceding coordinates are $(b, a-1, n-1,
- n-1)$, where $n$ is a number of variables |nv|. So we find the left
- most element which is equal to the last element, decrease it by one,
- and then set all elements to the right to $n-1$. */
+/* Decrementing of coordinates of folded index is not that easy. Note that if a
+ trailing part of coordinates is (b,a,a,a) (for instance) with b<a, then a
+ preceding coordinates are (b,a−1,n−1,n−1), where n is a number of variables
+ ‘nv’. So we find the left most element which is equal to the last element,
+ decrease it by one, and then set all elements to the right to n−1. */
void
FTensor::decrement(IntSequence &v, int nv)
diff --git a/dynare++/tl/cc/twod_matrix.hh b/dynare++/tl/cc/twod_matrix.hh
index 7fc51be945b9fc808fdd40fcac745bb8982d9fc6..83e1b27a210dc08d45392425f0601db9d340954b 100644
--- a/dynare++/tl/cc/twod_matrix.hh
+++ b/dynare++/tl/cc/twod_matrix.hh
@@ -20,15 +20,14 @@
// Matrix interface.
-/* Here we make an interface to 2-dimensional matrix defined in the
- Sylvester module. That abstraction provides an interface to BLAS. The
- main purpose of this file is to only make its subclass in order to
- keep the tensor library and Sylvester module independent. So here is
- mainly renaming of methods.
+/* Here we make an interface to 2-dimensional matrix defined in the Sylvester
+ module. That abstraction provides an interface to BLAS. The main purpose of
+ this file is to only make its subclass in order to keep the tensor library
+ and Sylvester module independent. So here is mainly renaming of methods.
- Similarly as in the Sylvester module we declare two classes
- |TwoDMatrix| and |ConstTwoDMatrix|. The only purpose of the latter is
- to allow submatrix construction from const reference arguments. */
+ Similarly as in the Sylvester module we declare two classes TwoDMatrix and
+ ConstTwoDMatrix. The only purpose of the latter is to allow submatrix
+ construction from const reference arguments. */
#ifndef TWOD_MATRIX_H
#define TWOD_MATRIX_H
@@ -42,10 +41,6 @@
class TwoDMatrix;
-/* We make two obvious constructors, and then a constructor making
- submatrix of subsequent columns. We also rename
- |GeneralMatrix::nrows()| and |GeneralMatrix::ncols()|. */
-
class ConstTwoDMatrix : public ConstGeneralMatrix
{
public:
@@ -55,12 +50,18 @@ public:
}
ConstTwoDMatrix(const ConstTwoDMatrix &m) = default;
ConstTwoDMatrix(ConstTwoDMatrix &&m) = default;
+
// Implicit conversion from TwoDMatrix is ok, since it's cheap
ConstTwoDMatrix(const TwoDMatrix &m);
+
+ // Constructors creating a submatrix of consecutive columns
ConstTwoDMatrix(const TwoDMatrix &m, int first_col, int num);
ConstTwoDMatrix(const ConstTwoDMatrix &m, int first_col, int num);
+
+ // Constructors creating a submatrix of consecutive rows
ConstTwoDMatrix(int first_row, int num, const TwoDMatrix &m);
ConstTwoDMatrix(int first_row, int num, const ConstTwoDMatrix &m);
+
ConstTwoDMatrix(const ConstTwoDMatrix &m, int first_row, int first_col, int rows, int cols)
: ConstGeneralMatrix(m, first_row, first_col, rows, cols)
{
@@ -78,12 +79,6 @@ public:
void writeMat(mat_t *fd, const std::string &vname) const;
};
-/* Here we do the same as for |ConstTwoDMatrix| plus define
- methods for copying and adding rows and columns.
-
- Also we have |save| method which dumps the matrix to a file with a
- given name. The file can be read by Scilab {\tt fscanfMat} function. */
-
class TwoDMatrix : public GeneralMatrix
{
public:
@@ -154,7 +149,7 @@ public:
TwoDMatrix &operator=(TwoDMatrix &&m) = default;
TwoDMatrix &operator=(const ConstTwoDMatrix &m);
- // |TwoDMatrix| row methods declarations
+ // TwoDMatrix row methods declarations
void copyRow(int from, int to);
void copyRow(const ConstTwoDMatrix &m, int from, int to);
void
@@ -179,7 +174,7 @@ public:
addRow(d, ConstTwoDMatrix(m), from, to);
}
- // |TwoDMatrix| column methods declarations
+ // TwoDMatrix column methods declarations
void copyColumn(int from, int to);
void copyColumn(const ConstTwoDMatrix &m, int from, int to);
void
@@ -204,7 +199,9 @@ public:
addColumn(d, ConstTwoDMatrix(m), from, to);
}
+ // Saves the matrix to a text file
void save(const std::string &fname) const;
+
void
writeMat(mat_t *fd, const std::string &vname) const
{
diff --git a/dynare++/tl/testing/factory.hh b/dynare++/tl/testing/factory.hh
index 87007b650ff40f478131637ea9f43f705c1d9418..f5eab5d4943209206481470f131c8fb3db31e0b9 100644
--- a/dynare++/tl/testing/factory.hh
+++ b/dynare++/tl/testing/factory.hh
@@ -41,7 +41,7 @@ class Factory
void fillMatrix(TwoDMatrix &m);
public:
double get();
- // this can be used with UGSTensor, FGSTensor
+ // This can be used with UGSTensor, FGSTensor
template <class _Ttype>
std::unique_ptr<_Ttype>
make(int r, const Symmetry &s, const IntSequence &nvs)
@@ -52,7 +52,7 @@ public:
return res;
}
- // this can be used with FFSTensor, UFSTensor, FRTensor, URTensor
+ // This can be used with FFSTensor, UFSTensor, FRTensor, URTensor
template <class _Ttype>
std::unique_ptr<_Ttype>
make(int r, int nv, int dim)
@@ -71,10 +71,10 @@ public:
_Ctype res(symnum);
for (int dim = 1; dim <= maxdim; dim++)
if (symnum == 1)
- // full symmetry
+ // Full symmetry
res.insert(make<_Ttype>(r, Symmetry{dim}, nvs));
else
- // general symmetry
+ // General symmetry
for (int i = 0; i <= dim; i++)
res.insert(make<_Ttype>(r, Symmetry{i, dim-i}, nvs));
return res;
diff --git a/dynare++/tl/testing/monoms.cc b/dynare++/tl/testing/monoms.cc
index 1f31e656bb844cca2e941e3c439afa31dcc35e01..e0690018c4c081f7d394d4de40a84733a6140d43 100644
--- a/dynare++/tl/testing/monoms.cc
+++ b/dynare++/tl/testing/monoms.cc
@@ -168,7 +168,7 @@ Monom2Vector::Monom2Vector(const Monom1Vector &g, const Monom2Vector &xmon)
}
for (int i = 0; i < len; i++)
- // multiply from xmon
+ // Multiply from xmon
for (int j = 0; j < g.nx; j++)
{
int ex = g.x[i].operator[](j);
@@ -285,7 +285,7 @@ Monom4Vector::Monom4Vector(const Monom4Vector &f, const Monom4Vector &bigg,
for (int i = 0; i < len; i++)
{
- // multiply from G (first argument)
+ // Multiply from G (first argument)
for (int j = 0; j < f.nx1; j++)
{
int ex = f.x1[i].operator[](j);
@@ -294,7 +294,7 @@ Monom4Vector::Monom4Vector(const Monom4Vector &f, const Monom4Vector &bigg,
x3[i].multiplyWith(ex, bigg.x3[j]);
x4[i].multiplyWith(ex, bigg.x4[j]);
}
- // multiply from g (second argument)
+ // Multiply from g (second argument)
for (int j = 0; j < f.nx2; j++)
{
int ex = f.x2[i].operator[](j);
@@ -302,9 +302,9 @@ Monom4Vector::Monom4Vector(const Monom4Vector &f, const Monom4Vector &bigg,
x2[i].multiplyWith(ex, g.x2[j]);
x4[i].multiplyWith(ex, g.x4[j]);
}
- // add y as third argument of f
+ // Add y as third argument of f
x1[i].add(1, f.x3[i]);
- // add u as fourth argument of f
+ // Add u as fourth argument of f
x2[i].add(1, f.x4[i]);
}
}
diff --git a/dynare++/tl/testing/monoms.hh b/dynare++/tl/testing/monoms.hh
index 0c1fffe929a7f115784df51edea9e236f4aa7c13..96db16501fe66a72053bbe75a57ddc1564192aa2 100644
--- a/dynare++/tl/testing/monoms.hh
+++ b/dynare++/tl/testing/monoms.hh
@@ -51,7 +51,7 @@ public:
Monom(int len); // generate a random monom
Monom(int len, int item); // generate monom whose items are the given item
double deriv(const IntSequence &vars) const;
- // this = this*m^ex (in monomial sense)
+ // this = this·mᵉˣ (in monomial sense)
void multiplyWith(int ex, const Monom &m);
void print() const;
};
@@ -71,16 +71,15 @@ public:
void print() const;
};
-//class Monom3Vector;
class Monom2Vector
{
int ny, nu;
int len;
std::vector<Monom> y, u;
public:
- // generate random vector of monom two vector
+ // Generate random vector of monom two vector
Monom2Vector(int nyy, int nuu, int l);
- // calculate g(x(y,u))
+ // Calculate g(x(y,u))
Monom2Vector(const Monom1Vector &g, const Monom2Vector &xmon);
~Monom2Vector() = default;
void deriv(const Symmetry &s, const IntSequence &c, Vector &out) const;
@@ -95,13 +94,13 @@ class Monom4Vector
int nx1, nx2, nx3, nx4;
std::vector<Monom> x1, x2, x3, x4;
public:
- /* random for g(y,u,sigma) */
+ /* Random for g(y,u,σ) */
Monom4Vector(int l, int ny, int nu);
- /* random for G(y,u,u',sigma) */
+ /* Random for G(y,u,u′,σ) */
Monom4Vector(int l, int ny, int nu, int nup);
- /* random for f(y+,y,y-,u) */
+ /* Random for f(y⁺,y,y⁻,u) */
Monom4Vector(int l, int nbigg, int ng, int ny, int nu);
- /* substitution f(G(y,u,u',sigma),g(y,u,sigma),y,u) */
+ /* Substitution f(G(y,u,u′,σ),g(y,u,σ),y,u) */
Monom4Vector(const Monom4Vector &f, const Monom4Vector &bigg,
const Monom4Vector &g);
~Monom4Vector() = default;
diff --git a/dynare++/tl/testing/tests.cc b/dynare++/tl/testing/tests.cc
index 1bce947d43f193f16dc92f25b33625a7eac45639..d2f01fe5c1013b622c8a3c8b0c24711de8f2d544 100644
--- a/dynare++/tl/testing/tests.cc
+++ b/dynare++/tl/testing/tests.cc
@@ -400,7 +400,7 @@ TestRunnable::unfold_zcont(int nf, int ny, int nu, int nup, int nbigg,
IntSequence nvs{ny, nu, nup, 1};
double maxnorm = 0.0;
- // form ZContainer
+ // Form ZContainer
UnfoldedZContainer zc(&uG_cont, nbigg, &ug_cont, ng, ny, nu);
for (int d = 2; d <= dim; d++)
@@ -1124,12 +1124,12 @@ main()
all_tests.push_back(std::make_unique<UnfoldZCont>());
/* Initialize library. We do it for each test individually instead of
- computing the maximum dimension and number of variables, because otherwise it
- does not pass the check on maximum problem size. */
+ computing the maximum dimension and number of variables, because otherwise
+ it does not pass the check on maximum problem size. */
for (const auto &test : all_tests)
TLStatic::init(test->dim, test->nvar);
- // launch the tests
+ // Launch the tests
int success = 0;
for (const auto &test : all_tests)
{
diff --git a/dynare++/utils/cc/pascal_triangle.cc b/dynare++/utils/cc/pascal_triangle.cc
index d48922db0b0666b048da36222efb734cd90f8b90..dff2d6b7e7367224c56da3741cb1566d2aace7c9 100644
--- a/dynare++/utils/cc/pascal_triangle.cc
+++ b/dynare++/utils/cc/pascal_triangle.cc
@@ -54,7 +54,7 @@ PascalRow::prolong(const PascalRow &prev)
void
PascalRow::prolongFirst(int n)
{
- // todo: check n = 1;
+ // TODO: check n = 1
for (int i = static_cast<int>(size())+2; i <= n; i++)
push_back(i);
}
@@ -71,7 +71,7 @@ PascalRow::print() const
namespace PascalTriangle
{
- namespace // Anonymous namespace that is a functional equivalent of "private"
+ namespace // Anonymous namespace that is a functional equivalent of “private”
{
std::vector<PascalRow> tr(1);
std::mutex mut; // For protecting the triangle from concurrent modifications
@@ -92,7 +92,7 @@ namespace PascalTriangle
void
ensure(int n, int k)
{
- // add along n
+ // Add along n
if (n > max_n())
{
std::lock_guard<std::mutex> lk{mut};
@@ -115,7 +115,7 @@ namespace PascalTriangle
int
noverk(int n, int k)
{
- // todo: rais if out of bounds
+ // TODO: raise exception if out of bounds
if (n-k < k)
k = n-k;
if (k == 0)
diff --git a/dynare++/utils/cc/sthread.cc b/dynare++/utils/cc/sthread.cc
index 3b8eaec0da5ce09b92d0debfc4f3d38db57b2aa7..303cd7b0009f0dbfc22e090cc2f707731dc3b313 100644
--- a/dynare++/utils/cc/sthread.cc
+++ b/dynare++/utils/cc/sthread.cc
@@ -24,7 +24,7 @@
namespace sthread
{
- /* We set the default value for |max_parallel_threads| to half the number of
+ /* We set the default value for ‘max_parallel_threads’ to half the number of
logical CPUs */
int
default_threads_number()
@@ -34,10 +34,10 @@ namespace sthread
int detach_thread_group::max_parallel_threads = default_threads_number();
- /* We cycle through all threads in the group, and in each cycle we wait
- for the change in the |counter|. If the counter indicates less than
- maximum parallel threads running, then a new thread is run, and the
- iterator in the list is moved.
+ /* We cycle through all threads in the group, and in each cycle we wait for
+ the change in the ‘counter’. If the counter indicates less than maximum
+ parallel threads running, then a new thread is run, and the iterator in
+ the list is moved.
At the end we have to wait for all thread to finish. */
void
@@ -49,7 +49,7 @@ namespace sthread
{
counter++;
std::thread th{[&, it] {
- // The "it" variable is captured by value, because otherwise the iterator may move
+ // The ‘it’ variable is captured by value, because otherwise the iterator may move
(*it)->operator()(mut_threads);
std::unique_lock<std::mutex> lk2{mut_cv};
counter--;
diff --git a/dynare++/utils/cc/sthread.hh b/dynare++/utils/cc/sthread.hh
index a4d01705693a2b7d9ce958b73e94bf60136bdfad..fe52e3d9c6b1be3f891cf1b39fd7efa1512c54eb 100644
--- a/dynare++/utils/cc/sthread.hh
+++ b/dynare++/utils/cc/sthread.hh
@@ -20,25 +20,23 @@
// Simple threads.
-/* This file defines types making a simple interface to
- multi-threading.
+/* This file defines types making a simple interface to multi-threading.
The file provides the following interfaces:
- \unorderedlist
- \li |detach_thread| is a pure virtual class, which must be inherited and a
- method |operator()()| be implemented as the running code of the
- thread.
- \li |detach_thread_group| allows insertion of |detach_thread|s and running
- all of them simultaneously. The threads
- are not joined, they are synchronized by means of a counter counting
- running threads. A change of the counter is checked by waiting on an
- associated condition. The number of maximum parallel
- threads can be controlled. See below. The group also provides a mutex to be
- shared between the workers for their own synchronization purposes.
- \endunorderedlist
- The number of maximum parallel threads is controlled via a static
- member of the |detach_thread_group| class. */
+ — detach_thread is a pure virtual class, which must be inherited and a
+ method operator()() be implemented as the running code of the thread.
+
+ — detach_thread_group allows insertion of detach_thread’s and running all of
+ them simultaneously. The threads are not joined, they are synchronized by
+ means of a counter counting running threads. A change of the counter is
+ checked by waiting on an associated condition. The number of maximum
+ parallel threads can be controlled. See below. The group also provides a
+ mutex to be shared between the workers for their own synchronization
+ purposes.
+
+ The number of maximum parallel threads is controlled via a static member of
+ the detach_thread_group class. */
#ifndef STHREAD_H
#define STHREAD_H
@@ -53,11 +51,6 @@
namespace sthread
{
- /* The detached thread is the same as joinable |thread|. We only
- re-implement |run| method to call |thread_traits::detach_run|, and add
- a method which installs a counter. The counter is increased and
- decreased on the body of the new thread. */
-
class detach_thread
{
public:
@@ -65,10 +58,6 @@ namespace sthread
virtual void operator()(std::mutex &mut) = 0;
};
- /* The detach thread group is (by interface) the same as
- |thread_group|. The extra thing we have here is the |counter|. The
- implementation of |insert| and |run| is different. */
-
class detach_thread_group
{
std::vector<std::unique_ptr<detach_thread>> tlist;