diff --git a/dynare++/tl/cc/pyramid_prod.cc b/dynare++/tl/cc/pyramid_prod.cc index 9d29a8258..83458845f 100644 --- a/dynare++/tl/cc/pyramid_prod.cc +++ b/dynare++/tl/cc/pyramid_prod.cc @@ -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 &ts, diff --git a/dynare++/tl/cc/pyramid_prod.hh b/dynare++/tl/cc/pyramid_prod.hh index fe37671ac..86fe4d003 100644 --- a/dynare++/tl/cc/pyramid_prod.hh +++ b/dynare++/tl/cc/pyramid_prod.hh @@ -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 -/* 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 { diff --git a/dynare++/tl/cc/pyramid_prod2.cc b/dynare++/tl/cc/pyramid_prod2.cc index 87d4ae87b..1dbbd34a6 100644 --- a/dynare++/tl/cc/pyramid_prod2.cc +++ b/dynare++/tl/cc/pyramid_prod2.cc @@ -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 &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) diff --git a/dynare++/tl/cc/pyramid_prod2.hh b/dynare++/tl/cc/pyramid_prod2.hh index bb96019b0..fe23ccd21 100644 --- a/dynare++/tl/cc/pyramid_prod2.hh +++ b/dynare++/tl/cc/pyramid_prod2.hh @@ -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| and coordinate - |c| of the column. + input of this helper class is StackProduct 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 { diff --git a/dynare++/tl/cc/stack_container.hh b/dynare++/tl/cc/stack_container.hh index 5b691c509..3ee5df59b 100644 --- a/dynare++/tl/cc/stack_container.hh +++ b/dynare++/tl/cc/stack_container.hh @@ -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.