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 // Permutations.
22 
23 /* The permutation class is useful when describing a permutation of
24    indices in permuted symmetry tensor. This tensor comes to existence,
25    for instance, as a result of the following tensor multiplication:
26 
27     [g_y³]_γ₁γ₂γ₃ [g_yu]^γ₁_α₁β₃ [g_yu]^γ₂_α₂β₁ [g_u]^γ₃_β₂
28 
29    If this operation is done by a Kronecker product of unfolded tensors, the
30    resulting tensor has permuted indices. So, in this case the permutation is
31    implied by the equivalence: { {0,4}, {1,3}, {2} }. This results in a
32    permutation which maps indices (0,1,2,3,4)↦(0,2,4,3,1).
33 
34    The other application of Permutation class is to permute indices with the
35    same permutation as done during sorting.
36 
37    Here we only define an abstraction for the permutation defined by an
38    equivalence. Its basic operation is to apply the permutation to the integer
39    sequence. The application is right (or inner), in sense that it works on
40    indices of the sequence not items of the sequence. More formally s∘m ≠ m∘s.
41    In here, the application of the permutation defined by map m is s∘m.
42 
43    Also, we need PermutationSet class which contains all permutations
44    of an n-element set, and a bundle of permutations PermutationBundle
45    which contains all permutation sets up to a given number.
46 */
47 
48 #ifndef PERMUTATION_H
49 #define PERMUTATION_H
50 
51 #include "int_sequence.hh"
52 #include "equivalence.hh"
53 
54 #include <vector>
55 
56 /* The permutation object will have a map, which defines mapping of indices
57    (0,1,…,n-1)↦(m₀,m₁,…, mₙ₋₁). The map is the sequence (m₀,m₁,…, mₙ₋₁). When
58    the permutation with the map m is applied on sequence s, it permutes its
59    indices: s∘id↦s∘m.
60 
61    So we have one constructor from equivalence, then a method apply(),
62    and finally a method tailIdentity() which returns a number of trailing
63    indices which yield identity. Also we have a constructor calculating
64    map, which corresponds to permutation in sort. This is, we want
65    (sorted s)∘m = s.
66 */
67 
68 class Permutation
69 {
70 protected:
71   IntSequence permap;
72 public:
Permutation(int len)73   explicit Permutation(int len)
74     : permap(len)
75   {
76     for (int i = 0; i < len; i++)
77       permap[i] = i;
78   }
Permutation(const Equivalence & e)79   explicit Permutation(const Equivalence &e)
80     : permap(e.getN())
81   {
82     e.trace(permap);
83   }
Permutation(const Equivalence & e,const Permutation & per)84   Permutation(const Equivalence &e, const Permutation &per)
85     : permap(e.getN())
86   {
87     e.trace(permap, per);
88   }
Permutation(const IntSequence & s)89   explicit Permutation(const IntSequence &s)
90     : permap(s.size())
91   {
92     computeSortingMap(s);
93   };
Permutation(const Permutation & p1,const Permutation & p2)94   Permutation(const Permutation &p1, const Permutation &p2)
95     : permap(p2.permap)
96   {
97     p1.apply(permap);
98   }
Permutation(const Permutation & p,int i)99   Permutation(const Permutation &p, int i)
100     : permap(p.permap.insert(p.size(), i))
101   {
102   }
103   bool
operator ==(const Permutation & p)104   operator==(const Permutation &p)
105   {
106     return permap == p.permap;
107   }
108   int
size() const109   size() const
110   {
111     return permap.size();
112   }
113   void
print() const114   print() const
115   {
116     permap.print();
117   }
118   void apply(const IntSequence &src, IntSequence &tar) const;
119   void apply(IntSequence &tar) const;
120   void inverse();
121   int tailIdentity() const;
122   const IntSequence &
getMap() const123   getMap() const
124   {
125     return permap;
126   }
127   IntSequence &
getMap()128   getMap()
129   {
130     return permap;
131   }
132 protected:
133   void computeSortingMap(const IntSequence &s);
134 };
135 
136 /* The PermutationSet maintains an array of of all permutations. The default
137    constructor constructs one element permutation set of one element sets. The
138    second constructor constructs a new permutation set over n from all
139    permutations over n-1. The parameter n needs not to be provided, but it
140    serves to distinguish the constructor from copy constructor.
141 
142    The method getPreserving() returns a factor subgroup of permutations, which
143    are invariants with respect to the given sequence. This are all permutations
144    p yielding p∘s = s, where s is the given sequence. */
145 
146 class PermutationSet
147 {
148   int order{1};
149   int size{1};
150   std::vector<Permutation> pers;
151 public:
152   PermutationSet();
153   PermutationSet(const PermutationSet &ps, int n);
154   int
getNum() const155   getNum() const
156   {
157     return size;
158   }
159   const Permutation &
get(int i) const160   get(int i) const
161   {
162     return pers[i];
163   }
164   std::vector<Permutation> getPreserving(const IntSequence &s) const;
165 };
166 
167 /* The permutation bundle encapsulates all permutations sets up to some
168    given dimension. */
169 
170 class PermutationBundle
171 {
172   std::vector<PermutationSet> bundle;
173 public:
174   PermutationBundle(int nmax);
175   const PermutationSet &get(int n) const;
176   void generateUpTo(int nmax);
177 };
178 
179 #endif
180