diff --git a/dynare++/tl/cc/normal_moments.cc b/dynare++/tl/cc/normal_moments.cc index 48f6f14e4..33be11ba3 100644 --- a/dynare++/tl/cc/normal_moments.cc +++ b/dynare++/tl/cc/normal_moments.cc @@ -14,14 +14,12 @@ UNormalMoments::UNormalMoments(int maxdim, const TwoDMatrix &v) generateMoments(maxdim, v); } -/* Here we fill up the container with the tensors for $d=2,4,6,\ldots$ - up to the given dimension. Each tensor of moments is equal to - $F_n\left(\otimes^nv\right).$ This has a dimension equal to - $2n$. See the header file for proof and details. - - Here we sequentially construct the Kronecker power - $\otimes^nv$, and apply $F_n$. */ +/* Here we fill up the container with the tensors for d=2,4,6,… up to the given + dimension. Each tensor of moments is equal to Fₙ(⊗ⁿv). This has a dimension + equal to 2n. See the header file for proof and details. + Here we sequentially construct the Kronecker power ⊗ⁿv and apply Fₙ. + */ void UNormalMoments::generateMoments(int maxdim, const TwoDMatrix &v) { @@ -42,15 +40,15 @@ UNormalMoments::generateMoments(int maxdim, const TwoDMatrix &v) newkronv->getData()); kronv = std::move(newkronv); auto mom = std::make_unique(nv, d); - // apply $F_n$ to |kronv| + // apply Fₙ to ‘kronv’ /* Here we go through all equivalences, select only those having 2 - elements in each class, then go through all elements in |kronv| and - add to permuted location of |mom|. + elements in each class, then go through all elements in ‘kronv’ and + add to permuted location of ‘mom’. The permutation must be taken as inverse of the permutation implied by the equivalence, since we need a permutation which after application to identity of indices yileds indices in the equivalence classes. Note - how the |Equivalence::apply| method works. */ + how the Equivalence::apply() method works. */ mom->zeros(); const EquivalenceSet eset = TLStatic::getEquiv(d); for (const auto &cit : eset) @@ -70,7 +68,7 @@ UNormalMoments::generateMoments(int maxdim, const TwoDMatrix &v) } } -/* We return |true| for an equivalence whose each class has 2 elements. */ +/* We return true for an equivalence whose each class has 2 elements. */ bool UNormalMoments::selectEquiv(const Equivalence &e) diff --git a/dynare++/tl/cc/normal_moments.hh b/dynare++/tl/cc/normal_moments.hh index 65ec279f6..d148c3706 100644 --- a/dynare++/tl/cc/normal_moments.hh +++ b/dynare++/tl/cc/normal_moments.hh @@ -11,94 +11,89 @@ ∂f(t)/∂t = f(t)·Vt ∂²f(t)/∂t² = f(t)·(Vt⊗v) - ∂³f(t)/∂t³ = f(t)·(Vt⊗Vt⊗Vt + P(v⊗Vt) + P(Vt⊗v) + v⊗Vt) - ∂⁴f(t)/∂t⁴ = f(t)·(Vt⊗Vt⊗Vt⊗Vt + S(v⊗Vt⊗Vt) + S(Vt⊗v⊗Vt) + S(Vt⊗Vt⊗v) + S(v⊗v)) + ∂³f(t)/∂t³ = f(t)·(Vt⊗Vt⊗Vt + P_?(v⊗Vt) + P_?(Vt⊗v) + v⊗Vt) + ∂⁴f(t)/∂t⁴ = f(t)·(Vt⊗Vt⊗Vt⊗Vt + S_?(v⊗Vt⊗Vt) + S_?(Vt⊗v⊗Vt) + S_?(Vt⊗Vt⊗v) + S_?(v⊗v)) - where v is vectorized V (v=vec(V)), and P is a suitable row permutation + where v is vectorized V (v=vec(V)), and P_? is a suitable row permutation (corresponds to permutation of multidimensional indices) which permutes the tensor data, so that the index of a variable being derived would be the last. This ensures that all (permuted) tensors can be summed yielding a tensor whose indices have some order (in here we chose the order that more - recent derivating variables are to the right). Finally, S is a suitable sum - of various P. + recent derivating variables are to the right). Finally, S_? is a suitable sum + of various P_?. - We are interested in S multiplying the Kronecker powers ⊗ⁿv. The S is a + We are interested in S_? multiplying the Kronecker powers ⊗ⁿv. The S_? is a (possibly) multi-set of permutations of even order. Note that we know the - number of permutations in S. The above formulas for f(t) derivatives are + number of permutations in S_?. The above formulas for f(t) derivatives are valid also for monomial u, and from literature we know that 2n-th moment is - (2n!)/(n!2ⁿ)·σ². So there are (2n!)/(n!2ⁿ) permutations in S. + (2n!)/(n!2ⁿ)·σ². So there are (2n!)/(n!2ⁿ) permutations in S_?. - In order to find the $S_?$ we need to define a couple of - things. First we define a sort of equivalence between the permutations - applicable to even number of indices. We write $P_1\equiv P_2$ - whenever $P_1^{-1}\circ P_2$ permutes only whole pairs, or items - within pairs, but not indices across the pairs. For instance the - permutations $(0,1,2,3)$ and $(3,2,0,1)$ are equivalent, but - $(0,2,1,3)$ is not equivalent with the two. Clearly, the $\equiv$ is - an equivalence. + In order to find the S_? permutation we need to define a couple of things. + First we define a sort of equivalence between the permutations applicable to + even number of indices. We write P₁ ≡ P₂ whenever P₁⁻¹∘P₂ permutes only whole + pairs, or items within pairs, but not indices across the pairs. For instance + the permutations (0,1,2,3) and (3,2,0,1) are equivalent, but (0,2,1,3) is + not equivalent with the two. Clearly, the ≡ relationship is an equivalence. - This allows to define a relation $\sqsubseteq$ between the permutation - multi-sets $S$, which is basically the subset relation $\subseteq$ but - with respect to the equivalence $\equiv$, more formally: - $$S_1\sqsubseteq S_2\quad\hbox{iff}\quad P\in S_1 - \Rightarrow\exists Q\in S_2:P\equiv Q$$ - This induces an equivalence $S_1\equiv S_2$. + This allows to define a relation ⊑ between the permutation multi-sets S, + which is basically the subset relation ⊆ but with respect to the equivalence + relation ≡, more formally: - Now let $F_n$ denote a set of permutations on $2n$ indices which is - maximal with respect to $\sqsubseteq$, and minimal with respect to - $\equiv$. (In other words, it contains everything up to the - equivalence $\equiv$.) It is straightforward to calculate a number of - permutations in $F_n$. This is a total number of all permutations of - $2n$ divided by permutations of pairs divided by permutations within - the pairs. This is ${(2n!)\over n!2^n}$. + S₁ ⊑ S₂ iff P ∈ S₁ ⇒ ∃Q ∈ S₂ : P ≡ Q - We prove that $S_?\equiv F_n$. Clearly $S_?\sqsubseteq F_n$, since - $F_n$ is maximal. In order to prove that $F_n\sqsubseteq S_?$, let us - assert that for any permutation $P$ and for any (semi)positive - definite matrix $V$ we have $PS_?\otimes^nv=S_?\otimes^nv$. Below we - show that there is a positive definite matrix $V$ of some dimension - that for any two permutation multi-sets $S_1$, $S_2$, we have - $$S_1\not\equiv S_2\Rightarrow S_1(\otimes^nv)\neq S_2(\otimes^nv)$$ - So it follows that for any permutation $P$, we have $PS_?\equiv - S_?$. For a purpose of contradiction let $P\in F_n$ be a permutation - which is not equivalent to any permutation from $S_?$. Since $S_?$ is - non-empty, let us pick $P_0\in S_?$. Now assert that - $P_0^{-1}S_?\not\equiv P^{-1}S_?$ since the first contains an identity - and the second does not contain a permutation equivalent to - identity. Thus we have $(P\circ P_0^{-1})S_?\not\equiv S_?$ which - gives the contradiction and we have proved that $F_n\sqsubseteq - S_?$. Thus $F_n\equiv S_?$. Moreover, we know that $S_?$ and $F_n$ - have the same number of permutations, hence the minimality of $S_?$ - with respect to $\equiv$. + This induces an equivalence S₁ ≡ S₂. - Now it suffices to prove that there exists a positive definite $V$ - such that for any two permutation multi-sets $S_1$, and $S_2$ holds - $S_1\not\equiv S_2\Rightarrow S_1(\otimes^nv)\neq S_2(\otimes^nv)$. If - $V$ is $n\times n$ matrix, then $S_1\not\equiv S_2$ implies that there - is identically nonzero polynomial of elements from $V$ of order $n$ - over integers. If $V=A^TA$ then there is identically non-zero - polynomial of elements from $A$ of order $2n$. This means, that we - have to find $n(n+1)/2$ tuple $x$ of real numbers such that all - identically non-zero polynomials $p$ of order $2n$ over integers yield - $p(x)\neq 0$. + Now let Fₙ denote a set of permutations on 2n indices which is maximal with + respect to ⊑, and minimal with respect to ≡ (in other words, it contains + everything up to the equivalence ≡). It is straightforward to calculate a + number of permutations in Fₙ. This is a total number of all permutations + of 2n divided by permutations of pairs divided by permutations within the + pairs. This is (2n!)/(n!2ⁿ). - The $x$ is constructed as follows: $x_i = \pi^{\log{r_i}}$, where $r_i$ - is $i$-th prime. Let us consider monom $x_1^{j_1}\cdot\ldots\cdot - x_k^{j_k}$. When the monom is evaluated, we get - $$\pi^{\log{r_1^{j_1}}+\ldots+\log{r_k^{j_k}}}= - \pi^{\log{\left(r_1^{j_1}\cdot\ldots\cdot r_k^{j_k}\right)}}$$ - Now it is easy to see that if an integer combination of such terms is - zero, then the combination must be either trivial or sum to $0$ and - all monoms must be equal. Both cases imply a polynomial identically - equal to zero. So, any non-trivial integer polynomial evaluated at $x$ - must be non-zero. + We now prove that S_? ≡ Fₙ. Clearly S_? ⊑ Fₙ, since Fₙ is maximal. In order + to prove that Fₙ ⊑ S_?, let us assert that for any permutation P and for + any (semi-)positive definite matrix V we have + P·S_?(⊗ⁿv)=S_?(⊗ⁿv). Below we show that there is a positive + definite matrix V of some dimension that for any two permutation + multi-sets S₁, S₂, we have + + S₁ ≢ S₂ ⇒ S₁(⊗ⁿv) ≠ S₂(⊗ⁿv) + + So it follows that for any permutation P, we have P·S_? ≡ S_?. For a purpose + of contradiction let P ∈ Fₙ be a permutation which is not equivalent to any + permutation from S_?. Since S_? is non-empty, let us pick P₀ ∈ S_?. Now + assert that P₀⁻¹S_? ≢ P⁻¹S_? since the first contains an identity and the + second does not contain a permutation equivalent to identity. Thus we have + (P∘P₀⁻¹)S_? ≢ S_? which gives the contradiction and we have proved that + Fₙ ⊑ S_?. Thus Fₙ ≡ S_?. Moreover, we know that S_? and Fₙ have the same + number of permutations, hence the minimality of S_? with respect to ≡. + + Now it suffices to prove that there exists a positive definite V such that + for any two permutation multi-sets S₁ and S₂ holds S₁ ≢ S₂ ⇒ S₁(⊗ⁿv) ≠ S₂⊗ⁿv. + If V is a n×n matrix, then S₁ ≢ S₂ implies that there is + identically nonzero polynomial of elements from V of order n over + integers. If V=AᵀA then there is identically non-zero polynomial of + elements from A of order 2n. This means, that we have to find n(n+1)/2 + tuple x of real numbers such that all identically non-zero polynomials p + of order 2n over integers yield p(x)≠0. + + The x is constructed as follows: x_i = π^log(rᵢ), where rᵢ is i-th prime. + Let us consider the monom x₁^j₁·…·xₖ^jₖ. When the monom is evaluated, we get + + π^{log(r₁^j₁)+…+log(rₖ^jₖ)}= π^log(r₁^j₁+…+rₖ^jₖ) + + Now it is easy to see that if an integer combination of such terms is zero, + then the combination must be either trivial or sum to 0 and all monoms + must be equal. Both cases imply a polynomial identically equal to zero. So, + any non-trivial integer polynomial evaluated at x must be non-zero. So, having this result in hand, now it is straightforward to calculate - higher moments of normal distribution. Here we define a container, - which does the job. In its constructor, we simply calculate Kronecker - powers of $v$ and apply $F_n$ to $\otimes^nv$. $F_n$ is, in fact, a - set of all equivalences in sense of class |Equivalence| over $2n$ - elements, having $n$ classes each of them having exactly 2 elements. */ + higher moments of normal distribution. Here we define a container, which + does the job. In its constructor, we simply calculate Kronecker powers of v + and apply Fₙ to ⊗ⁿv. Fₙ is, in fact, a set of all equivalences in sense of + class Equivalence over 2n elements, having n classes each of them having + exactly 2 elements. +*/ #ifndef NORMAL_MOMENTS_H #define NORMAL_MOMENTS_H diff --git a/dynare++/tl/cc/permutation.cc b/dynare++/tl/cc/permutation.cc index 0afbe6c3d..831b7186a 100644 --- a/dynare++/tl/cc/permutation.cc +++ b/dynare++/tl/cc/permutation.cc @@ -3,7 +3,7 @@ #include "permutation.hh" #include "tl_exception.hh" -/* This is easy, we simply apply the map in the fashion $s\circ m$.. */ +/* This is easy, we simply apply the map in the fashion s∘m */ void Permutation::apply(const IntSequence &src, IntSequence &tar) const @@ -42,13 +42,13 @@ Permutation::tailIdentity() const } /* This calculates a map which corresponds to sorting in the following - sense: $(\hbox{sorted }s)\circ m = s$, where $s$ is a given sequence. + sense: $(sorted s)∘m = s, where s is a given sequence. - We go through |s| and find an the same item in sorted |s|. We - construct the |permap| from the found pair of indices. We have to be - careful, to not assign to two positions in |s| the same position in - sorted |s|, so we maintain a bitmap |flag|, in which we remember - indices from the sorted |s| already assigned. */ + We go through ‘s’ and find an the same item in sorted ‘s’. We + construct the ‘permap’ from the found pair of indices. We have to be + careful, to not assign to two positions in ‘s’ the same position in + sorted ‘s’, so we maintain a bitmap ‘flag’, in which we remember + indices from the sorted ‘s’ already assigned. */ void Permutation::computeSortingMap(const IntSequence &s) diff --git a/dynare++/tl/cc/permutation.hh b/dynare++/tl/cc/permutation.hh index 36b4a3f08..53d772621 100644 --- a/dynare++/tl/cc/permutation.hh +++ b/dynare++/tl/cc/permutation.hh @@ -5,30 +5,27 @@ /* The permutation class is useful when describing a permutation of indices in permuted symmetry tensor. This tensor comes to existence, for instance, as a result of the following tensor multiplication: - $$\left[g_{y^3}\right]_{\gamma_1\gamma_2\gamma_3} - \left[g_{yu}\right]^{\gamma_1}_{\alpha_1\beta_3} - \left[g_{yu}\right]^{\gamma_2}_{\alpha_2\beta_1} - \left[g_u\right]^{\gamma_3}_{\beta_2} - $$ - If this operation is done by a Kronecker product of unfolded tensors, - the resulting tensor has permuted indices. So, in this case the - permutation is implied by the equivalence: - $\{\{0,4\},\{1,3\},\{2\}\}$. This results in a permutation which maps - indices $(0,1,2,3,4)\mapsto(0,2,4,3,1)$. - The other application of |Permutation| class is to permute indices - with the same permutation as done during sorting. + [g_y³]_γ₁γ₂γ₃ [g_yu]^γ₁_α₁β₃ [g_yu]^γ₂_α₂β₁ [g_u]^γ₃_β₂ + + If this operation is done by a Kronecker product of unfolded tensors, the + resulting tensor has permuted indices. So, in this case the permutation is + implied by the equivalence: { {0,4}, {1,3}, {2} }. This results in a + permutation which maps indices (0,1,2,3,4)↦(0,2,4,3,1). + + The other application of Permutation class is to permute indices with the + same permutation as done during sorting. Here we only define an abstraction for the permutation defined by an - equivalence. Its basic operation is to apply the permutation to the - integer sequence. The application is right (or inner), in sense that - it works on indices of the sequence not items of the sequence. More - formally $s\circ m \not=m\circ s$. In here, the application of the - permutation defined by map $m$ is $s\circ m$. + equivalence. Its basic operation is to apply the permutation to the integer + sequence. The application is right (or inner), in sense that it works on + indices of the sequence not items of the sequence. More formally s∘m ≠ m∘s. + In here, the application of the permutation defined by map m is s∘m. - Also, we need |PermutationSet| class which contains all permutations - of $n$ element set, and a bundle of permutations |PermutationBundle| - which contains all permutation sets up to a given number. */ + Also, we need PermutationSet class which contains all permutations + of an n-element set, and a bundle of permutations PermutationBundle + which contains all permutation sets up to a given number. + */ #ifndef PERMUTATION_H #define PERMUTATION_H @@ -38,17 +35,17 @@ #include -/* The permutation object will have a map, which defines mapping of - indices $(0,1,\ldots,n-1)\mapsto(m_0,m_1,\ldots, m_{n-1})$. The map is - the sequence $(m_0,m_1,\ldots, m_{n-1}$. When the permutation with the - map $m$ is applied on sequence $s$, it permutes its indices: - $s\circ\hbox{id}\mapsto s\circ m$. +/* The permutation object will have a map, which defines mapping of indices + (0,1,…,n-1)↦(m₀,m₁,…, mₙ₋₁). The map is the sequence (m₀,m₁,…, mₙ₋₁). When + the permutation with the map m is applied on sequence s, it permutes its + indices: s∘id↦s∘m. - So we have one constructor from equivalence, then a method |apply|, - and finally a method |tailIdentity| which returns a number of trailing + So we have one constructor from equivalence, then a method apply(), + and finally a method tailIdentity() which returns a number of trailing indices which yield identity. Also we have a constructor calculating map, which corresponds to permutation in sort. This is, we want - $(\hbox{sorted }s)\circ m = s$. */ + (sorted s)∘m = s. + */ class Permutation { @@ -118,17 +115,15 @@ protected: void computeSortingMap(const IntSequence &s); }; -/* The |PermutationSet| maintains an array of of all permutations. The - default constructor constructs one element permutation set of one - element sets. The second constructor constructs a new permutation set - over $n$ from all permutations over $n-1$. The parameter $n$ need not - to be provided, but it serves to distinguish the constructor from copy - constructor. +/* The PermutationSet maintains an array of of all permutations. The default + constructor constructs one element permutation set of one element sets. The + second constructor constructs a new permutation set over n from all + permutations over n-1. The parameter n needs not to be provided, but it + serves to distinguish the constructor from copy constructor. - The method |getPreserving| returns a factor subgroup of permutations, - which are invariants with respect to the given sequence. This are all - permutations $p$ yielding $p\circ s = s$, where $s$ is the given - sequence. */ + The method getPreserving() returns a factor subgroup of permutations, which + are invariants with respect to the given sequence. This are all permutations + p yielding p∘s = s, where s is the given sequence. */ class PermutationSet { diff --git a/dynare++/tl/cc/ps_tensor.cc b/dynare++/tl/cc/ps_tensor.cc index 7905d6954..35975ecc1 100644 --- a/dynare++/tl/cc/ps_tensor.cc +++ b/dynare++/tl/cc/ps_tensor.cc @@ -8,11 +8,10 @@ #include -/* Here we decide, what method for filling a slice in slicing - constructor to use. A few experiments suggest, that if the tensor is - more than 8\% filled, the first method (|fillFromSparseOne|) is - better. For fill factors less than 1\%, the second can be 3 times - quicker. */ +/* Here we decide, what method for filling a slice in slicing constructor to + use. A few experiments suggest, that if the tensor is more than 8% filled, + the first method (fillFromSparseOne()) is better. For fill factors less than + 1%, the second can be 3 times quicker. */ UPSTensor::fill_method UPSTensor::decideFillMethod(const FSSparseTensor &t) @@ -100,17 +99,17 @@ UPSTensor::addTo(FGSTensor &out) const iterator, which iterates over tensor at some higher levels. So we simulate it by the following code. - First we set |cols| to the length of the data chunk and |off| to its - dimension. Then we need a front part of |nvmax| of |out|, which is - |nvmax_part|. Our iterator here is an integer sequence |outrun| with - full length, and |outrun_part| its front part. The |outrun| is - initialized to zeros. In each step we need to increment |outrun| - |cols|-times, this is done by incrementing its prefix |outrun_part|. + First we set ‘cols’ to the length of the data chunk and ‘off’ to its + dimension. Then we need a front part of ‘nvmax’ of ‘out’, which is + ‘nvmax_part’. Our iterator here is an integer sequence ‘outrun’ with + full length, and ‘outrun_part’ its front part. The ‘outrun’ is + initialized to zeros. In each step we need to increment ‘outrun’ + ‘cols’-times, this is done by incrementing its prefix ‘outrun_part’. - So we loop over all |cols|wide partitions of |out|, permute |outrun| - to obtain |perrun| to obtain column of this matrix. (note that the - trailing part of |perrun| is the same as of |outrun|. Then we - construct submatrices, add them, and increment |outrun|. */ + So we loop over all ‘cols’-wide partitions of ‘out’, permute ‘outrun’ + to obtain ‘perrun’ to obtain column of this matrix. (note that the + trailing part of ‘perrun’ is the same as of ‘outrun’. Then we + construct submatrices, add them, and increment ‘outrun’. */ void UPSTensor::addTo(UGSTensor &out) const @@ -124,7 +123,7 @@ UPSTensor::addTo(UGSTensor &out) const IntSequence nvmax_part(out.getDims().getNVX(), 0, out.dimen()-off); for (int out_col = 0; out_col < out.ncols(); out_col += cols) { - // permute |outrun| + // permute ‘outrun’ IntSequence perrun(out.dimen()); tdims.getPer().apply(outrun, perrun); index from(*this, perrun); @@ -133,12 +132,12 @@ UPSTensor::addTo(UGSTensor &out) const TwoDMatrix subout(out, out_col, cols); // add subout.add(1, subfrom); - // increment |outrun| by cols + // increment ‘outrun’ by cols UTensor::increment(outrun_part, nvmax_part); } } -/* This returns a product of all items in |nvmax| which make up the +/* This returns a product of all items in ‘nvmax’ which make up the trailing identity part. */ int @@ -147,7 +146,7 @@ UPSTensor::tailIdentitySize() const return tdims.getNVX().mult(dimen()-tdims.tailIdentity(), dimen()); } -/* This fill method is pretty dumb. We go through all columns in |this| +/* This fill method is pretty dumb. We go through all columns in ‘this’ tensor, translate coordinates to sparse tensor, sort them and find an item in the sparse tensor. There are many not successful lookups for really sparse tensor, that is why the second method works better for @@ -181,20 +180,18 @@ UPSTensor::fillFromSparseOne(const FSSparseTensor &t, const IntSequence &ss, } } -/* This is the second way of filling the slice. For instance, let the - slice correspond to partitions $abac$. In here we first calculate - lower and upper bounds for index of the sparse tensor for the - slice. These are |lb_srt| and |ub_srt| respectively. They corresponds - to ordering $aabc$. Then we go through that interval, and select items - which are really between the bounds. Then we take the index, subtract - the lower bound to get it to coordinates of the slice. We get - something like $(i_a,j_a,k_b,l_c)$. Then we apply the inverse - permutation as of the sorting form $abac\mapsto aabc$ to get index - $(i_a,k_b,j_a,l_c)$. Recall that the slice is unfolded, so we have to - apply all permutations preserving the stack coordinates $abac$. In our - case we get list of indices $(i_a,k_b,j_a,l_c)$ and - $(j_a,k_b,i_a,l_c)$. For all we copy the item of the sparse tensor to - the appropriate column. */ +/* This is the second way of filling the slice. For instance, let the slice + correspond to partitions “abac”. In here we first calculate lower and upper + bounds for index of the sparse tensor for the slice. These are ‘lb_srt’ and + ‘ub_srt’ respectively. They corresponds to ordering “aabc”. Then we go + through that interval, and select items which are really between the bounds. + Then we take the index, subtract the lower bound to get it to coordinates of + the slice. We get something like (i_a,j_a,k_b,l_c). Then we apply the + inverse permutation as of the sorting form abac ↦ aabc to get index + (i_a,k_b,j_a,l_c). Recall that the slice is unfolded, so we have to apply + all permutations preserving the stack coordinates “abac”. In our case we get + list of indices (i_a,k_b,j_a,l_c) and (j_a,k_b,i_a,l_c). For all we copy the + item of the sparse tensor to the appropriate column. */ void UPSTensor::fillFromSparseTwo(const FSSparseTensor &t, const IntSequence &ss, @@ -242,7 +239,7 @@ UPSTensor::fillFromSparseTwo(const FSSparseTensor &t, const IntSequence &ss, } /* Here we calculate the maximum offsets in each folded dimension - (dimension sizes, hence |ds|). */ + (dimension sizes, hence ‘ds’). */ void PerTensorDimens2::setDimensionSizes() @@ -255,12 +252,12 @@ PerTensorDimens2::setDimensionSizes() } } -/* If there are two folded dimensions, the offset in such a dimension - is offset of the second plus offset of the first times the maximum - offset of the second. If there are $n+1$ dimensions, the offset is a - sum of offsets of the last dimension plus the offset in the first $n$ - dimensions multiplied by the maximum offset of the last - dimension. This is exactly what the following code does. */ +/* If there are two folded dimensions, the offset in such a dimension is offset + of the second plus offset of the first times the maximum offset of the + second. If there are n+1 dimensions, the offset is a sum of offsets of the + last dimension plus the offset in the first n dimensions multiplied by the + maximum offset of the last dimension. This is exactly what the following + code does. */ int PerTensorDimens2::calcOffset(const IntSequence &coor) const @@ -295,9 +292,8 @@ PerTensorDimens2::print() const } /* Here we increment the given integer sequence. It corresponds to - |UTensor::increment| of the whole sequence, and then partial - monotonizing of the subsequences with respect to the - symmetries of each dimension. */ + UTensor::increment() of the whole sequence, and then partial monotonizing of + the subsequences with respect to the symmetries of each dimension. */ void FPSTensor::increment(IntSequence &v) const @@ -326,7 +322,7 @@ FPSTensor::unfold() const TL_RAISE("Unfolding of FPSTensor not implemented"); } -/* We only call |calcOffset| of the |PerTensorDimens2|. */ +/* We only call calcOffset() on the PerTensorDimens2. */ int FPSTensor::getOffset(const IntSequence &v) const @@ -334,8 +330,8 @@ FPSTensor::getOffset(const IntSequence &v) const return tdims.calcOffset(v); } -/* Here we add the tensor to |out|. We go through all columns of the - |out|, apply the permutation to get index in the tensor, and add the +/* Here we add the tensor to ‘out’. We go through all columns of the + ‘out’, apply the permutation to get index in the tensor, and add the column. Note that if the permutation is identity, then the dimensions of the tensors might not be the same (since this tensor is partially folded). */ @@ -353,11 +349,11 @@ FPSTensor::addTo(FGSTensor &out) const } /* Here is the constructor which multiplies the Kronecker product with - the general symmetry sparse tensor |GSSparseTensor|. The main idea is + the general symmetry sparse tensor GSSparseTensor. The main idea is to go through items in the sparse tensor (each item selects rows in the matrices form the Kornecker product), then to Kronecker multiply the rows and multiply with the item, and to add the resulting row to - the appropriate row of the resulting |FPSTensor|. + the appropriate row of the resulting FPSTensor. The realization of this idea is a bit more complicated since we have to go through all items, and each item must be added as many times as @@ -365,12 +361,12 @@ FPSTensor::addTo(FGSTensor &out) const order of rows in their Kronecker product. So, we through all unfolded indices in a tensor with the same - dimensions as the |GSSparseTensor| (sparse slice). For each such index + dimensions as the GSSparseTensor (sparse slice). For each such index we calculate its folded version (corresponds to ordering of subsequences within symmetries), we test if there is an item in the sparse slice with such coordinates, and if there is, we construct the Kronecker product of the rows, and go through all of items with the - coordinates, and add to appropriate rows of |this| tensor. */ + coordinates, and add to appropriate rows of ‘this’ tensor. */ FPSTensor::FPSTensor(const TensorDimens &td, const Equivalence &e, const Permutation &p, const GSSparseTensor &a, const KronProdAll &kp) diff --git a/dynare++/tl/cc/ps_tensor.hh b/dynare++/tl/cc/ps_tensor.hh index 9d4472dbf..45d9c73b4 100644 --- a/dynare++/tl/cc/ps_tensor.hh +++ b/dynare++/tl/cc/ps_tensor.hh @@ -2,43 +2,39 @@ // Even more general symmetry tensor. -/* Here we define an abstraction for a tensor, which has a general - symmetry, but the symmetry is not of what is modelled by - |Symmetry|. This kind of tensor comes to existence when we evaluate - something like: - $$\left[B_{y^2u^3}\right]_{\alpha_1\alpha_2\beta_1\beta_2\beta_3}= - \cdots+\left[g_{y^3}\right]_{\gamma_1\gamma_2\gamma_3} - \left[g_{yu}\right]^{\gamma_1}_{\alpha_1\beta_3} - \left[g_{yu}\right]^{\gamma_2}_{\alpha_2\beta_1} - \left[g_u\right]^{\gamma_3}_{\beta_2}+\cdots - $$ +/* Here we define an abstraction for a tensor, which has a general symmetry, + but the symmetry is not of what is modelled by Symmetry. This kind of tensor + comes to existence when we evaluate something like: + + [B_y²u³]_α₁α₂β₁β₂β₃ = … + [g_y³]_γ₁γ₂γ₃ [g_yu]^γ₁_α₁β₃ [g_yu]^γ₂_α₂β₁ + [g_u]^γ₃_β₂ + … + If the tensors are unfolded, we obtain a tensor - $$g_{y^3}\cdot\left(g_{yu}\otimes g_{yu}\otimes g_{u}\right)$$ - Obviously, this tensor can have a symmetry not compatible with - ordering $\alpha_1\alpha_2\beta_1\beta_2\beta_3$, (in other words, not - compatible with symmetry $y^2u^3$). In fact, the indices are permuted. + g_y³·(g_yu ⊗ g_yu ⊗ g_u) - This kind of tensor must be added to $\left[B_{y^2u^3}\right]$. Its - dimensions are the same as of $\left[B_{y^2u^3}\right]$, but some - coordinates are permuted. The addition is the only action we need to - do with the tensor. + Obviously, this tensor can have a symmetry not compatible with ordering + α₁α₂β₁β₂β₃, (in other words, not compatible with symmetry y²u³). In fact, + the indices are permuted. - Another application where this permuted symmetry tensor appears is a - slice of a fully symmetric tensor. If the symmetric dimension of the - tensor is partitioned to continuous parts, and we are interested only - in data with a given symmetry (permuted) of the partitions, then we - have the permuted symmetry tensor. For instance, if $x$ is partitioned - $x=[a,b,c,d]$, and having tensor $\left[f_{x^3}\right]$, one can d a - slice (subtensor) $\left[f_{aca}\right]$. The data of this tensor are - permuted of $\left[f_{a^c}\right]$. + This kind of tensor must be added to [B_y²u³]. Its dimensions are the same + as of [B_y²u³], but some coordinates are permuted. The addition is the only + action we need to do with the tensor. + + Another application where this permuted symmetry tensor appears is a slice + of a fully symmetric tensor. If the symmetric dimension of the tensor is + partitioned to continuous parts, and we are interested only in data with a + given symmetry (permuted) of the partitions, then we have the permuted + symmetry tensor. For instance, if x is partitioned x=(a,b,c,d), and having + tensor [f_x³], one can define a slice (subtensor) [f_aca]. The data of this + tensor are permuted of [f_a²c]. Here we also define the folded version of permuted symmetry tensor. It has permuted symmetry and is partially folded. One can imagine it as a product of a few dimensions, each of them is folded and having a few variables. The underlying variables are permuted. The product of such - dimensions is described by |PerTensorDimens2|. The tensor holding the - underlying data is |FPSTensor|. */ + dimensions is described by PerTensorDimens2. The tensor holding the + underlying data is FPSTensor. */ #ifndef PS_TENSOR_H #define PS_TENSOR_H @@ -50,25 +46,24 @@ #include "kron_prod.hh" #include "sparse_tensor.hh" -/* Here we declare a class describing dimensions of permuted symmetry - tensor. It inherits from |TensorDimens| and adds a permutation which - permutes |nvmax|. It has two constructors, each corresponds to a - context where the tensor appears. +/* Here we declare a class describing dimensions of permuted symmetry tensor. + It inherits from TensorDimens and adds a permutation which permutes ‘nvmax’. + It has two constructors, each corresponds to a context where the tensor + appears. The first constructor calculates the permutation from a given equivalence. - The second constructor corresponds to dimensions of a slice. Let us - take $\left[f_{aca}\right]$ as an example. First it calculates - |TensorDimens| of $\left[f_{a^c}\right]$, then it calculates a - permutation corresponding to ordering of $aca$ to $a^2c$, and applies - this permutation on the dimensions as the first constructor. The - constructor takes only stack sizes (lengths of $a$, $b$, $c$, and - $d$), and coordinates of picked partitions. + The second constructor corresponds to dimensions of a slice. Let us take + [f_aca] as an example. First it calculates TensorDimens of [f_a²c], then it + calculates a permutation corresponding to ordering of “aca” to “a²c”, and + applies this permutation on the dimensions as the first constructor. The + constructor takes only stack sizes (lengths of a, b, c, and d), and + coordinates of picked partitions. - Note that inherited methods |calcUnfoldColumns| and |calcFoldColumns| - work, since number of columns is independent on the permutation, and - |calcFoldColumns| does not use changed |nvmax|, it uses |nvs|, so it - is OK. */ + Note that inherited methods calcUnfoldColumns() and calcFoldColumns() work, + since number of columns is independent on the permutation, and + calcFoldColumns() does not use changed ‘nvmax’, it uses ‘nvs’, so it is + OK. */ class PerTensorDimens : public TensorDimens { @@ -119,25 +114,25 @@ public: }; /* Here we declare the permuted symmetry unfolded tensor. It has - |PerTensorDimens| as a member. It inherits from |UTensor| which - requires to implement |fold| method. There is no folded counterpart, - so in our implementation we raise unconditional exception, and return - some dummy object (just to make it compilable without warnings). + PerTensorDimens as a member. It inherits from UTensor which requires to + implement fold() method. There is no folded counterpart, so in our + implementation we raise unconditional exception, and return some dummy + object (just to make it compilable without warnings). The class has two sorts of constructors corresponding to a context where it - appears. The first constructs object from a given matrix, and - Kronecker product. Within the constructor, all the calculations are - performed. Also we need to define dimensions, these are the same of - the resulting matrix (in our example $\left[B_{y^2u^3}\right]$) but - permuted. The permutation is done in |PerTensorDimens| constructor. + appears. The first constructs object from a given matrix, and Kronecker + product. Within the constructor, all the calculations are performed. Also we + need to define dimensions, these are the same of the resulting matrix (in + our example [B_y²u³]) but permuted. The permutation is done in + PerTensorDimens constructor. The second type of constructor is slicing. It makes a slice from - |FSSparseTensor|. The slice is given by stack sizes, and coordinates of + FSSparseTensor. The slice is given by stack sizes, and coordinates of picked stacks. There are two algorithms for filling a slice of a sparse tensor. The - first |fillFromSparseOne| works well for more dense tensors, the - second |fillFromSparseTwo| is better for very sparse tensors. We + first fillFromSparseOne() works well for more dense tensors, the + second fillFromSparseTwo() is better for very sparse tensors. We provide a static method, which decides what of the two algorithms is better. */ @@ -145,18 +140,12 @@ class UPSTensor : public UTensor { const PerTensorDimens tdims; public: - // |UPSTensor| constructors from Kronecker product - /* Here we have four constructors making an |UPSTensor| from a product - of matrix and Kronecker product. The first constructs the tensor from - equivalence classes of the given equivalence in an order given by the - equivalence. The second does the same but with optimized - |KronProdAllOptim|, which has a different order of matrices than given - by the classes in the equivalence. This permutation is projected to - the permutation of the |UPSTensor|. The third, is the same as the - first, but the classes of the equivalence are permuted by the given - permutation. Finally, the fourth is the most general combination. It - allows for a permutation of equivalence classes, and for optimized - |KronProdAllOptim|, which permutes the permuted equivalence classes. */ + // UPSTensor constructors from Kronecker product + /* Here we have four constructors making an UPSTensor from a product + of matrix and Kronecker product. */ + + /* Constructs the tensor from equivalence classes of the given equivalence in + an order given by the equivalence */ UPSTensor(TensorDimens td, const Equivalence &e, const ConstTwoDMatrix &a, const KronProdAll &kp) : UTensor(indor::along_col, PerTensorDimens(td, e).getNVX(), @@ -165,6 +154,10 @@ public: { kp.mult(a, *this); } + + /* Same as the previous one but with optimized KronProdAllOptim, which has a + different order of matrices than given by the classes in the equivalence. + This permutation is projected to the permutation of the UPSTensor. */ UPSTensor(TensorDimens td, const Equivalence &e, const ConstTwoDMatrix &a, const KronProdAllOptim &kp) : UTensor(indor::along_col, PerTensorDimens(td, Permutation(e, kp.getPer())).getNVX(), @@ -173,6 +166,9 @@ public: { kp.mult(a, *this); } + + /* Same as the first constructor, but the classes of the equivalence are + permuted by the given permutation. */ UPSTensor(TensorDimens td, const Equivalence &e, const Permutation &p, const ConstTwoDMatrix &a, const KronProdAll &kp) : UTensor(indor::along_col, PerTensorDimens(td, Permutation(e, p)).getNVX(), @@ -181,6 +177,10 @@ public: { kp.mult(a, *this); } + + /* Most general constructor. It allows for a permutation of equivalence + classes, and for optimized KronProdAllOptim, which permutes the permuted + equivalence classes. */ UPSTensor(TensorDimens td, const Equivalence &e, const Permutation &p, const ConstTwoDMatrix &a, const KronProdAllOptim &kp) : UTensor(indor::along_col, PerTensorDimens(td, Permutation(e, Permutation(p, kp.getPer()))).getNVX(), @@ -189,6 +189,7 @@ public: { kp.mult(a, *this); } + UPSTensor(const FSSparseTensor &t, const IntSequence &ss, const IntSequence &coor, PerTensorDimens ptd); UPSTensor(const UPSTensor &) = default; @@ -212,17 +213,16 @@ private: const IntSequence &coor); }; -/* Here we define an abstraction for the tensor dimension with the - symmetry like $xuv\vert uv\vert xu\vert y\vert y\vert x\vert x\vert - y$. These symmetries come as induces symmetries of equivalence and - some outer symmetry. Thus the underlying variables are permuted. One - can imagine the dimensions as an unfolded product of dimensions which - consist of folded products of variables. +/* Here we define an abstraction for the tensor dimension with the symmetry + like xuv|uv|xu|y|y|x|x|y. These symmetries come as induces symmetries of + equivalence and some outer symmetry. Thus the underlying variables are + permuted. One can imagine the dimensions as an unfolded product of + dimensions which consist of folded products of variables. - We inherit from |PerTensorDimens| since we need the permutation - implied by the equivalence. The new member are the induced symmetries - (symmetries of each folded dimensions) and |ds| which are sizes of the - dimensions. The number of folded dimensions is return by |numSyms|. + We inherit from PerTensorDimens since we need the permutation implied by the + equivalence. The new member are the induced symmetries (symmetries of each + folded dimensions) and ‘ds’ which are sizes of the dimensions. The number of + folded dimensions is return by ‘numSyms’. The object is constructed from outer tensor dimensions and from equivalence with optionally permuted classes. */ @@ -268,44 +268,38 @@ protected: void setDimensionSizes(); }; -/* Here we define an abstraction of the permuted symmetry folded - tensor. It is needed in context of the Faà Di Bruno formula for folded - stack container multiplied with container of dense folded tensors, or - multiplied by one full symmetry sparse tensor. +/* Here we define an abstraction of the permuted symmetry folded tensor. It is + needed in context of the Faà Di Bruno formula for folded stack container + multiplied with container of dense folded tensors, or multiplied by one full + symmetry sparse tensor. For example, if we perform the Faà Di Bruno for $F=f(z)$, where - $z=[g(x,y,u,v), h(x,y,u), x, y]^T$, we get for one concrete - equivalence: - $$ - \left[F_{x^4y^3u^3v^2}\right]=\ldots+ - \left[f_{g^2h^2x^2y}\right]\left( - [g]_{xv}\otimes[g]_{u^2v}\otimes - [h]_{xu}\otimes[h]_{y^2}\otimes - \left[\vphantom{\sum}[I]_x\otimes[I]_x\right]\otimes - \left[\vphantom{\sum}[I]_y\right] - \right) - +\ldots - $$ + ⎛g(x,y,u,v)⎞ + ⎢ h(x,y,u) ⎥ + z=⎢ x ⎥ + ⎝ y ⎠ + we get for one concrete equivalence: - The class |FPSTensor| represents the tensor at the right. Its - dimension corresponds to a product of 7 dimensions with the following - symmetries: $xv\vert u^v\vert xu\vert y^2\vert x\vert x\vert y$. Such - the dimension is described by |PerTensorDimens2|. + [F_x⁴y³u³v²] = … + [f_g²h²x²y]·([g]_xv ⊗ [g]_u²v ⊗ [h]_xu ⊗ [h]_y² ⊗ + [I]_x ⊗ [I]_x ⊗ [I]_y) + + … - The tensor is constructed in a context of stack container - multiplication, so, it is constructed from dimensions |td| (dimensions - of the output tensor), stack product |sp| (implied symmetries picking - tensors from a stack container, here it is $z$), then a sorted integer - sequence of the picked stacks of the stack product (it is always - sorted, here it is $(0,0,1,1,2,2,3)$), then the tensor - $\left[f_{g^2h^2x^2y}\right]$ (its symmetry must be the same as - symmetry given by the |istacks|), and finally from the equivalence - with permuted classes. + The class FPSTensor represents the tensor on the right. Its dimension + corresponds to a product of 7 dimensions with the following symmetries: + xv|u²v|xu|y²|x|x|y. Such a dimension is described by PerTensorDimens2. - We implement |increment| and |getOffset| methods, |decrement| and - |unfold| raise an exception. Also, we implement |addTo| method, which - adds the tensor data (partially unfolded) to folded general symmetry - tensor. */ + The tensor is constructed in a context of stack container multiplication, + so, it is constructed from dimensions ‘td’ (dimensions of the output + tensor), stack product ‘sp’ (implied symmetries picking tensors from a stack + container, here it is z), then a sorted integer sequence of the picked + stacks of the stack product (it is always sorted, here it is + (0,0,1,1,2,2,3)), then the tensor [f_g²h²x²y] (its symmetry must be the same + as symmetry given by the ‘istacks’), and finally from the equivalence with + permuted classes. + + We implement increment() and getOffset() methods, decrement() and unfold() + raise an exception. Also, we implement addTo() method, which adds the tensor + data (partially unfolded) to folded general symmetry tensor. */ template class StackProduct; @@ -314,9 +308,9 @@ class FPSTensor : public FTensor { const PerTensorDimens2 tdims; public: - /* As for |UPSTensor|, we provide four constructors allowing for + /* As for UPSTensor, we provide four constructors allowing for combinations of permuting equivalence classes, and optimization of - |KronProdAllOptim|. These constructors multiply with dense general + KronProdAllOptim. These constructors multiply with dense general symmetry tensor (coming from the dense container, or as a dense slice of the full symmetry sparse tensor). In addition to these 4 constructors, we have one constructor multiplying with general