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 "tensor.hh"
22 #include "tl_exception.hh"
23 #include "tl_static.hh"
24 #include "pascal_triangle.hh"
25
26 /* Here we increment a given sequence within full symmetry given by ‘nv’, which
27 is number of variables in each dimension. The underlying tensor is unfolded,
28 so we increase the rightmost by one, and if it is ‘nv’ we zero it and
29 increase the next one to the left. */
30
31 void
increment(IntSequence & v,int nv)32 UTensor::increment(IntSequence &v, int nv)
33 {
34 if (v.size() == 0)
35 return;
36 int i = v.size()-1;
37 v[i]++;
38 while (i > 0 && v[i] == nv)
39 {
40 v[i] = 0;
41 v[--i]++;
42 }
43 }
44
45 /* This is dual to UTensor::increment(IntSequence& v, int nv). */
46
47 void
decrement(IntSequence & v,int nv)48 UTensor::decrement(IntSequence &v, int nv)
49 {
50 if (v.size() == 0)
51 return;
52 int i = v.size()-1;
53 v[i]--;
54 while (i > 0 && v[i] == -1)
55 {
56 v[i] = nv -1;
57 v[--i]--;
58 }
59 }
60
61 /* Here we increment index for general symmetry for unfolded storage. The
62 sequence ‘nvmx’ assigns for each coordinate a number of variables. Since the
63 storage is unfolded, we do not need information about what variables are
64 symmetric, everything necessary is given by ‘nvmx’. */
65
66 void
increment(IntSequence & v,const IntSequence & nvmx)67 UTensor::increment(IntSequence &v, const IntSequence &nvmx)
68 {
69 if (v.size() == 0)
70 return;
71 int i = v.size()-1;
72 v[i]++;
73 while (i > 0 && v[i] == nvmx[i])
74 {
75 v[i] = 0;
76 v[--i]++;
77 }
78 }
79
80 /* This is a dual code to UTensor::increment(IntSequence& v, const IntSequence&
81 nvmx). */
82
83 void
decrement(IntSequence & v,const IntSequence & nvmx)84 UTensor::decrement(IntSequence &v, const IntSequence &nvmx)
85 {
86 if (v.size() == 0)
87 return;
88 int i = v.size()-1;
89 v[i]--;
90 while (i > 0 && v[i] == -1)
91 {
92 v[i] = nvmx[i] -1;
93 v[--i]--;
94 }
95 }
96
97 /* Here we return an offset for a given coordinates of unfolded full symmetry
98 tensor. This is easy. */
99
100 int
getOffset(const IntSequence & v,int nv)101 UTensor::getOffset(const IntSequence &v, int nv)
102 {
103 int res = 0;
104 for (int i = 0; i < v.size(); i++)
105 {
106 res *= nv;
107 res += v[i];
108 }
109 return res;
110 }
111
112 /* Also easy. */
113
114 int
getOffset(const IntSequence & v,const IntSequence & nvmx)115 UTensor::getOffset(const IntSequence &v, const IntSequence &nvmx)
116 {
117 int res = 0;
118 for (int i = 0; i < v.size(); i++)
119 {
120 res *= nvmx[i];
121 res += v[i];
122 }
123 return res;
124 }
125
126 /* Decrementing of coordinates of folded index is not that easy. Note that if a
127 trailing part of coordinates is (b,a,a,a) (for instance) with b<a, then a
128 preceding coordinates are (b,a−1,n−1,n−1), where n is a number of variables
129 ‘nv’. So we find the left most element which is equal to the last element,
130 decrease it by one, and then set all elements to the right to n−1. */
131
132 void
decrement(IntSequence & v,int nv)133 FTensor::decrement(IntSequence &v, int nv)
134 {
135 int i = v.size()-1;
136 while (i > 0 && v[i-1] == v[i])
137 i--;
138 v[i]--;
139 for (int j = i+1; j < v.size(); j++)
140 v[j] = nv-1;
141 }
142
143 /* This calculates order of the given index of our ordering of indices. In
144 order to understand how it works, let us take number of variables n and
145 dimension k, and write down all the possible combinations of indices in our
146 ordering. For example for n=4 and k=3, the sequence looks as the following
147 (in column-major order):
148
149 000 111 222 333
150 001 112 223
151 002 113 233
152 003 122
153 011 123
154 012 133
155 013
156 022
157 023
158 033
159
160 Now observe, that a number of sequences starting with zero is the same as
161 the total number of sequences with the same number of variables but with
162 dimension minus one. More generally, if Sₙ,ₖ denotes the number of indices
163 of n variables and dimension k, then the number of indices beginning with m
164 is exactly Sₙ₋ₘ,ₖ₋₁. This is because m can be subtracted from all items, and
165 we obtain the sequence of indices of n−m variables. So we have the formula:
166
167 Sₙ,ₖ = Sₙ,ₖ₋₁ + Sₙ₋₁,ₖ₋₁ + … + S₁,ₖ₋₁
168
169 Now it is easy to calculate the offset of an index of the form (m,…,m). It
170 is the sum of all above it, that is Sₙ,ₖ₋₁ + … + Sₙ₋ₘ,ₖ₋₁. Also, since Sₙ,ₖ
171 is the number of ways to choose k elements from a set of n elements with
172 repetitions allowed, it is well known that:
173 ⎛n+k−1⎞
174 Sₙ,ₖ = ⎝ k ⎠.
175 Using the above formula, one can therefore calculate the offset of (m,…,m):
176 ⎛n+k−1⎞ ⎛n−m+k−1⎞
177 ⎝ k ⎠−⎝ k ⎠
178
179 The offset of general index (m₁,m₂,…,mₖ) is calculated recursively, since it
180 is the offset of (m₁,…,m₁) for n variables plus the offset of
181 (m₂−m₁,m₃−m₁,…,mₖ−m₁) for n−m₁ variables. */
182
183 int
getOffsetRecurse(IntSequence & v,int nv)184 FTensor::getOffsetRecurse(IntSequence &v, int nv)
185 {
186 if (v.size() == 0)
187 return 0;
188 int prefix = v.getPrefixLength();
189 int m = v[0];
190 int k = v.size();
191 int s1 = PascalTriangle::noverk(nv+k-1, k) - PascalTriangle::noverk(nv-m+k-1, k);
192 IntSequence subv(v, prefix, k);
193 subv.add(-m);
194 int s2 = getOffsetRecurse(subv, nv-m);
195 return s1+s2;
196 }
197