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 "permutation.hh"
22 #include "tl_exception.hh"
23
24 /* This is easy, we simply apply the map in the fashion s∘m */
25
26 void
apply(const IntSequence & src,IntSequence & tar) const27 Permutation::apply(const IntSequence &src, IntSequence &tar) const
28 {
29 TL_RAISE_IF(src.size() != permap.size() || tar.size() != permap.size(),
30 "Wrong sizes of input or output in Permutation::apply");
31 for (int i = 0; i < permap.size(); i++)
32 tar[i] = src[permap[i]];
33 }
34
35 void
apply(IntSequence & tar) const36 Permutation::apply(IntSequence &tar) const
37 {
38 IntSequence tmp(tar);
39 apply(tmp, tar);
40 }
41
42 void
inverse()43 Permutation::inverse()
44 {
45 IntSequence former(permap);
46 for (int i = 0; i < size(); i++)
47 permap[former[i]] = i;
48 }
49
50 /* Here we find a number of trailing indices which are identical with
51 the permutation. */
52
53 int
tailIdentity() const54 Permutation::tailIdentity() const
55 {
56 int i = permap.size();
57 while (i > 0 && permap[i-1] == i-1)
58 i--;
59 return permap.size() - i;
60 }
61
62 /* This calculates a map which corresponds to sorting in the following
63 sense: $(sorted s)∘m = s, where s is a given sequence.
64
65 We go through ‘s’ and find an the same item in sorted ‘s’. We
66 construct the ‘permap’ from the found pair of indices. We have to be
67 careful, to not assign to two positions in ‘s’ the same position in
68 sorted ‘s’, so we maintain a bitmap ‘flag’, in which we remember
69 indices from the sorted ‘s’ already assigned. */
70
71 void
computeSortingMap(const IntSequence & s)72 Permutation::computeSortingMap(const IntSequence &s)
73 {
74 IntSequence srt(s);
75 srt.sort();
76 IntSequence flags(s.size(), 0);
77
78 for (int i = 0; i < s.size(); i++)
79 {
80 int j = 0;
81 while (j < s.size() && (flags[j] || srt[j] != s[i]))
82 j++;
83 TL_RAISE_IF(j == s.size(),
84 "Internal algorithm error in Permutation::computeSortingMap");
85 flags[j] = 1;
86 permap[i] = j;
87 }
88 }
89
PermutationSet()90 PermutationSet::PermutationSet()
91 {
92 pers.emplace_back(1);
93 }
94
PermutationSet(const PermutationSet & sp,int n)95 PermutationSet::PermutationSet(const PermutationSet &sp, int n)
96 : order(n), size(n*sp.size)
97 {
98 TL_RAISE_IF(n != sp.order+1,
99 "Wrong new order in PermutationSet constructor");
100
101 for (int i = 0; i < sp.size; i++)
102 for (int j = 0; j < order; j++)
103 pers.emplace_back(sp.pers[i], j);
104 }
105
106 std::vector<Permutation>
getPreserving(const IntSequence & s) const107 PermutationSet::getPreserving(const IntSequence &s) const
108 {
109 TL_RAISE_IF(s.size() != order,
110 "Wrong sequence length in PermutationSet::getPreserving");
111
112 std::vector<Permutation> res;
113 IntSequence tmp(s.size());
114 for (int i = 0; i < size; i++)
115 {
116 pers[i].apply(s, tmp);
117 if (s == tmp)
118 res.push_back(pers[i]);
119 }
120
121 return res;
122 }
123
PermutationBundle(int nmax)124 PermutationBundle::PermutationBundle(int nmax)
125 {
126 nmax = std::max(nmax, 1);
127 generateUpTo(nmax);
128 }
129
130 const PermutationSet &
get(int n) const131 PermutationBundle::get(int n) const
132 {
133 if (n > static_cast<int>(bundle.size()) || n < 1)
134 {
135 TL_RAISE("Permutation set not found in PermutationSet::get");
136 return bundle[0];
137 }
138 else
139 return bundle[n-1];
140 }
141
142 void
generateUpTo(int nmax)143 PermutationBundle::generateUpTo(int nmax)
144 {
145 if (bundle.size() == 0)
146 bundle.emplace_back();
147
148 int curmax = bundle.size();
149 for (int n = curmax+1; n <= nmax; n++)
150 bundle.emplace_back(bundle.back(), n);
151 }
152