stepan
/
dynare
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Dynare++: improvements to comments 

[skip ci]
time-shift
Sébastien Villemot 2019-04-05 18:39:47 +02:00
parent bdc95f23aa
commit e67c172000
No known key found for this signature in database
GPG Key ID: 2CECE9350ECEBE4A
5 changed files with 125 additions and 141 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

42
dynare++/tl/cc/pyramid_prod.cc
View File

@ -4,17 +4,16 @@
#include "permutation.hh"
#include "tl_exception.hh"
/* Here we construct the |USubTensor| object. We allocate space via the
parent |URTensor|. Number of columns is a length of the list of
indices |lst|, number of variables and dimensions are of the tensor
$h$, this is given by |hdims|.
/* Here we construct the USubTensor object. We allocate space via the parent
URTensor. Number of columns is a length of the list of indices ‘lst’, number
of variables and dimensions are of the tensor $h$, this is given by ‘hdims’.
We go through all equivalences with number of classes equal to
dimension of $B$. For each equivalence we make a permutation
|per|. Then we fetch all the necessary tensors $g$ with symmetries
implied by symmetry of $B$ and the equivalence. Then we go through the
list of indices, permute them by the permutation and add the Kronecker
product of the selected columns. This is done by |addKronColumn|. */
We go through all equivalences with number of classes equal to dimension of
B. For each equivalence we make a permutation ‘per’. Then we fetch all the
necessary tensors g with symmetries implied by symmetry of B and the
equivalence. Then we go through the list of indices, permute them by the
permutation and add the Kronecker product of the selected columns. This is
done by addKronColumn(). */
USubTensor::USubTensor(const TensorDimens &bdims,
const TensorDimens &hdims,
@ -43,19 +42,18 @@ USubTensor::USubTensor(const TensorDimens &bdims,
}
}
/* This makes a Kronecker product of appropriate columns from tensors
in |fs| and adds such data to |i|-th column of this matrix. The
appropriate columns are defined by |pindex| sequence. A column of a
tensor has index created from a corresponding part of |pindex|. The
sizes of these parts are given by dimensions of the tensors in |ts|.
/* This makes a Kronecker product of appropriate columns from tensors in ‘fs’
and adds such data to i-th column of this matrix. The appropriate columns
are defined by ‘pindex’ sequence. A column of a tensor has index created
from a corresponding part of ‘pindex’. The sizes of these parts are given by
dimensions of the tensors in ‘ts’.
Here we break the given index |pindex| according to the dimensions of
the tensors in |ts|, and for each subsequence of the |pindex| we find
an index of the folded tensor, which involves calling |getOffset| for
folded tensor, which might be costly. We gather all columns to a
vector |tmpcols| which are Kronecker multiplied in constructor of
|URSingleTensor|. Finally we add data of |URSingleTensor| to the
|i|-th column. */
Here we break the given index ‘pindex’ according to the dimensions of the
tensors in ‘ts’, and for each subsequence of the ‘pindex’ we find an index
of the folded tensor, which involves calling getOffset() for folded tensor,
which might be costly. We gather all columns to a vector ‘tmpcols’ which are
Kronecker multiplied in constructor of URSingleTensor. Finally we add data
of URSingleTensor to the i-th column. */
void
USubTensor::addKronColumn(int i, const std::vector<const FGSTensor *> &ts,

67
dynare++/tl/cc/pyramid_prod.hh
View File

@ -4,45 +4,49 @@
/* In here, we implement the Faà Di Bruno for folded
tensors. Recall, that one step of the Faà Di Bruno is a formula:
$$\left[B_{s^k}\right]_{\alpha_1\ldots\alpha_k}=
[h_{y^l}]_{\gamma_1\ldots\gamma_l}
\prod_{m=1}^l\left[g_{s^{\vert c_m\vert}}\right]^{\gamma_m}_{c_m(\alpha)}
$$
In contrast to unfolded implementation of |UGSContainer::multAndAdd|
with help of |KronProdAll| and |UPSTensor|, we take a completely
ₗ
[B_sᵏ]_α₁…αₗ = [h_yˡ]_γ₁…γₗ ∏ [g_(s^|cₘ|)]_cₘ(α)^γₘ
ᵐ⁼¹
In contrast to unfolded implementation of UGSContainer::multAndAdd()
with help of KronProdAll and UPSTensor, we take a completely
different strategy. We cannot afford full instantiation of
$$\sum_{c\in M_{l,k}}
\prod_{m=1}^l\left[g_{s^{\vert c_m\vert}}\right]^{\gamma_m}_{c_m(\alpha)}$$
ₗ
∑ ∏ [g_(s^|cₘ|)]_cₘ(α)^γₘ
c∈ℳₗ,ₖ ᵐ⁼¹
and therefore we do it per partes. We select some number of columns,
for instance 10, calculate 10 continuous iterators of tensor $B$. Then we
for instance 10, calculate 10 continuous iterators of tensor B. Then we
form unfolded tensor
$$[G]_S^{\gamma_1\ldots\gamma_l}=\left[\sum_{c\in M_{l,k}}
\prod_{m=1}^l\left[g_{s^{\vert c_m\vert}}\right]^{\gamma_m}_{c_m(\alpha)}
\right]_S$$
where $S$ is the selected set of 10 indices. This is done as Kronecker
⎡ ₗ ⎤
[G]_S^γ₁…γₗ = ⎢ ∑ ∏ [g_(s^|cₘ|)]_cₘ(α)^γₘ⎥
⎣c∈ℳₗ,ₖ ᵐ⁼¹ ⎦_S
where S is the selected set of 10 indices. This is done as Kronecker
product of vectors corresponding to selected columns. Note that, in
general, there is no symmetry in $G$, its type is special class for
general, there is no symmetry in G, its type is special class for
this purpose.
If $g$ is folded, then we have to form folded version of $G$. There is
no symmetry in $G$ data, so we sum all unfolded indices corresponding
If g is folded, then we have to form the folded version of G. There is
no symmetry in G data, so we sum all unfolded indices corresponding
to folded index together. This is perfectly OK, since we multiply
these groups of (equivalent) items with the same number in fully
symmetric $g$.
symmetric g.
After this, we perform ordinary matrix multiplication to obtain a
selected set of columns of $B$.
selected set of columns of B.
In here, we define a class for forming and representing
$[G]_S^{\gamma_1\ldots\gamma_l}$. Basically, this tensor is
row-oriented (multidimensional index is along rows), and it is fully
symmetric. So we inherit from |URTensor|. If we need its folded
version, we simply use a suitable conversion. The new abstraction will
have only a new constructor allowing a construction from the given set
of indices $S$, and given set of tensors $g$. The rest of the process
is implemented in |@<|FGSContainer::multAndAdd| unfolded code@>| or
|@<|FGSContainer::multAndAdd| folded code@>|. */
In here, we define a class for forming and representing [G]_S^γ₁…γₗ.
Basically, this tensor is row-oriented (multidimensional index is along
rows), and it is fully symmetric. So we inherit from URTensor. If we need
its folded version, we simply use a suitable conversion. The new abstraction
will have only a new constructor allowing a construction from the given set
of indices S, and given set of tensors g. The rest of the process is
implemented in FGSContainer::multAndAdd() unfolded code or
FGSContainer::multAndAdd() folded code. */
#ifndef PYRAMID_PROD_H
#define PYRAMID_PROD_H
@ -54,11 +58,10 @@
#include <vector>
/* Here we define the new tensor for representing
$[G]_S^{\gamma_1\ldots\gamma_l}$. It allows a construction from
container of folded general symmetry tensors |cont|, and set of
indices |ts|. Also we have to supply dimensions of resulting tensor
$B$, and dimensions of tensor $h$. */
/* Here we define the new tensor for representing [G]_S^γ₁…γₗ. It allows a
construction from container of folded general symmetry tensors ‘cont’, and
set of indices ‘ts’. Also we have to supply dimensions of resulting tensor
B, and dimensions of tensor h. */
class USubTensor : public URTensor
{

16
dynare++/tl/cc/pyramid_prod2.cc
View File

@ -3,9 +3,9 @@
#include "pyramid_prod2.hh"
#include "rfs_tensor.hh"
/* Here we only call |sp.createPackedColumns(c, cols, unit_flag)| which
fills |cols| and |unit_flag| for the given column |c|. Then we set
|end_seq| according to |unit_flag| and columns lengths. */
/* Here we only call sp.createPackedColumns(c, cols, unit_flag) which
fills ‘cols’ and ‘unit_flag’ for the given column ‘c’. Then we set
‘end_seq’ according to ‘unit_flag’ and columns lengths. */
IrregTensorHeader::IrregTensorHeader(const StackProduct<FGSTensor> &sp,
const IntSequence &c)
@ -36,10 +36,10 @@ IrregTensorHeader::increment(IntSequence &v) const
return;
int i = v.size()-1;
// increment |i|-th item in coordinate |v|
// increment i-th item in coordinate ‘v’
/* Here we increment item of coordinates. Whenever we reached end of
column coming from matrices, and |unit_flag| is not $-1$, we have to
jump to that |unit_flag|. */
column coming from matrices, and ‘unit_flag’ is not -1, we have to
jump to that ‘unit_flag’. */
v[i]++;
if (unit_flag[i] != -1 && v[i] == cols[i]->length()-1)
v[i] = unit_flag[i];
@ -48,7 +48,7 @@ IrregTensorHeader::increment(IntSequence &v) const
{
v[i] = 0;
i--;
// increment |i|-th item in coordinate |v|
// increment i-th item in coordinate ‘v’
/* Same code as above */
v[i]++;
if (unit_flag[i] != -1 && v[i] == cols[i]->length()-1)
@ -67,7 +67,7 @@ IrregTensorHeader::calcMaxOffset() const
return res;
}
/* Everything is done in |IrregTensorHeader|, only we have to Kronecker
/* Everything is done in IrregTensorHeader, only we have to Kronecker
multiply all columns of the header. */
IrregTensor::IrregTensor(const IrregTensorHeader &h)

127
dynare++/tl/cc/pyramid_prod2.hh
View File

@ -3,67 +3,51 @@
// Multiplying stacked tensor columns.
/* We need to calculate the following tensor product:
$$\left[f_{s^j}\right]_{\alpha_1\ldots\alpha_j}=
\sum_{l=1}^j\left[f_{z^l}\right]_{\beta_1\ldots\beta_l}
\sum_{c\in M_{l,j}}\prod_{m=1}^l\left[z_{c_m}\right]^{\beta_m}_{c_m(\alpha)}
$$
where $s=[y,u,u',\sigma]$, and $z$ is a composition of four variables,
say $[v,w,y,u]$. Note that $z$ ends with $y$ and $u$, and the only
non-zero derivative of the trailing part of $z$ involving $y$ or $u$
is the first derivative and is the unit matrix $y_y=[1]$ or
$u_u=[1]$. Also, we suppose that the dependence of $v$, and $w$ on $s$
is such that whenever derivative of $w$ is nonzero, then also of
$v$. This means that there for any derivative and any index there is a
continuous part of derivatives of $v$ and optionally of $w$ followed by
column of zeros containing at most one $1$.
ⱼ ₗ
[f_sʲ]_α₁…αⱼ = ∑ [f_zˡ]_β₁…βₗ ∑ ∏ [z_(s^|cₘ|)]_cₘ(α)^βₘ
ˡ⁼¹ c∈ℳₗ,ₖ ᵐ⁼¹
where s=(y,u,u′,σ), and z is a composition of four variables, say (v,w,y,u).
Note that z ends with y and u, and the only non-zero derivative of the
trailing part of z involving y or u is the first derivative and is the unit
matrix y_y=(1) or uᵤ=(1). Also, we suppose that the dependence of v, and w
on s is such that whenever derivative of w is nonzero, then also of v. This
means that there for any derivative and any index there is a continuous part
of derivatives of v and optionally of w followed by column of zeros
containing at most one 1.
This structure can be modelled and exploited with some costs at
programming. For example, let us consider the following product:
$$\left[B_{y^2u^3}\right]_{\alpha_1\alpha_2\beta_1\beta_2\beta_3}=
\ldots
\left[f_{z^3}\right]_{\gamma_1\gamma_2\gamma_3}
\left[z_{yu}\right]^{\gamma_1}_{\alpha_1\beta_1}
\left[z_{y}\right]^{\gamma_2}_{\alpha_2}
\left[z_{uu}\right]^{\gamma_3}_{\beta_2\beta_3}
\ldots$$
The term corresponds to equivalence $\{\{0,2\},\{1\},\{3,4\}\}$. For
the fixed index $\alpha_1\alpha_2\beta_1\beta_2\beta_3$ we have to
make a Kronecker product of the columns
$$
\left[z_{yu}\right]_{\alpha_1\beta_1}\otimes
\left[z_{y}\right]_{\alpha_2}\otimes
\left[z_{uu}\right]_{\beta_2\beta_3}
$$
which can be written as
$$
\left[\matrix{\left[v_{yu}\right]_{\alpha_1\beta_1}\cr
\left[w_{yu}\right]_{\alpha_1\beta_1}\cr 0\cr 0}\right]\otimes
\left[\matrix{\left[v_y\right]_{\alpha_2\vphantom{(}}\cr
\left[w_y\right]_{\alpha_2}\cr 1_{\alpha_2}\cr 0}\right]\otimes
\left[\matrix{\left[v_{uu}\right]_{\beta_2\beta_3\vphantom{(}}\cr
\left[w_{uu}\right]_{\beta_2\beta_3}\cr 0\cr 0}\right]
$$
where $1_{\alpha_2}$ is a column of zeros having the only $1$ at
$\alpha_2$ index.
[B_y²u³]_α₁α₂β₁β₂β₃ = … + [f_z³]_γ₁γ₂γ₃ [z_yu]^γ₁_α₁β₁ [z_y]^γ₂_α₂ [zᵤᵤ]^γ₃_β₂β₃ + …
The term corresponds to equivalence { {0,2}, {1}, {3,4} }. For the fixed
index α₁α₂β₁β₂β₃ we have to make a Kronecker product of the columns
[z_yu]_α₁β₁ ⊗ [z_y]_α₂ ⊗ [zᵤᵤ]_β₂β₃
which can be written as:
⎛[v_yu]_α₁β₁⎞ ⎛[v_y]_α₂⎞ ⎛[vᵤᵤ]_β₂β₃⎞
⎢[w_yu]_α₁β₁⎥ ⎢[w_y]_α₂⎥ ⎢[wᵤᵤ]_β₂β₃⎥
⎢ 0 ⎥⊗⎢ 1_α₂ ⎥⊗⎢ 0 ⎥
⎝ 0 ⎠ ⎝ 0 ⎠ ⎝ 0 ⎠
where 1_α₂ is a column of zeros having the only 1 at α₂ index.
This file develops the abstraction for this Kronecker product column
without multiplication of the zeros at the top. Basically, it will be
a column which is a Kronecker product of the columns without the
zeros:
$$
\left[\matrix{\left[v_{yu}\right]_{\alpha_1\beta_1}\cr
\left[w_{yu}\right]_{\alpha_1\beta_1}}\right]\otimes
\left[\matrix{\left[v_y\right]_{\alpha_2}\cr
\left[w_y\right]_{\alpha_2}\cr 1}\right]\otimes
\left[\matrix{\left[v_{uu}\right]_{\beta_2\beta_3}\cr
\left[w_{uu}\right]_{\beta_2\beta_3}}\right]
$$
The class will have a tensor infrastructure introducing |index| which
iterates over all items in the column with $\gamma_1\gamma_2\gamma_3$
as coordinates in $\left[f_{z^3}\right]$. The data of such a tensor is
⎛[v_y]_α₂⎞
⎛[v_yu]_α₁β₁⎞⊗ ⎢[w_y]_α₂⎥⊗⎛[vᵤᵤ]_β₂β₃⎞
⎝[w_yu]_α₁β₁⎠ ⎝ 1 ⎠ ⎝[wᵤᵤ]_β₂β₃⎠
The class will have a tensor infrastructure introducing ‘index’ which
iterates over all items in the column with γ₁γ₂γ₃
as coordinates in [f_z³]. The data of such a tensor is
not suitable for any matrix operation and will have to be accessed
only through the |index|. Note that this does not matter, since
$\left[f_{z^l}\right]$ are sparse. */
only through the ‘index’. Note that this does not matter, since
[f_zˡ] are sparse. */
#ifndef PYRAMID_PROD2_H
#define PYRAMID_PROD2_H
@ -80,14 +64,14 @@
/* First we declare a helper class for the tensor. Its purpose is to
gather the columns which are going to be Kronecker multiplied. The
input of this helper class is |StackProduct<FGSTensor>| and coordinate
|c| of the column.
input of this helper class is StackProduct<FGSTensor> and coordinate
‘c’ of the column.
It maintains |unit_flag| array which says for what columns we must
stack 1 below $v$ and $w$. In this case, the value of |unit_flag| is
an index of the $1$, otherwise the value of |unit_flag| is -1.
It maintains ‘unit_flag’ array which says for what columns we must
stack 1 below v and w. In this case, the value of ‘unit_flag’ is
an index of the 1, otherwise the value of ‘unit_flag’ is -1.
Also we have storage for the stacked columns |cols|. The object is
Also we have storage for the stacked columns ‘cols’. The object is
responsible for memory management associated to this storage. That is
why we do not allow any copy constructor, since we need to be sure
that no accidental copies take place. We declare the copy constructor
@ -112,21 +96,20 @@ public:
int calcMaxOffset() const;
};
/* Here we declare the irregular tensor. There is no special logic
here. We inherit from |Tensor| and we must implement three methods,
|increment|, |decrement| and |getOffset|. The last two are not
implemented now, since they are not needed, and they raise an
exception. The first just calls |increment| of the header. Also we
declare a method |addTo| which adds this unfolded irregular single
column tensor to folded (regular) single column tensor.
/* Here we declare the irregular tensor. There is no special logic here. We
inherit from Tensor and we must implement three methods, increment(),
decrement() and getOffset(). The last two are not implemented now, since
they are not needed, and they raise an exception. The first just calls
increment() of the header. Also we declare a method addTo() which adds this
unfolded irregular single column tensor to folded (regular) single column
tensor.
The header |IrregTensorHeader| lives with an object by a
reference. This is dangerous. However, we will use this class only in
a simple loop and both |IrregTensor| and |IrregTensorHeader| will be
destructed at the end of a block. Since the super class |Tensor| must
be initialized before any member, we could do either a save copy of
|IrregTensorHeader|, or relatively dangerous the reference member. For
the reason above we chose the latter. */
The header IrregTensorHeader lives with an object by a reference. This is
dangerous. However, we will use this class only in a simple loop and both
IrregTensor and IrregTensorHeader will be destructed at the end of a block.
Since the super class Tensor must be initialized before any member, we could
do either a save copy of IrregTensorHeader, or relatively dangerous the
reference member. For the reason above we chose the latter. */
class IrregTensor : public Tensor
{

14
dynare++/tl/cc/stack_container.hh
View File

@ -5,16 +5,16 @@
/* Here we develop abstractions for stacked containers of tensors. For
instance, in perturbation methods for DSGE models, we need the function:
⎛ G(y*,u,u′,σ) ⎞
z(y*,u,u′,σ) = ⎢ g(y*,u,σ) ⎥
⎢ y* ⎥
⎝ u ⎠
⎛ G(y*,u,u′,σ) ⎞
z(y*,u,u′,σ) = ⎢ g(y*,u,σ) ⎥
⎢ y* ⎥
⎝ u ⎠
and we need to calculate one step of Faà Di Bruno formula:
ₗ
[B_sᵏ]_α₁,…,αₗ = [f_zˡ]_β₁,…,βₗ ∑ ∏ [z_(s^|cₘ|)]_cₘ(α)^βₘ
c∈ℳₗ,ₖ ᵐ⁼¹
ₗ
[B_sᵏ]_α₁…αₗ = [f_zˡ]_β₁…βₗ ∑ ∏ [z_(s^|cₘ|)]_cₘ(α)^βₘ
c∈ℳₗ,ₖ ᵐ⁼¹
where we have containers for derivatives of G and g.