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 // Tensor concept. 22 23 /* Here we define a tensor class. A tensor is a mathematical object 24 corresponding to an (n+1)-dimensional array. An element of such array is 25 denoted [B]_α₁…αₙ^β, where β is a special index and α₁…αₙ are other indices. 26 The class Tensor and its subclasses view such array as a 2D matrix, where β 27 corresponds to one dimension, and α₁…αₙ unfold to the other dimension. 28 Whether β corresponds to rows or columns is decided by tensor subclasses, 29 however, most of our tensors will have rows indexed by β, and α₁…αₙ will 30 unfold column-wise. 31 32 There might be some symmetries in the tensor data. For instance, if α₁ is 33 interchanged with α₃ and the other elements remain equal for all possible αᵢ 34 and β, then there is a symmetry of α₁ and α₃. 35 36 For any symmetry, there are basically two possible storages of the 37 data. The first is unfolded storage, which stores all elements 38 regardless the symmetry. The other storage type is folded, which 39 stores only elements which do not repeat. We declare abstract classes 40 for unfolded and folded tensor alike. 41 42 Also, here we also define a concept of tensor index which is the n-tuple 43 α₁…αₙ. It is an iterator, which iterates in dependence of symmetry and 44 storage of the underlying tensor. 45 46 Although we do not decide about possible symmetries at this point, it is 47 worth noting that we implement two kinds of symmetries in subclasses. The 48 first one is a full symmetry where all indices are interchangeable. The 49 second one is a generalization of the first, where there are a few groups of 50 indices interchangeable within a group and not across. Moreover, the groups 51 are required to be consequent partitions of the index n-tuple. For example, 52 we do not allow α₁ to be interchangeable with α₃ and not with α₂ at the same 53 time. 54 55 However, some intermediate results are, in fact, tensors with a symmetry not 56 fitting to our concept. We develop the tensor abstraction for it, but these 57 objects are not used very often. They have limited usage due to their 58 specialized constructor. */ 59 60 #ifndef TENSOR_H 61 #define TENSOR_H 62 63 #include "int_sequence.hh" 64 #include "twod_matrix.hh" 65 66 #include <memory> 67 #include <iostream> 68 69 /* Here is the Tensor class, which is nothing else than a simple subclass of 70 TwoDMatrix. The unique semantically new member is ‘dim’ which is tensor 71 dimension (length of α₁…αₙ). We also declare increment(), decrement() and 72 getOffset() methods as pure virtual. 73 74 We also add members for index begin and index end. This is useful, since 75 begin() and end() methods do not return instance but only references, which 76 prevent making additional copy of index (for example in for cycles as ‘in != 77 end()’ which would do a copy of index for each cycle). The index begin 78 ‘in_beg’ is constructed as a sequence of all zeros, and ‘in_end’ is 79 constructed from the sequence ‘last’ passed to the constructor, since it 80 depends on subclasses. Also we have to say, along what coordinate is the 81 multidimensional index. This is used only for initialization of ‘in_end’. */ 82 83 class Tensor : public TwoDMatrix 84 { 85 public: 86 enum class indor { along_row, along_col }; 87 88 /* The index represents n-tuple α₁…αₙ. Since its movement is dependent on the 89 underlying tensor (with storage and symmetry), we maintain a reference to 90 that tensor, we maintain the n-tuple (or coordinates) as IntSequence and 91 also we maintain the offset number (column, or row) of the index in the 92 tensor. The reference is const, since we do not need to change data 93 through the index. 94 95 Here we require the Tensor to implement increment() and decrement() 96 methods, which calculate following and preceding n-tuple. Also, we need to 97 calculate offset number from the given coordinates, so the tensor must 98 implement method getOffset(). This method is used only in construction of 99 the index from the given coordinates. As the index is created, the offset 100 is automatically incremented, and decremented together with index. The 101 getOffset() method can be relatively computationally complex. This must be 102 kept in mind. Also we generally suppose that n-tuple of all zeros is the 103 first offset (first columns or row). 104 105 What follows is a definition of index class, the only interesting point is 106 operator==() which decides only according to offset, not according to the 107 coordinates. This is useful since there can be more than one of coordinate 108 representations of past-the-end index. */ 109 class index 110 { 111 const Tensor &tensor; 112 int offset; 113 IntSequence coor; 114 public: index(const Tensor & t,int n)115 index(const Tensor &t, int n) 116 : tensor(t), offset(0), coor(n, 0) 117 { 118 } index(const Tensor & t,IntSequence cr,int c)119 index(const Tensor &t, IntSequence cr, int c) 120 : tensor(t), offset(c), coor(std::move(cr)) 121 { 122 } index(const Tensor & t,IntSequence cr)123 index(const Tensor &t, IntSequence cr) 124 : tensor(t), offset(tensor.getOffset(cr)), coor(std::move(cr)) 125 { 126 } 127 index(const index &) = default; 128 index(index &&) = default; 129 index &operator=(const index &) = delete; 130 index &operator=(index &&) = delete; 131 index & operator ++()132 operator++() 133 { 134 tensor.increment(coor); 135 offset++; 136 return *this; 137 } 138 index & operator --()139 operator--() 140 { 141 tensor.decrement(coor); 142 offset--; 143 return *this; 144 } 145 int operator *() const146 operator*() const 147 { 148 return offset; 149 } 150 bool operator ==(const index & n) const151 operator==(const index &n) const 152 { 153 return offset == n.offset; 154 } 155 bool operator !=(const index & n) const156 operator!=(const index &n) const 157 { 158 return offset != n.offset; 159 } 160 const IntSequence & getCoor() const161 getCoor() const 162 { 163 return coor; 164 } 165 void print() const166 print() const 167 { 168 std::cout << offset << ": "; 169 coor.print(); 170 } 171 }; 172 173 protected: 174 const index in_beg; 175 const index in_end; 176 int dim; 177 public: Tensor(indor io,IntSequence last,int r,int c,int d)178 Tensor(indor io, IntSequence last, int r, int c, int d) 179 : TwoDMatrix(r, c), 180 in_beg(*this, d), 181 in_end(*this, std::move(last), (io == indor::along_row) ? r : c), 182 dim(d) 183 { 184 } Tensor(indor io,IntSequence first,IntSequence last,int r,int c,int d)185 Tensor(indor io, IntSequence first, IntSequence last, 186 int r, int c, int d) 187 : TwoDMatrix(r, c), 188 in_beg(*this, std::move(first), 0), 189 in_end(*this, std::move(last), (io == indor::along_row) ? r : c), 190 dim(d) 191 { 192 } Tensor(int first_row,int num,Tensor & t)193 Tensor(int first_row, int num, Tensor &t) 194 : TwoDMatrix(first_row, num, t), 195 in_beg(t.in_beg), 196 in_end(t.in_end), 197 dim(t.dim) 198 { 199 } Tensor(const Tensor & t)200 Tensor(const Tensor &t) 201 : TwoDMatrix(t), 202 in_beg(*this, t.in_beg.getCoor(), *(t.in_beg)), 203 in_end(*this, t.in_end.getCoor(), *(t.in_end)), 204 dim(t.dim) 205 { 206 } 207 Tensor(Tensor &&) = default; 208 ~Tensor() override = default; 209 210 Tensor &operator=(const Tensor &) = delete; 211 Tensor &operator=(Tensor &&) = delete; 212 213 virtual void increment(IntSequence &v) const = 0; 214 virtual void decrement(IntSequence &v) const = 0; 215 virtual int getOffset(const IntSequence &v) const = 0; 216 int dimen() const217 dimen() const 218 { 219 return dim; 220 } 221 222 const index & begin() const223 begin() const 224 { 225 return in_beg; 226 } 227 const index & end() const228 end() const 229 { 230 return in_end; 231 } 232 }; 233 234 /* Here is an abstraction for unfolded tensor. We provide a pure virtual method 235 fold() which returns a new instance of folded tensor of the same symmetry. 236 Also we provide static methods for incrementing and decrementing an index 237 with full symmetry and general symmetry as defined above. */ 238 239 class FTensor; 240 class UTensor : public Tensor 241 { 242 public: UTensor(indor io,IntSequence last,int r,int c,int d)243 UTensor(indor io, IntSequence last, int r, int c, int d) 244 : Tensor(io, std::move(last), r, c, d) 245 { 246 } 247 UTensor(const UTensor &) = default; 248 UTensor(UTensor &&) = default; UTensor(int first_row,int num,UTensor & t)249 UTensor(int first_row, int num, UTensor &t) 250 : Tensor(first_row, num, t) 251 { 252 } 253 ~UTensor() override = default; 254 virtual std::unique_ptr<FTensor> fold() const = 0; 255 256 UTensor &operator=(const UTensor &) = delete; 257 UTensor &operator=(UTensor &&) = delete; 258 259 static void increment(IntSequence &v, int nv); 260 static void decrement(IntSequence &v, int nv); 261 static void increment(IntSequence &v, const IntSequence &nvmx); 262 static void decrement(IntSequence &v, const IntSequence &nvmx); 263 static int getOffset(const IntSequence &v, int nv); 264 static int getOffset(const IntSequence &v, const IntSequence &nvmx); 265 }; 266 267 /* This is an abstraction for folded tensor. It only provides a method 268 unfold(), which returns the unfolded version of the same symmetry, and 269 static methods for decrementing indices. 270 271 We also provide static methods for decrementing the IntSequence in 272 folded fashion and also calculating an offset for a given 273 IntSequence. However, this is relatively complex calculation, so 274 this should be avoided if possible. */ 275 276 class FTensor : public Tensor 277 { 278 public: FTensor(indor io,IntSequence last,int r,int c,int d)279 FTensor(indor io, IntSequence last, int r, int c, int d) 280 : Tensor(io, std::move(last), r, c, d) 281 { 282 } 283 FTensor(const FTensor &) = default; 284 FTensor(FTensor &&) = default; FTensor(int first_row,int num,FTensor & t)285 FTensor(int first_row, int num, FTensor &t) 286 : Tensor(first_row, num, t) 287 { 288 } 289 ~FTensor() override = default; 290 virtual std::unique_ptr<UTensor> unfold() const = 0; 291 292 FTensor &operator=(const FTensor &) = delete; 293 FTensor &operator=(FTensor &&) = delete; 294 295 static void decrement(IntSequence &v, int nv); 296 static int getOffset(const IntSequence & v,int nv)297 getOffset(const IntSequence &v, int nv) 298 { 299 IntSequence vtmp(v); 300 return getOffsetRecurse(vtmp, nv); 301 } 302 private: 303 static int getOffsetRecurse(IntSequence &v, int nv); 304 }; 305 306 #endif 307