permlib  0.2.9
Library for permutation computations
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 
44 namespace permlib {
45 
47 
51 template<class BSGSIN>
53 public:
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 
80 private:
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 
123 template<class BSGSIN>
124 inline 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 
152 template<class BSGSIN>
153 inline 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 
259 template<class BSGSIN>
260 inline 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 
289 template<class BSGSIN>
290 inline 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 
320 template<class BSGSIN>
321 inline 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_
dset lexMin(const dset &element, const BSGSIN *stabilizer=NULL)
searches the lexicographically minimal element of an orbit
Definition: orbit_lex_min_search.h:124
OrbitLexMinSearch(const BSGSIN &bsgs, bool abortSearchIfSmallerFound=false)
constructor
Definition: orbit_lex_min_search.h:60
static bool isLexSmaller(const dset &a, const dset &b)
compares two sets given as dsets lexicographically
Definition: orbit_lex_min_search.h:321
algorithm to find the lexicographically minimal set in an orbit
Definition: orbit_lex_min_search.h:52
callable object to delete a pointer
Definition: common.h:82
Definition: abstract_bsgs.h:49