permlib 0.2.9
Library for permutation computations
Loading...
Searching...
No Matches
orbit_lex_min_search.h
1// ---------------------------------------------------------------------------
2//
3// This file is part of PermLib.
4//
5// Copyright (c) 2009-2011 Thomas Rehn <thomas@carmen76.de>
6// All rights reserved.
7//
8// Redistribution and use in source and binary forms, with or without
9// modification, are permitted provided that the following conditions
10// are met:
11// 1. Redistributions of source code must retain the above copyright
12// notice, this list of conditions and the following disclaimer.
13// 2. Redistributions in binary form must reproduce the above copyright
14// notice, this list of conditions and the following disclaimer in the
15// documentation and/or other materials provided with the distribution.
16// 3. The name of the author may not be used to endorse or promote products
17// derived from this software without specific prior written permission.
18//
19// THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
20// IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
21// OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
22// IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
23// INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
24// NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
25// DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
26// THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
27// (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
28// THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
29//
30// ---------------------------------------------------------------------------
31
32#ifndef ORBIT_LEX_MIN_SEARCH_H_
33#define ORBIT_LEX_MIN_SEARCH_H_
34
35#include <permlib/predicate/pointwise_stabilizer_predicate.h>
36#include <permlib/change/conjugating_base_change.h>
37#include <permlib/change/random_base_transpose.h>
38#include <permlib/search/dset.h>
39
40#include <vector>
41#include <limits>
42#include <boost/dynamic_bitset.hpp>
43
44namespace permlib {
45
47
51template<class BSGSIN>
53public:
55
60 OrbitLexMinSearch(const BSGSIN& bsgs, bool abortSearchIfSmallerFound = false)
61 : m_bsgs(bsgs), m_cbc(bsgs), m_dsetAction(bsgs.n), m_orb(m_bsgs.n), m_orbVector(m_bsgs.n, 0),
62 m_orbVectorIndex(0), m_abortSearchIfSmallerFound(abortSearchIfSmallerFound) {}
63
65
70 dset lexMin(const dset& element, const BSGSIN* stabilizer = NULL);
71
73
78 static bool isLexSmaller(const dset& a, const dset& b);
79
80private:
81 BSGSIN m_bsgs;
82 const BSGSIN* m_bsgsStabilizer;
83 typedef typename BSGSIN::PERMtype PERM;
84 typedef std::vector<boost::shared_ptr<PERM> > PERMvec;
85
86 struct Candidate {
87 dset D;
88 dset J;
89
90 Candidate(dset D_) : D(D_), J(D_.size()) {}
91 Candidate(dset D_, dset J_) : D(D_), J(J_) {}
92
93 void print(const char* prefix) const {
94 std::cout << prefix << ".J = " << J << " ; " << prefix << ".D = " << D << std::endl;
95 }
96 };
97
98 typedef Candidate* CandidatePtr;
99
100 ConjugatingBaseChange<PERM, typename BSGSIN::TRANStype, RandomBaseTranspose<PERM, typename BSGSIN::TRANStype> > m_cbc;
101 DSetAction<PERM> m_dsetAction;
102
103 bool lexMin(unsigned int i, unsigned int k, const dset& element, const BSGSIN* stabilizer, const std::list<CandidatePtr>& candidates, std::list<CandidatePtr>& candidatesNext, dset& M_i, std::list<unsigned long>& base, PERMvec& S_i);
105 unsigned long orbMin(unsigned long element, const PERMvec& generators);
106
108
113 dset* orbRepresentatives(dset element, const std::list<typename PERM::ptr>& generators);
114
115 // temporary variables for the orbMin calculation
116 dset m_orb;
117 std::vector<unsigned long> m_orbVector;
118 unsigned int m_orbVectorIndex;
119 const bool m_abortSearchIfSmallerFound;
120};
121
122
123template<class BSGSIN>
124inline dset OrbitLexMinSearch<BSGSIN>::lexMin(const dset& element, const BSGSIN* stabilizer) {
125 if (element.count() == element.size())
126 return element;
127 if (element.count() == 0)
128 return element;
129 CandidatePtr c0(new Candidate(element));
130
131 std::list<CandidatePtr> candList0, candList1;
132 std::list<CandidatePtr>* cand0 = &candList0;
133 std::list<CandidatePtr>* cand1 = &candList1;
134
135 cand0->push_back(c0);
136 dset M_i(element.size());
137 std::list<unsigned long> base;
138 PERMvec S_i;
139 S_i.reserve(m_bsgs.S.size());
140
141 for (unsigned int i = 0; i < element.count(); ++i) {
142 if (lexMin(i, element.count(), element, stabilizer, *cand0, *cand1, M_i, base, S_i))
143 break;
144 std::swap(cand0, cand1);
145 }
146 std::for_each(candList0.begin(), candList0.end(), delete_object());
147 std::for_each(candList1.begin(), candList1.end(), delete_object());
148
149 return M_i;
150}
151
152template<class BSGSIN>
153inline bool OrbitLexMinSearch<BSGSIN>::lexMin(unsigned int i, unsigned int k, const dset& element, const BSGSIN* stabilizer, const std::list<CandidatePtr>& candidates, std::list<CandidatePtr>& candidatesNext, dset& M_i, std::list<unsigned long>& base, PERMvec& S_i) {
154 PERMLIB_DEBUG(std::cout << "### START " << i << " with #" << candidates.size() << std::endl;)
155
156 // if current stabilizer in the stabilizer chain is trivial we may
157 // choose the minimal candidate and abort the search
158 bool allOne = true;
159 for (unsigned int j = i; j < m_bsgs.B.size(); ++j) {
160 if (m_bsgs.U[j].size() > 1) {
161 allOne = false;
162 break;
163 }
164 }
165 if (allOne) {
166 M_i = candidates.front()->D;
167 BOOST_FOREACH(const CandidatePtr& R, candidates) {
168 if (isLexSmaller(R->D, M_i)) {
169 M_i = R->D;
170 }
171 }
172 return true;
173 }
174
175 unsigned int m = m_bsgs.n + 1;
176 S_i.clear();
177 PointwiseStabilizerPredicate<PERM> stab_i(m_bsgs.B.begin(), m_bsgs.B.begin() + i);
178 std::copy_if(m_bsgs.S.begin(), m_bsgs.S.end(), std::back_inserter(S_i), stab_i);
179 const unsigned long UNDEFINED_ORBIT = std::numeric_limits<unsigned long>::max();
180 std::vector<unsigned long> orbitCache(m_bsgs.n, UNDEFINED_ORBIT);
181 std::list<CandidatePtr> pass;
182
183 BOOST_FOREACH(const CandidatePtr& R, candidates) {
184 unsigned long m_R = m;
185 if (m_abortSearchIfSmallerFound && isLexSmaller(R->D, element)) {
186 BOOST_ASSERT(R->D.count() == element.count());
187 M_i = R->D;
188 return true;
189 }
190 for (unsigned long j = 0; j < R->D.size(); ++j) {
191 if (R->J[j] || !R->D[j])
192 continue;
193
194 unsigned long val = orbitCache[j];
195 if (val == UNDEFINED_ORBIT) {
196 val = orbMin(j, S_i);
197 orbitCache[j] = val;
198 }
199 if (m_R > val)
200 m_R = val;
201 }
202
203 if (m_R < m) {
204 m = m_R;
205 pass.clear();
206 pass.push_back(R);
207 } else if (m_R == m) {
208 pass.push_back(R);
209 }
210 }
211
212 PERMLIB_DEBUG(std::cout << " found m = " << m << std::endl;)
213 M_i.set(m, 1);
214 if (i == k-1)
215 return true;
216
217 base.push_back(m);
218 m_cbc.change(m_bsgs, base.begin(), base.end());
219
220 std::for_each(candidatesNext.begin(), candidatesNext.end(), delete_object());
221 candidatesNext.clear();
222
223 PERM* UNDEFINED_TRANSVERSAL = reinterpret_cast<PERM*>(1L);
224 std::vector<PERM*> transversalCache(m_bsgs.n);
225 BOOST_FOREACH(PERM*& p, transversalCache) {
226 p = UNDEFINED_TRANSVERSAL;
227 }
228 BOOST_FOREACH(const CandidatePtr& R, pass) {
229 for (unsigned long j = 0; j < R->D.size(); ++j) {
230 if (!R->D[j])
231 continue;
232
233 PERM* perm = transversalCache[j];
234 if (perm == UNDEFINED_TRANSVERSAL) {
235 perm = m_bsgs.U[i].at(j);
236 if (perm) {
237 perm->invertInplace();
238 }
239 transversalCache[j] = perm;
240 }
241
242 if (!perm)
243 continue;
244
245 CandidatePtr c(new Candidate(R->D, R->J));
246 m_dsetAction.apply(*perm, R->D, c->D);
247 c->J.set(m);
248 candidatesNext.push_back(c);
249 }
250 }
251
252 BOOST_FOREACH(PERM* p, transversalCache) {
253 if (p != UNDEFINED_TRANSVERSAL)
254 delete p;
255 }
256 return false;
257}
258
259template<class BSGSIN>
260inline unsigned long OrbitLexMinSearch<BSGSIN>::orbMin(unsigned long element, const PERMvec& generators) {
261 if (element == 0)
262 return 0;
263
264 unsigned long minElement = element;
265 m_orb.reset();
266 m_orb.set(element, 1);
267 m_orbVectorIndex = 0;
268 m_orbVector[m_orbVectorIndex++] = element;
269
270 for (unsigned int i = 0; i < m_orbVectorIndex; ++i) {
271 const unsigned long &alpha = m_orbVector[i];
272 BOOST_FOREACH(const typename PERM::ptr& p, generators) {
273 unsigned long alpha_p = *p / alpha;
274 if (alpha_p == 0)
275 return 0;
276 if (!m_orb[alpha_p]) {
277 m_orbVector[m_orbVectorIndex++] = alpha_p;
278 m_orb.set(alpha_p);
279 if (alpha_p < minElement)
280 minElement = alpha_p;
281 }
282 }
283 }
284
285 return minElement;
286}
287
288
289template<class BSGSIN>
290inline dset* OrbitLexMinSearch<BSGSIN>::orbRepresentatives(dset element, const std::list<typename PERM::ptr>& generators) {
291 dset* ret = new dset(element.size());
292
293 for (unsigned int j = 0; j < element.size(); ++j) {
294 if (!element[j])
295 continue;
296
297 m_orb.set();
298 m_orb.set(j, 0);
299 m_orbVectorIndex = 0;
300 m_orbVector[m_orbVectorIndex++] = j;
301 for (unsigned int i = 0; i < m_orbVectorIndex; ++i) {
302 const unsigned long &alpha = m_orbVector[i];
303 BOOST_FOREACH(const typename PERM::ptr& p, generators) {
304 unsigned long alpha_p = *p / alpha;
305 if (m_orb[alpha_p]) {
306 m_orbVector[m_orbVectorIndex++] = alpha_p;
307 m_orb.reset(alpha_p);
308 }
309 }
310 }
311
312 element &= m_orb;
313 ret->set(j);
314 }
315
316 return ret;
317}
318
319
320template<class BSGSIN>
321inline bool OrbitLexMinSearch<BSGSIN>::isLexSmaller(const dset& a, const dset& b) {
322 dset::size_type i = a.find_first(), j = b.find_first();
323 while (i != dset::npos && j != dset::npos) {
324 if (i < j)
325 return true;
326 if (i > j)
327 return false;
328 i = a.find_next(i);
329 j = b.find_next(j);
330 }
331 return false;
332 }
333
334} // ns permlib
335
336#endif // ORBIT_LEX_MIN_SEARCH_H_
algorithm to find the lexicographically minimal set in an orbit
Definition: orbit_lex_min_search.h:52
OrbitLexMinSearch(const BSGSIN &bsgs, bool abortSearchIfSmallerFound=false)
constructor
Definition: orbit_lex_min_search.h:60
dset lexMin(const dset &element, const BSGSIN *stabilizer=NULL)
searches the lexicographically minimal element of an orbit
Definition: orbit_lex_min_search.h:124
static bool isLexSmaller(const dset &a, const dset &b)
compares two sets given as dsets lexicographically
Definition: orbit_lex_min_search.h:321
OutputIterator copy_if(InputIterator begin, InputIterator end, OutputIterator destBegin, Predicate p)
copies elements of (begin to end) to destBegin if they match the given predicate
Definition: common.h:49
callable object to delete a pointer
Definition: common.h:82