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