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 // Even more general symmetry tensor. 22 23 /* Here we define an abstraction for a tensor, which has a general symmetry, 24 but the symmetry is not of what is modelled by Symmetry. This kind of tensor 25 comes to existence when we evaluate something like: 26 27 [B_y²u³]_α₁α₂β₁β₂β₃ = … + [g_y³]_γ₁γ₂γ₃ [g_yu]^γ₁_α₁β₃ [g_yu]^γ₂_α₂β₁ 28 [g_u]^γ₃_β₂ + … 29 30 If the tensors are unfolded, we obtain a tensor 31 32 g_y³·(g_yu ⊗ g_yu ⊗ g_u) 33 34 Obviously, this tensor can have a symmetry not compatible with ordering 35 α₁α₂β₁β₂β₃, (in other words, not compatible with symmetry y²u³). In fact, 36 the indices are permuted. 37 38 This kind of tensor must be added to [B_y²u³]. Its dimensions are the same 39 as of [B_y²u³], but some coordinates are permuted. The addition is the only 40 action we need to do with the tensor. 41 42 Another application where this permuted symmetry tensor appears is a slice 43 of a fully symmetric tensor. If the symmetric dimension of the tensor is 44 partitioned to continuous parts, and we are interested only in data with a 45 given symmetry (permuted) of the partitions, then we have the permuted 46 symmetry tensor. For instance, if x is partitioned x=(a,b,c,d), and having 47 tensor [f_x³], one can define a slice (subtensor) [f_aca]. The data of this 48 tensor are permuted of [f_a²c]. 49 50 Here we also define the folded version of permuted symmetry tensor. It 51 has permuted symmetry and is partially folded. One can imagine it as a 52 product of a few dimensions, each of them is folded and having a few 53 variables. The underlying variables are permuted. The product of such 54 dimensions is described by PerTensorDimens2. The tensor holding the 55 underlying data is FPSTensor. */ 56 57 #ifndef PS_TENSOR_H 58 #define PS_TENSOR_H 59 60 #include "tensor.hh" 61 #include "gs_tensor.hh" 62 #include "equivalence.hh" 63 #include "permutation.hh" 64 #include "kron_prod.hh" 65 #include "sparse_tensor.hh" 66 67 /* Here we declare a class describing dimensions of permuted symmetry tensor. 68 It inherits from TensorDimens and adds a permutation which permutes ‘nvmax’. 69 It has two constructors, each corresponds to a context where the tensor 70 appears. 71 72 The first constructor calculates the permutation from a given equivalence. 73 74 The second constructor corresponds to dimensions of a slice. Let us take 75 [f_aca] as an example. First it calculates TensorDimens of [f_a²c], then it 76 calculates a permutation corresponding to ordering of “aca” to “a²c”, and 77 applies this permutation on the dimensions as the first constructor. The 78 constructor takes only stack sizes (lengths of a, b, c, and d), and 79 coordinates of picked partitions. 80 81 Note that inherited methods calcUnfoldColumns() and calcFoldColumns() work, 82 since number of columns is independent on the permutation, and 83 calcFoldColumns() does not use changed ‘nvmax’, it uses ‘nvs’, so it is 84 OK. */ 85 86 class PerTensorDimens : public TensorDimens 87 { 88 private: 89 static IntSequence sortIntSequence(IntSequence s)90 sortIntSequence(IntSequence s) 91 { 92 s.sort(); 93 return s; 94 } 95 protected: 96 Permutation per; 97 public: PerTensorDimens(Symmetry s,IntSequence nvars,const Equivalence & e)98 PerTensorDimens(Symmetry s, IntSequence nvars, const Equivalence &e) 99 : TensorDimens(std::move(s), std::move(nvars)), per(e) 100 { 101 per.apply(nvmax); 102 } PerTensorDimens(TensorDimens td,const Equivalence & e)103 PerTensorDimens(TensorDimens td, const Equivalence &e) 104 : TensorDimens(std::move(td)), per(e) 105 { 106 per.apply(nvmax); 107 } PerTensorDimens(TensorDimens td,Permutation p)108 PerTensorDimens(TensorDimens td, Permutation p) 109 : TensorDimens(std::move(td)), per(std::move(p)) 110 { 111 per.apply(nvmax); 112 } PerTensorDimens(IntSequence ss,IntSequence coor)113 PerTensorDimens(IntSequence ss, IntSequence coor) 114 : TensorDimens(std::move(ss), sortIntSequence(coor)), per(std::move(coor)) 115 { 116 per.apply(nvmax); 117 } 118 bool operator ==(const PerTensorDimens & td)119 operator==(const PerTensorDimens &td) 120 { 121 return TensorDimens::operator==(td) && per == td.per; 122 } 123 int tailIdentity() const124 tailIdentity() const 125 { 126 return per.tailIdentity(); 127 } 128 const Permutation & getPer() const129 getPer() const 130 { 131 return per; 132 } 133 }; 134 135 /* Here we declare the permuted symmetry unfolded tensor. It has 136 PerTensorDimens as a member. It inherits from UTensor which requires to 137 implement fold() method. There is no folded counterpart, so in our 138 implementation we raise unconditional exception, and return some dummy 139 object (just to make it compilable without warnings). 140 141 The class has two sorts of constructors corresponding to a context where it 142 appears. The first constructs object from a given matrix, and Kronecker 143 product. Within the constructor, all the calculations are performed. Also we 144 need to define dimensions, these are the same of the resulting matrix (in 145 our example [B_y²u³]) but permuted. The permutation is done in 146 PerTensorDimens constructor. 147 148 The second type of constructor is slicing. It makes a slice from 149 FSSparseTensor. The slice is given by stack sizes, and coordinates of 150 picked stacks. 151 152 There are two algorithms for filling a slice of a sparse tensor. The 153 first fillFromSparseOne() works well for more dense tensors, the 154 second fillFromSparseTwo() is better for very sparse tensors. We 155 provide a static method, which decides what of the two algorithms is 156 better. */ 157 158 class UPSTensor : public UTensor 159 { 160 const PerTensorDimens tdims; 161 public: 162 // UPSTensor constructors from Kronecker product 163 /* Here we have four constructors making an UPSTensor from a product 164 of matrix and Kronecker product. */ 165 166 /* Constructs the tensor from equivalence classes of the given equivalence in 167 an order given by the equivalence */ UPSTensor(TensorDimens td,const Equivalence & e,const ConstTwoDMatrix & a,const KronProdAll & kp)168 UPSTensor(TensorDimens td, const Equivalence &e, 169 const ConstTwoDMatrix &a, const KronProdAll &kp) 170 : UTensor(indor::along_col, PerTensorDimens(td, e).getNVX(), 171 a.nrows(), kp.ncols(), td.dimen()), 172 tdims(std::move(td), e) 173 { 174 kp.mult(a, *this); 175 } 176 177 /* Same as the previous one but with optimized KronProdAllOptim, which has a 178 different order of matrices than given by the classes in the equivalence. 179 This permutation is projected to the permutation of the UPSTensor. */ UPSTensor(TensorDimens td,const Equivalence & e,const ConstTwoDMatrix & a,const KronProdAllOptim & kp)180 UPSTensor(TensorDimens td, const Equivalence &e, 181 const ConstTwoDMatrix &a, const KronProdAllOptim &kp) 182 : UTensor(indor::along_col, PerTensorDimens(td, Permutation(e, kp.getPer())).getNVX(), 183 a.nrows(), kp.ncols(), td.dimen()), 184 tdims(std::move(td), Permutation(e, kp.getPer())) 185 { 186 kp.mult(a, *this); 187 } 188 189 /* Same as the first constructor, but the classes of the equivalence are 190 permuted by the given permutation. */ UPSTensor(TensorDimens td,const Equivalence & e,const Permutation & p,const ConstTwoDMatrix & a,const KronProdAll & kp)191 UPSTensor(TensorDimens td, const Equivalence &e, const Permutation &p, 192 const ConstTwoDMatrix &a, const KronProdAll &kp) 193 : UTensor(indor::along_col, PerTensorDimens(td, Permutation(e, p)).getNVX(), 194 a.nrows(), kp.ncols(), td.dimen()), 195 tdims(std::move(td), Permutation(e, p)) 196 { 197 kp.mult(a, *this); 198 } 199 200 /* Most general constructor. It allows for a permutation of equivalence 201 classes, and for optimized KronProdAllOptim, which permutes the permuted 202 equivalence classes. */ UPSTensor(TensorDimens td,const Equivalence & e,const Permutation & p,const ConstTwoDMatrix & a,const KronProdAllOptim & kp)203 UPSTensor(TensorDimens td, const Equivalence &e, const Permutation &p, 204 const ConstTwoDMatrix &a, const KronProdAllOptim &kp) 205 : UTensor(indor::along_col, PerTensorDimens(td, Permutation(e, Permutation(p, kp.getPer()))).getNVX(), 206 a.nrows(), kp.ncols(), td.dimen()), 207 tdims(std::move(td), Permutation(e, Permutation(p, kp.getPer()))) 208 { 209 kp.mult(a, *this); 210 } 211 212 UPSTensor(const FSSparseTensor &t, const IntSequence &ss, 213 const IntSequence &coor, PerTensorDimens ptd); 214 UPSTensor(const UPSTensor &) = default; 215 UPSTensor(UPSTensor &&) = default; 216 217 void increment(IntSequence &v) const override; 218 void decrement(IntSequence &v) const override; 219 std::unique_ptr<FTensor> fold() const override; 220 221 int getOffset(const IntSequence &v) const override; 222 void addTo(FGSTensor &out) const; 223 void addTo(UGSTensor &out) const; 224 225 enum class fill_method { first, second }; 226 static fill_method decideFillMethod(const FSSparseTensor &t); 227 private: 228 int tailIdentitySize() const; 229 void fillFromSparseOne(const FSSparseTensor &t, const IntSequence &ss, 230 const IntSequence &coor); 231 void fillFromSparseTwo(const FSSparseTensor &t, const IntSequence &ss, 232 const IntSequence &coor); 233 }; 234 235 /* Here we define an abstraction for the tensor dimension with the symmetry 236 like xuv|uv|xu|y|y|x|x|y. These symmetries come as induces symmetries of 237 equivalence and some outer symmetry. Thus the underlying variables are 238 permuted. One can imagine the dimensions as an unfolded product of 239 dimensions which consist of folded products of variables. 240 241 We inherit from PerTensorDimens since we need the permutation implied by the 242 equivalence. The new member are the induced symmetries (symmetries of each 243 folded dimensions) and ‘ds’ which are sizes of the dimensions. The number of 244 folded dimensions is return by ‘numSyms’. 245 246 The object is constructed from outer tensor dimensions and from 247 equivalence with optionally permuted classes. */ 248 249 class PerTensorDimens2 : public PerTensorDimens 250 { 251 InducedSymmetries syms; 252 IntSequence ds; 253 public: PerTensorDimens2(const TensorDimens & td,const Equivalence & e,const Permutation & p)254 PerTensorDimens2(const TensorDimens &td, const Equivalence &e, 255 const Permutation &p) 256 : PerTensorDimens(td, Permutation(e, p)), 257 syms(e, p, td.getSym()), 258 ds(syms.size()) 259 { 260 setDimensionSizes(); 261 } PerTensorDimens2(const TensorDimens & td,const Equivalence & e)262 PerTensorDimens2(const TensorDimens &td, const Equivalence &e) 263 : PerTensorDimens(td, e), 264 syms(e, td.getSym()), 265 ds(syms.size()) 266 { 267 setDimensionSizes(); 268 } 269 int numSyms() const270 numSyms() const 271 { 272 return static_cast<int>(syms.size()); 273 } 274 const Symmetry & getSym(int i) const275 getSym(int i) const 276 { 277 return syms[i]; 278 } 279 int calcMaxOffset() const280 calcMaxOffset() const 281 { 282 return ds.mult(); 283 } 284 int calcOffset(const IntSequence &coor) const; 285 void print() const; 286 protected: 287 void setDimensionSizes(); 288 }; 289 290 /* Here we define an abstraction of the permuted symmetry folded tensor. It is 291 needed in context of the Faà Di Bruno formula for folded stack container 292 multiplied with container of dense folded tensors, or multiplied by one full 293 symmetry sparse tensor. 294 295 For example, if we perform the Faà Di Bruno for $F=f(z)$, where 296 ⎛g(x,y,u,v)⎞ 297 ⎢ h(x,y,u) ⎥ 298 z=⎢ x ⎥ 299 ⎝ y ⎠ 300 we get for one concrete equivalence: 301 302 [F_x⁴y³u³v²] = … + [f_g²h²x²y]·([g]_xv ⊗ [g]_u²v ⊗ [h]_xu ⊗ [h]_y² ⊗ 303 [I]_x ⊗ [I]_x ⊗ [I]_y) 304 + … 305 306 The class FPSTensor represents the tensor on the right. Its dimension 307 corresponds to a product of 7 dimensions with the following symmetries: 308 xv|u²v|xu|y²|x|x|y. Such a dimension is described by PerTensorDimens2. 309 310 The tensor is constructed in a context of stack container multiplication, 311 so, it is constructed from dimensions ‘td’ (dimensions of the output 312 tensor), stack product ‘sp’ (implied symmetries picking tensors from a stack 313 container, here it is z), then a sorted integer sequence of the picked 314 stacks of the stack product (it is always sorted, here it is 315 (0,0,1,1,2,2,3)), then the tensor [f_g²h²x²y] (its symmetry must be the same 316 as symmetry given by the ‘istacks’), and finally from the equivalence with 317 permuted classes. 318 319 We implement increment() and getOffset() methods, decrement() and unfold() 320 raise an exception. Also, we implement addTo() method, which adds the tensor 321 data (partially unfolded) to folded general symmetry tensor. */ 322 323 template<typename _Ttype> 324 class StackProduct; 325 326 class FPSTensor : public FTensor 327 { 328 const PerTensorDimens2 tdims; 329 public: 330 /* As for UPSTensor, we provide four constructors allowing for 331 combinations of permuting equivalence classes, and optimization of 332 KronProdAllOptim. These constructors multiply with dense general 333 symmetry tensor (coming from the dense container, or as a dense slice 334 of the full symmetry sparse tensor). In addition to these 4 335 constructors, we have one constructor multiplying with general 336 symmetry sparse tensor (coming as a sparse slice of the full symmetry 337 sparse tensor). */ FPSTensor(const TensorDimens & td,const Equivalence & e,const ConstTwoDMatrix & a,const KronProdAll & kp)338 FPSTensor(const TensorDimens &td, const Equivalence &e, 339 const ConstTwoDMatrix &a, const KronProdAll &kp) 340 : FTensor(indor::along_col, PerTensorDimens(td, e).getNVX(), 341 a.nrows(), kp.ncols(), td.dimen()), 342 tdims(td, e) 343 { 344 kp.mult(a, *this); 345 } FPSTensor(const TensorDimens & td,const Equivalence & e,const ConstTwoDMatrix & a,const KronProdAllOptim & kp)346 FPSTensor(const TensorDimens &td, const Equivalence &e, 347 const ConstTwoDMatrix &a, const KronProdAllOptim &kp) 348 : FTensor(indor::along_col, PerTensorDimens(td, Permutation(e, kp.getPer())).getNVX(), 349 a.nrows(), kp.ncols(), td.dimen()), 350 tdims(td, e, kp.getPer()) 351 { 352 kp.mult(a, *this); 353 } FPSTensor(const TensorDimens & td,const Equivalence & e,const Permutation & p,const ConstTwoDMatrix & a,const KronProdAll & kp)354 FPSTensor(const TensorDimens &td, const Equivalence &e, const Permutation &p, 355 const ConstTwoDMatrix &a, const KronProdAll &kp) 356 : FTensor(indor::along_col, PerTensorDimens(td, Permutation(e, p)).getNVX(), 357 a.nrows(), kp.ncols(), td.dimen()), 358 tdims(td, e, p) 359 { 360 kp.mult(a, *this); 361 } FPSTensor(const TensorDimens & td,const Equivalence & e,const Permutation & p,const ConstTwoDMatrix & a,const KronProdAllOptim & kp)362 FPSTensor(const TensorDimens &td, const Equivalence &e, const Permutation &p, 363 const ConstTwoDMatrix &a, const KronProdAllOptim &kp) 364 : FTensor(indor::along_col, PerTensorDimens(td, Permutation(e, Permutation(p, kp.getPer()))).getNVX(), 365 a.nrows(), kp.ncols(), td.dimen()), 366 tdims(td, e, Permutation(p, kp.getPer())) 367 { 368 kp.mult(a, *this); 369 } 370 371 FPSTensor(const TensorDimens &td, const Equivalence &e, const Permutation &p, 372 const GSSparseTensor &t, const KronProdAll &kp); 373 374 FPSTensor(const FPSTensor &) = default; 375 FPSTensor(FPSTensor &&) = default; 376 377 void increment(IntSequence &v) const override; 378 void decrement(IntSequence &v) const override; 379 std::unique_ptr<UTensor> unfold() const override; 380 381 int getOffset(const IntSequence &v) const override; 382 void addTo(FGSTensor &out) const; 383 }; 384 385 #endif 386