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 #include "fs_tensor.hh"
22 #include "gs_tensor.hh"
23 #include "sparse_tensor.hh"
24 #include "rfs_tensor.hh"
25 #include "tl_exception.hh"
26 #include "pascal_triangle.hh"
27
28 /* This constructs a fully symmetric tensor as given by the contraction:
29
30 [g_yⁿ]_α₁…αₙ = [t_yⁿ⁺¹]_α₁…αₙβ [x]^β
31
32 We go through all columns of output tensor [g] and for each column
33 we cycle through all variables, insert a variable to the column
34 coordinates obtaining a column of tensor [t]. The column is multiplied
35 by an appropriate item of x and added to the column of [g] tensor. */
36
FFSTensor(const FFSTensor & t,const ConstVector & x)37 FFSTensor::FFSTensor(const FFSTensor &t, const ConstVector &x)
38 : FTensor(indor::along_col, IntSequence(t.dimen()-1, t.nvar()),
39 t.nrows(), calcMaxOffset(t.nvar(), t.dimen()-1), t.dimen()-1),
40 nv(t.nvar())
41 {
42 TL_RAISE_IF(t.dimen() < 1,
43 "Wrong dimension for tensor contraction of FFSTensor");
44 TL_RAISE_IF(t.nvar() != x.length(),
45 "Wrong number of variables for tensor contraction of FFSTensor");
46
47 zeros();
48
49 for (Tensor::index to = begin(); to != end(); ++to)
50 for (int i = 0; i < nvar(); i++)
51 {
52 IntSequence from_ind(to.getCoor().insert(i));
53 Tensor::index from(t, from_ind);
54 addColumn(x[i], t, *from, *to);
55 }
56 }
57
58 /* This returns number of indices for folded tensor with full
59 symmetry. Let n be a number of variables and d the
60 ⎛n+d-1⎞
61 dimension dim. Then the number of indices is ⎝ d ⎠.
62 */
63
64 int
calcMaxOffset(int nvar,int d)65 FFSTensor::calcMaxOffset(int nvar, int d)
66 {
67 if (nvar == 0 && d == 0)
68 return 1;
69 if (nvar == 0 && d > 0)
70 return 0;
71 return PascalTriangle::noverk(nvar + d - 1, d);
72 }
73
74 /* The conversion from sparse tensor is clear. We go through all the
75 tensor and write to the dense what is found. */
FFSTensor(const FSSparseTensor & t)76 FFSTensor::FFSTensor(const FSSparseTensor &t)
77 : FTensor(indor::along_col, IntSequence(t.dimen(), t.nvar()),
78 t.nrows(), calcMaxOffset(t.nvar(), t.dimen()), t.dimen()),
79 nv(t.nvar())
80 {
81 zeros();
82 for (const auto &it : t.getMap())
83 {
84 index ind(*this, it.first);
85 get(it.second.first, *ind) = it.second.second;
86 }
87 }
88
89 /* The conversion from unfolded copies only columns of respective
90 coordinates. So we go through all the columns in the folded tensor
91 (this), make an index of the unfolded vector from coordinates, and
92 copy the column. */
93
FFSTensor(const UFSTensor & ut)94 FFSTensor::FFSTensor(const UFSTensor &ut)
95 : FTensor(indor::along_col, IntSequence(ut.dimen(), ut.nvar()),
96 ut.nrows(), calcMaxOffset(ut.nvar(), ut.dimen()), ut.dimen()),
97 nv(ut.nvar())
98 {
99 for (index in = begin(); in != end(); ++in)
100 {
101 index src(ut, in.getCoor());
102 copyColumn(ut, *src, *in);
103 }
104 }
105
106 /* Here just make a new instance and return the reference. */
107 std::unique_ptr<UTensor>
unfold() const108 FFSTensor::unfold() const
109 {
110 return std::make_unique<UFSTensor>(*this);
111 }
112
113 /* Incrementing is easy. We have to increment by calling static method
114 UTensor::increment() first. In this way, we have coordinates of
115 unfolded tensor. Then we have to skip to the closest folded index
116 which corresponds to monotonizeing the integer sequence. */
117
118 void
increment(IntSequence & v) const119 FFSTensor::increment(IntSequence &v) const
120 {
121 TL_RAISE_IF(v.size() != dimen(),
122 "Wrong input/output vector size in FFSTensor::increment");
123
124 UTensor::increment(v, nv);
125 v.monotone();
126 }
127
128 /* Decrement calls static FTensor::decrement(). */
129
130 void
decrement(IntSequence & v) const131 FFSTensor::decrement(IntSequence &v) const
132 {
133 TL_RAISE_IF(v.size() != dimen(),
134 "Wrong input/output vector size in FFSTensor::decrement");
135
136 FTensor::decrement(v, nv);
137 }
138
139 int
getOffset(const IntSequence & v) const140 FFSTensor::getOffset(const IntSequence &v) const
141 {
142 TL_RAISE_IF(v.size() != dimen(),
143 "Wrong input vector size in FFSTensor::getOffset");
144
145 return FTensor::getOffset(v, nv);
146 }
147
148 /* Here we add a general symmetry tensor to the (part of) full symmetry
149 tensor provided that the unique variable of the full symmetry tensor
150 is a stack of variables from the general symmetry tensor.
151
152 We check for the dimensions and number of variables. Then we calculate
153 a shift of coordinates when going from the general symmetry tensor to
154 full symmetry (it corresponds to shift of coordinates induces by
155 stacking the variables). Then we add the appropriate columns by going
156 through the columns in general symmetry, adding the shift and sorting. */
157
158 void
addSubTensor(const FGSTensor & t)159 FFSTensor::addSubTensor(const FGSTensor &t)
160 {
161 TL_RAISE_IF(dimen() != t.getDims().dimen(),
162 "Wrong dimensions for FFSTensor::addSubTensor");
163 TL_RAISE_IF(nvar() != t.getDims().getNVS().sum(),
164 "Wrong nvs for FFSTensor::addSubTensor");
165
166 // set shift for addSubTensor()
167 /* Code shared with UFSTensor::addSubTensor() */
168 IntSequence shift_pre(t.getSym().num(), 0);
169 for (int i = 1; i < t.getSym().num(); i++)
170 shift_pre[i] = shift_pre[i-1]+t.getDims().getNVS()[i-1];
171 IntSequence shift(shift_pre.unfold(t.getSym()));
172
173 for (Tensor::index ind = t.begin(); ind != t.end(); ++ind)
174 {
175 IntSequence c(ind.getCoor());
176 c.add(1, shift);
177 c.sort();
178 Tensor::index tar(*this, c);
179 addColumn(t, *ind, *tar);
180 }
181 }
182
183 // UFSTensor contraction constructor
184 /* This is a bit more straightforward than FFSTensor contraction constructor.
185 We do not add column by column but we do it by submatrices due to
186 regularity of the unfolded tensor. */
187
UFSTensor(const UFSTensor & t,const ConstVector & x)188 UFSTensor::UFSTensor(const UFSTensor &t, const ConstVector &x)
189 : UTensor(indor::along_col, IntSequence(t.dimen()-1, t.nvar()),
190 t.nrows(), calcMaxOffset(t.nvar(), t.dimen()-1), t.dimen()-1),
191 nv(t.nvar())
192 {
193 TL_RAISE_IF(t.dimen() < 1,
194 "Wrong dimension for tensor contraction of UFSTensor");
195 TL_RAISE_IF(t.nvar() != x.length(),
196 "Wrong number of variables for tensor contraction of UFSTensor");
197
198 zeros();
199
200 for (int i = 0; i < ncols(); i++)
201 {
202 ConstTwoDMatrix tpart(t, i *nvar(), nvar());
203 Vector outcol{getCol(i)};
204 tpart.multaVec(outcol, x);
205 }
206 }
207
208 /* Here we convert folded full symmetry tensor to unfolded. We copy all
209 columns of folded tensor, and then call unfoldData(). */
210
UFSTensor(const FFSTensor & ft)211 UFSTensor::UFSTensor(const FFSTensor &ft)
212 : UTensor(indor::along_col, IntSequence(ft.dimen(), ft.nvar()),
213 ft.nrows(), calcMaxOffset(ft.nvar(), ft.dimen()), ft.dimen()),
214 nv(ft.nvar())
215 {
216 for (index src = ft.begin(); src != ft.end(); ++src)
217 {
218 index in(*this, src.getCoor());
219 copyColumn(ft, *src, *in);
220 }
221 unfoldData();
222 }
223
224 std::unique_ptr<FTensor>
fold() const225 UFSTensor::fold() const
226 {
227 return std::make_unique<FFSTensor>(*this);
228 }
229
230 // UFSTensor increment and decrement
231 /* Here we just call UTensor respective static methods. */
232 void
increment(IntSequence & v) const233 UFSTensor::increment(IntSequence &v) const
234 {
235 TL_RAISE_IF(v.size() != dimen(),
236 "Wrong input/output vector size in UFSTensor::increment");
237
238 UTensor::increment(v, nv);
239 }
240
241 void
decrement(IntSequence & v) const242 UFSTensor::decrement(IntSequence &v) const
243 {
244 TL_RAISE_IF(v.size() != dimen(),
245 "Wrong input/output vector size in UFSTensor::decrement");
246
247 UTensor::decrement(v, nv);
248 }
249
250 int
getOffset(const IntSequence & v) const251 UFSTensor::getOffset(const IntSequence &v) const
252 {
253 TL_RAISE_IF(v.size() != dimen(),
254 "Wrong input vector size in UFSTensor::getOffset");
255
256 return UTensor::getOffset(v, nv);
257 }
258
259 /* This is very similar to FFSTensor::addSubTensor(). The
260 only difference is the addition. We go through all columns in the full
261 symmetry tensor and cancel the shift. If the coordinates after the
262 cancellation are positive, we find the column in the general symmetry
263 tensor, and add it. */
264
265 void
addSubTensor(const UGSTensor & t)266 UFSTensor::addSubTensor(const UGSTensor &t)
267 {
268 TL_RAISE_IF(dimen() != t.getDims().dimen(),
269 "Wrong dimensions for UFSTensor::addSubTensor");
270 TL_RAISE_IF(nvar() != t.getDims().getNVS().sum(),
271 "Wrong nvs for UFSTensor::addSubTensor");
272
273 // set shift for addSubTensor()
274 /* Code shared with FFSTensor::addSubTensor() */
275 IntSequence shift_pre(t.getSym().num(), 0);
276 for (int i = 1; i < t.getSym().num(); i++)
277 shift_pre[i] = shift_pre[i-1]+t.getDims().getNVS()[i-1];
278 IntSequence shift(shift_pre.unfold(t.getSym()));
279
280 for (Tensor::index tar = begin(); tar != end(); ++tar)
281 {
282 IntSequence c(tar.getCoor());
283 c.sort();
284 c.add(-1, shift);
285 if (c.isPositive() && c.less(t.getDims().getNVX()))
286 {
287 Tensor::index from(t, c);
288 addColumn(t, *from, *tar);
289 }
290 }
291 }
292
293 /* Here we go through all columns, find a column of folded index, and
294 then copy the column data. Finding the index is done by sorting the
295 integer sequence. */
296
297 void
unfoldData()298 UFSTensor::unfoldData()
299 {
300 for (index in = begin(); in != end(); ++in)
301 {
302 IntSequence v(in.getCoor());
303 v.sort();
304 copyColumn(*index(*this, v), *in);
305 }
306 }
307