1 /*
2  * Copyright © 2004 Ondra Kamenik
3  * Copyright © 2019 Dynare Team
4  *
5  * This file is part of Dynare.
6  *
7  * Dynare is free software: you can redistribute it and/or modify
8  * it under the terms of the GNU General Public License as published by
9  * the Free Software Foundation, either version 3 of the License, or
10  * (at your option) any later version.
11  *
12  * Dynare is distributed in the hope that it will be useful,
13  * but WITHOUT ANY WARRANTY; without even the implied warranty of
14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
15  * GNU General Public License for more details.
16  *
17  * You should have received a copy of the GNU General Public License
18  * along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
19  */
20 
21 // Row-wise full symmetry tensor.
22 
23 /* Here we define classes for full symmetry tensors with the multidimensional
24    index identified with rows. The primary usage is for storage of data coming
25    from (or from a sum of)
26       ₗ
27       ∏  [g_{s^|cₘ|}]_cₘ(α)^γₘ
28      ᵐ⁼¹
29    where α comes from a multidimensional index that goes through some set S,
30    and cₘ is some equivalence class. So we model a tensor of the form:
31 
32      ⎡ ₗ                       ⎤
33      ⎢ ∏  [g_{s^|cₘ|}]_cₘ(α)^γₘ⎥
34      ⎣ᵐ⁼¹                      ⎦S^γ₁…γₗ
35 
36    Since all γ₁…γₗ correspond to the same variable, the tensor is fully
37    symmetric. The set of indices S cannot be very large and sometimes it is
38    only one element. This case is handled in a special subclass.
39 
40    We provide both folded and unfolded versions. Their logic is perfectly the
41    same as in UFSTensor and FFSTensor with two exceptions. One has been already
42    mentioned, the multidimensional index is along the rows. The second are
43    conversions between the two types. Since this kind of tensor is used to
44    multiply (from the right) a tensor whose multidimensional index is
45    identified with columns, we will need a different way of a conversion. If
46    the multiplication of two folded tensors is to be equivalent with
47    multiplication of two unfolded, the folding of the right tensor must sum all
48    equivalent elements since they are multiplied with the same number from the
49    folded tensor. (Equivalent here means all elements of unfolded tensor
50    corresponding to one element in folded tensor.) For this reason, it is
51    necessary to calculate a column number from the given sequence, so we
52    implement getOffset(). Process of unfolding is not used, so we implemented
53    it so that unfolding and then folding a tensor would yield the same data. */
54 
55 #ifndef RFS_TENSOR_H
56 #define RFS_TENSOR_H
57 
58 #include "tensor.hh"
59 #include "fs_tensor.hh"
60 #include "symmetry.hh"
61 
62 /* This is straightforward and very similar to UFSTensor. */
63 
64 class FRTensor;
65 class URTensor : public UTensor
66 {
67   int nv;
68 public:
URTensor(int c,int nvar,int d)69   URTensor(int c, int nvar, int d)
70     : UTensor(indor::along_row, IntSequence(d, nvar),
71               UFSTensor::calcMaxOffset(nvar, d), c, d), nv(nvar)
72   {
73   }
74   URTensor(const URTensor &) = default;
75   URTensor(URTensor &&) = default;
76   explicit URTensor(const FRTensor &ft);
77 
78   ~URTensor() override = default;
79 
80   void increment(IntSequence &v) const override;
81   void decrement(IntSequence &v) const override;
82   std::unique_ptr<FTensor> fold() const override;
83 
84   int getOffset(const IntSequence &v) const override;
85   int
nvar() const86   nvar() const
87   {
88     return nv;
89   }
90   Symmetry
getSym() const91   getSym() const
92   {
93     return Symmetry{dimen()};
94   }
95 };
96 
97 /* This is straightforward and very similar to FFSTensor. */
98 
99 class FRTensor : public FTensor
100 {
101   int nv;
102 public:
FRTensor(int c,int nvar,int d)103   FRTensor(int c, int nvar, int d)
104     : FTensor(indor::along_row, IntSequence(d, nvar),
105               FFSTensor::calcMaxOffset(nvar, d), c, d), nv(nvar)
106   {
107   }
108   FRTensor(const FRTensor &) = default;
109   FRTensor(FRTensor &&) = default;
110   explicit FRTensor(const URTensor &ut);
111 
112   ~FRTensor() override = default;
113 
114   void increment(IntSequence &v) const override;
115   void decrement(IntSequence &v) const override;
116   std::unique_ptr<UTensor> unfold() const override;
117 
118   int
nvar() const119   nvar() const
120   {
121     return nv;
122   }
123   int
getOffset(const IntSequence & v) const124   getOffset(const IntSequence &v) const override
125   {
126     return FTensor::getOffset(v, nv);
127   }
128   Symmetry
getSym() const129   getSym() const
130   {
131     return Symmetry{dimen()};
132   }
133 };
134 
135 /* The following class represents specialization of URTensor coming from
136    Kronecker multiplication of a few vectors. So the resulting row-oriented
137    tensor has one column. We provide two constructors, one constructs the
138    tensor from a few vectors stored as std::vector<ConstVector>. The second
139    makes the Kronecker power of one given vector. */
140 
141 class URSingleTensor : public URTensor
142 {
143 public:
URSingleTensor(int nvar,int d)144   URSingleTensor(int nvar, int d)
145     : URTensor(1, nvar, d)
146   {
147   }
148   URSingleTensor(const std::vector<ConstVector> &cols);
149   URSingleTensor(const ConstVector &v, int d);
150   URSingleTensor(const URSingleTensor &) = default;
151   URSingleTensor(URSingleTensor &&) = default;
152   ~URSingleTensor() override = default;
153   std::unique_ptr<FTensor> fold() const override;
154 };
155 
156 /* This class represents one column row-oriented tensor. The only way to
157    construct it is from URSingleTensor or from scratch. The folding algorithm
158    is the same as folding of general URTensor. Only its implementation is
159    different, since we do not copy rows, but only elements. */
160 
161 class FRSingleTensor : public FRTensor
162 {
163 public:
FRSingleTensor(int nvar,int d)164   FRSingleTensor(int nvar, int d)
165     : FRTensor(1, nvar, d)
166   {
167   }
168   explicit FRSingleTensor(const URSingleTensor &ut);
169   FRSingleTensor(const FRSingleTensor &) = default;
170   FRSingleTensor(FRSingleTensor &&) = default;
171   ~FRSingleTensor() override = default;
172 };
173 
174 #endif
175