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 "int_sequence.hh"
22 #include "symmetry.hh"
23 #include "tl_exception.hh"
24 #include "pascal_triangle.hh"
25
26 #include <iostream>
27 #include <limits>
28 #include <numeric>
29
30 IntSequence
unfold(const Symmetry & sy) const31 IntSequence::unfold(const Symmetry &sy) const
32 {
33 IntSequence r(sy.dimen());
34 int k = 0;
35 for (int i = 0; i < sy.num(); i++)
36 for (int j = 0; j < sy[i]; j++, k++)
37 r[k] = operator[](i);
38 return r;
39 }
40
41 Symmetry
getSymmetry() const42 IntSequence::getSymmetry() const
43 {
44 Symmetry r(getNumDistinct());
45 int p = 0;
46 if (size() > 0)
47 r[p] = 1;
48 for (int i = 1; i < size(); i++)
49 {
50 if (operator[](i) != operator[](i-1))
51 p++;
52 r[p]++;
53 }
54 return r;
55 }
56
57 IntSequence
insert(int i) const58 IntSequence::insert(int i) const
59 {
60 IntSequence r(size()+1);
61 int j;
62 for (j = 0; j < size() && operator[](j) < i; j++)
63 r[j] = operator[](j);
64 r[j] = i;
65 for (; j < size(); j++)
66 r[j+1] = operator[](j);
67 return r;
68 }
69
70 IntSequence
insert(int i,int pos) const71 IntSequence::insert(int i, int pos) const
72 {
73 TL_RAISE_IF(pos < 0 || pos > size(),
74 "Wrong position for IntSequence::insert()");
75 IntSequence r(size()+1);
76 int j;
77 for (j = 0; j < pos; j++)
78 r[j] = operator[](j);
79 r[j] = i;
80 for (; j < size(); j++)
81 r[j+1] = operator[](j);
82 return r;
83 }
84
85 IntSequence &
operator =(const IntSequence & s)86 IntSequence::operator=(const IntSequence &s)
87 {
88 TL_RAISE_IF(length != s.length, "Wrong length for in-place IntSequence::operator=");
89 std::copy_n(s.data, length, data);
90 return *this;
91 }
92
93 IntSequence &
operator =(IntSequence && s)94 IntSequence::operator=(IntSequence &&s)
95 {
96 TL_RAISE_IF(length != s.length, "Wrong length for in-place IntSequence::operator=");
97 std::copy_n(s.data, length, data);
98 return *this;
99 }
100
101 bool
operator ==(const IntSequence & s) const102 IntSequence::operator==(const IntSequence &s) const
103 {
104 return std::equal(data, data+length,
105 s.data, s.data+s.length);
106 }
107
108 bool
operator <(const IntSequence & s) const109 IntSequence::operator<(const IntSequence &s) const
110 {
111 return std::lexicographical_compare(data, data+length,
112 s.data, s.data+s.length);
113 }
114
115 bool
lessEq(const IntSequence & s) const116 IntSequence::lessEq(const IntSequence &s) const
117 {
118 TL_RAISE_IF(size() != s.size(),
119 "Sequence with different lengths in IntSequence::lessEq");
120
121 int i = 0;
122 while (i < size() && operator[](i) <= s[i])
123 i++;
124 return (i == size());
125 }
126
127 bool
less(const IntSequence & s) const128 IntSequence::less(const IntSequence &s) const
129 {
130 TL_RAISE_IF(size() != s.size(),
131 "Sequence with different lengths in IntSequence::less");
132
133 int i = 0;
134 while (i < size() && operator[](i) < s[i])
135 i++;
136 return (i == size());
137 }
138
139 void
sort()140 IntSequence::sort()
141 {
142 std::sort(data, data+length);
143 }
144
145 void
monotone()146 IntSequence::monotone()
147 {
148 for (int i = 1; i < length; i++)
149 if (operator[](i-1) > operator[](i))
150 operator[](i) = operator[](i-1);
151 }
152
153 void
pmonotone(const Symmetry & s)154 IntSequence::pmonotone(const Symmetry &s)
155 {
156 int cum = 0;
157 for (int i = 0; i < s.num(); i++)
158 {
159 for (int j = cum + 1; j < cum + s[i]; j++)
160 if (operator[](j-1) > operator[](j))
161 operator[](j) = operator[](j-1);
162 cum += s[i];
163 }
164 }
165
166 int
sum() const167 IntSequence::sum() const
168 {
169 return std::accumulate(data, data+length, 0);
170 }
171
172 int
mult(int i1,int i2) const173 IntSequence::mult(int i1, int i2) const
174 {
175 return std::accumulate(data+i1, data+i2,
176 1, std::multiplies<int>());
177 }
178
179 int
getPrefixLength() const180 IntSequence::getPrefixLength() const
181 {
182 int i = 0;
183 while (i+1 < size() && operator[](i+1) == operator[](0))
184 i++;
185 return i+1;
186 }
187
188 int
getNumDistinct() const189 IntSequence::getNumDistinct() const
190 {
191 int res = 0;
192 if (length > 0)
193 res++;
194 for (int i = 1; i < length; i++)
195 if (operator[](i) != operator[](i-1))
196 res++;
197 return res;
198 }
199
200 int
getMax() const201 IntSequence::getMax() const
202 {
203 if (length == 0)
204 return std::numeric_limits<int>::min();
205 return *std::max_element(data, data+length);
206 }
207
208 void
add(int i)209 IntSequence::add(int i)
210 {
211 for (int j = 0; j < size(); j++)
212 operator[](j) += i;
213 }
214
215 void
add(int f,const IntSequence & s)216 IntSequence::add(int f, const IntSequence &s)
217 {
218 TL_RAISE_IF(size() != s.size(),
219 "Wrong sequence length in IntSequence::add");
220 for (int j = 0; j < size(); j++)
221 operator[](j) += f*s[j];
222 }
223
224 bool
isPositive() const225 IntSequence::isPositive() const
226 {
227 return std::all_of(data, data+length,
228 [](int x) { return x >= 0; });
229 }
230
231 bool
isConstant() const232 IntSequence::isConstant() const
233 {
234 if (length < 2)
235 return true;
236 return std::all_of(data+1, data+length,
237 [this](int x) { return x == operator[](0); });
238 }
239
240 bool
isSorted() const241 IntSequence::isSorted() const
242 {
243 return std::is_sorted(data, data+length);
244 }
245
246 /* Debug print. */
247
248 void
print() const249 IntSequence::print() const
250 {
251 std::cout << '[';
252 for (int i = 0; i < size(); i++)
253 std::cout << operator[](i) << ' ';
254 std::cout << ']' << std::endl;
255 }
256
257 /* Here we calculate the multinomial coefficients
258 ⎛ a ⎞
259 ⎝b₁,…,bₙ⎠ where a=b₁+…+bₙ
260
261 (the notation gives the name to the function: “n over seq(ence)”)
262
263 See:
264 https://en.wikipedia.org/wiki/Binomial_coefficient#Generalization_to_multinomials
265 https://en.wikipedia.org/wiki/Multinomial_theorem
266
267 For n=1, the coefficient is equal to 1.
268 For n=2, the multinomial coeffs correspond to the binomial coeffs, i.e. the binomial
269 ⎛a⎞ ⎛ a ⎞
270 ⎝b⎠ is equal to the multinomial ⎝b,a−b⎠
271
272 For n≥3, we have the identity
273 ⎛ a ⎞ ⎛b₁+b₂⎞ ⎛ a ⎞
274 ⎝b₁,…,bₙ⎠=⎝ b₁ ⎠·⎝b₁+b₂,b₃,…,bₙ⎠
275 (where the first factor on the right hand side is to be interpreted as a
276 binomial coefficient)
277
278 This number is exactly a number of unfolded indices corresponding to one
279 folded index, where the sequence b₁,…,bₙ is the symmetry of the index. This
280 can be easily seen if the multinomial coefficient is interpreted as the
281 number of unique permutations of a word, where ‘a’ is the length of the
282 word, ‘n’ is the number of distinct letters, and the bᵢ are the number of
283 repetitions of each letter. For example, for a symmetry of the form y⁴u²v³,
284 we want to compute the number of permutations of the word ‘yyyyuuvvv’.
285 ⎛ 9 ⎞
286 This is equal to the multinomial coefficient ⎝4,2,3⎠.
287 */
288 int
noverseq()289 IntSequence::noverseq()
290 {
291 if (size() == 0 || size() == 1)
292 return 1;
293 data[1] += data[0];
294 return PascalTriangle::noverk(data[1], data[0]) * IntSequence(*this, 1, size()).noverseq();
295 }
296