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 // Equivalences.
22 
23 /* Here we define an equivalence of a set of integers {0, 1, …, k-1}.
24    The purpose is clear, in the tensor library we often iterate
25    through all equivalences and sum matrices. We need an abstraction for
26    an equivalence class, equivalence and a set of all equivalences.
27 
28    The equivalence class (which is basically a set of integers) is here
29    implemented as ordered integer sequence. The ordered sequence is not
30    implemented via IntSequence, but via vector<int> since we need
31    insertions. The equivalence is implemented as an ordered list of
32    equivalence classes, and equivalence set is a list of equivalences.
33 
34    The ordering of the equivalence classes within an equivalence is very
35    important. For instance, if we iterate through equivalences for k=5
36    and pickup some equivalence class, say { {0,4}, {1,2}, {3} }, we
37    then evaluate something like:
38 
39    [B_y²u³]_α₁α₂β₁β₂β₃ = … + [f_z³]_γ₁γ₂γ₃ [z_yu]^γ₁_α₁β₃ [z_yu]^γ₂_α₂β₁ [zᵤ]^γ₃_β₂ + …
40 
41    If the tensors are unfolded, we can evaluate this expression as
42    f_z³·(z_yu ⊗ z_yu ⊗ zᵤ)·P, where P is a suitable permutation of columns of
43    the expressions, which permutes them so that the index
44    (α₁,β₃,α₂,β₁,β₂) would go to (α₁,α₂,β₁,β₂,β₃).
45 
46    The permutation P can be very ineffective (copying great amount of
47    small chunks of data) if the equivalence class ordering is chosen
48    badly. However, we do not provide any heuristic minimizing a total
49    time spent in all permutations. We choose an ordering which orders the
50    classes according to their averages, and according to the smallest
51    equivalence class element if the averages are the same. */
52 
53 #ifndef EQUIVALENCE_H
54 #define EQUIVALENCE_H
55 
56 #include "int_sequence.hh"
57 
58 #include <vector>
59 #include <list>
60 #include <string>
61 
62 /* Here is the abstraction for an equivalence class. We implement it as
63    vector<int>. We have a constructor for empty class, copy
64    constructor. What is important here is the ordering operator
65    operator<() and methods for addition of an integer, and addition of
66    another sequence. Also we provide method has() which returns true if a
67    given integer is contained. */
68 
69 class OrdSequence
70 {
71   std::vector<int> data;
72 public:
OrdSequence()73   OrdSequence() : data()
74   {
75   }
76   bool operator==(const OrdSequence &s) const;
77   int operator[](int i) const;
78   bool operator<(const OrdSequence &s) const;
79   const std::vector<int> &
getData() const80   getData() const
81   {
82     return data;
83   }
84   int
length() const85   length() const
86   {
87     return data.size();
88   }
89   void add(int i);
90   void add(const OrdSequence &s);
91   bool has(int i) const;
92   void print(const std::string &prefix) const;
93 private:
94   double average() const;
95 };
96 
97 /* Here is the abstraction for the equivalence. It is a list of
98    equivalence classes. Also we remember n, which is a size of
99    underlying set {0, 1, …, n-1}.
100 
101    Method trace() “prints” the equivalence into the integer sequence. */
102 
103 class Permutation;
104 class Equivalence
105 {
106 private:
107   int n;
108   std::list<OrdSequence> classes;
109 public:
110   using const_seqit = std::list<OrdSequence>::const_iterator;
111   using seqit = std::list<OrdSequence>::iterator;
112 
113   // Constructs { {0}, {1}, …, {n-1} }
114   explicit Equivalence(int num);
115   // Constructs { {0,1,…,n-1 } }
116   Equivalence(int num, const std::string &dummy);
117   // Copy constructor plus gluing i1 and i2 in one class
118   Equivalence(const Equivalence &e, int i1, int i2);
119 
120   bool operator==(const Equivalence &e) const;
121   bool
operator !=(const Equivalence & e) const122   operator!=(const Equivalence &e) const
123   {
124     return !operator==(e);
125   }
126   int
getN() const127   getN() const
128   {
129     return n;
130   }
131   int
numClasses() const132   numClasses() const
133   {
134     return classes.size();
135   }
136   void trace(IntSequence &out, int n) const;
137   void
trace(IntSequence & out) const138   trace(IntSequence &out) const
139   {
140     trace(out, numClasses());
141   }
142   void trace(IntSequence &out, const Permutation &per) const;
143   void print(const std::string &prefix) const;
144   seqit
begin()145   begin()
146   {
147     return classes.begin();
148   }
149   const_seqit
begin() const150   begin() const
151   {
152     return classes.begin();
153   }
154   seqit
end()155   end()
156   {
157     return classes.end();
158   }
159   const_seqit
end() const160   end() const
161   {
162     return classes.end();
163   }
164   const_seqit find(int i) const;
165   seqit find(int i);
166 protected:
167   /* Here we have find methods. We can find an equivalence class having a
168      given number or we can find an equivalence class of a given index within
169      the ordering.
170 
171      We have also an insert() method which inserts a given class
172      according to the class ordering. */
173   const_seqit findHaving(int i) const;
174   seqit findHaving(int i);
175   void insert(const OrdSequence &s);
176 
177 };
178 
179 /* The EquivalenceSet is a list of equivalences. The unique
180    constructor constructs a set of all equivalences over an n-elements
181    set. The equivalences are sorted in the list so that equivalences with
182    fewer number of classes are in the end.
183 
184    The two methods has() and addParents() are useful in the constructor. */
185 
186 class EquivalenceSet
187 {
188   int n;
189   std::list<Equivalence> equis;
190 public:
191   explicit EquivalenceSet(int num);
192   void print(const std::string &prefix) const;
193   auto
begin() const194   begin() const
195   {
196     return equis.begin();
197   }
198   auto
end() const199   end() const
200   {
201     return equis.end();
202   }
203 private:
204   bool has(const Equivalence &e) const;
205   void addParents(const Equivalence &e, std::list<Equivalence> &added);
206 };
207 
208 /* The equivalence bundle class only encapsulates EquivalenceSet·s
209    from 1 up to a given number. It is able to retrieve the equivalence set
210    over n-element set for a given n, and also it can generate some more
211    sets on request.
212 
213    It is fully responsible for storage needed for EquivalenceSet·s. */
214 
215 class EquivalenceBundle
216 {
217   std::vector<EquivalenceSet> bundle;
218 public:
219   explicit EquivalenceBundle(int nmax);
220   const EquivalenceSet &get(int n) const;
221   void generateUpTo(int nmax);
222 };
223 
224 #endif
225