1 /*
2 * Copyright © 2004-2011 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 "symmetry.hh"
22 #include "permutation.hh"
23 #include "tl_exception.hh"
24
25 #include <iostream>
26
27 /* This constructs an implied symmetry from a more general symmetry and
28 equivalence class. For example, let the general symmetry be y³u² and the
29 equivalence class is {0,4} picking up first and fifth variable, we
30 calculate symmetry corresponding to the picked variables. These are “yu”.
31 Thus the constructed sequence must be (1,1), meaning that we picked one
32 y and one u. */
33
Symmetry(const Symmetry & sy,const OrdSequence & cl)34 Symmetry::Symmetry(const Symmetry &sy, const OrdSequence &cl)
35 : IntSequence(sy.num())
36 {
37 const std::vector<int> &se = cl.getData();
38 TL_RAISE_IF(sy.dimen() <= se[se.size()-1],
39 "Sequence is not reachable by symmetry in IntSequence()");
40 for (int i = 0; i < size(); i++)
41 operator[](i) = 0;
42
43 for (int i : se)
44 operator[](sy.findClass(i))++;
45 }
46
47 /* Find a class of the symmetry containing a given index. */
48
49 int
findClass(int i) const50 Symmetry::findClass(int i) const
51 {
52 int j = 0;
53 int sum = 0;
54 do
55 {
56 sum += operator[](j);
57 j++;
58 }
59 while (j < size() && sum <= i);
60
61 return j-1;
62 }
63
64 /* The symmetry is full if it allows for any permutation of indices. It
65 means, that there is at most one non-zero index. */
66
67 bool
isFull() const68 Symmetry::isFull() const
69 {
70 int count = 0;
71 for (int i = 0; i < num(); i++)
72 if (operator[](i) != 0)
73 count++;
74 return count <= 1;
75 }
76
77 /* Construct a symiterator of given dimension, starting from the given
78 symmetry. */
79
symiterator(int dim_arg,Symmetry run_arg)80 symiterator::symiterator(int dim_arg, Symmetry run_arg)
81 : dim{dim_arg}, run(std::move(run_arg))
82 {
83 }
84
85 /* Here we move to the next symmetry. We do so only, if we are not at
86 the end. If length is 2, we increase lower index and decrease upper
87 index, otherwise we increase the subordinal symmetry. If we got to the
88 end, we recreate the subordinal symmetry set and set the subordinal
89 iterator to the beginning. */
90
91 symiterator &
operator ++()92 symiterator::operator++()
93 {
94 if (run[0] == dim)
95 run[0]++; // Jump to the past-the-end iterator
96 else if (run.size() == 2)
97 {
98 run[0]++;
99 run[1]--;
100 }
101 else
102 {
103 symiterator subit{dim-run[0], Symmetry(run, run.size()-1)};
104 ++subit;
105 if (run[1] == dim-run[0]+1)
106 {
107 run[0]++;
108 run[1] = 0;
109 /* subit is equal to the past-the-end iterator, so the range
110 2…(size()−1) is already set to 0 */
111 run[run.size()-1] = dim-run[0];
112 }
113 }
114 return *this;
115 }
116
InducedSymmetries(const Equivalence & e,const Symmetry & s)117 InducedSymmetries::InducedSymmetries(const Equivalence &e, const Symmetry &s)
118 {
119 for (const auto &i : e)
120 emplace_back(s, i);
121 }
122
123 // InducedSymmetries permuted constructor code
InducedSymmetries(const Equivalence & e,const Permutation & p,const Symmetry & s)124 InducedSymmetries::InducedSymmetries(const Equivalence &e, const Permutation &p,
125 const Symmetry &s)
126 {
127 for (int i = 0; i < e.numClasses(); i++)
128 {
129 auto it = e.find(p.getMap()[i]);
130 emplace_back(s, *it);
131 }
132 }
133
134 /* Debug print. */
135
136 void
print() const137 InducedSymmetries::print() const
138 {
139 std::cout << "Induced symmetries: " << size() << std::endl;
140 for (unsigned int i = 0; i < size(); i++)
141 operator[](i).print();
142 }
143