32 #ifndef ORBIT_LEX_MIN_SEARCH_H_ 33 #define ORBIT_LEX_MIN_SEARCH_H_ 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> 42 #include <boost/dynamic_bitset.hpp> 51 template<
class BSGSIN>
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) {}
70 dset
lexMin(
const dset& element,
const BSGSIN* stabilizer = NULL);
82 const BSGSIN* m_bsgsStabilizer;
83 typedef typename BSGSIN::PERMtype PERM;
84 typedef std::vector<boost::shared_ptr<PERM> > PERMvec;
90 Candidate(dset D_) : D(D_), J(D_.size()) {}
91 Candidate(dset D_, dset J_) : D(D_), J(J_) {}
93 void print(
const char* prefix)
const {
94 std::cout << prefix <<
".J = " << J <<
" ; " << prefix <<
".D = " << D << std::endl;
98 typedef Candidate* CandidatePtr;
100 ConjugatingBaseChange<PERM, typename BSGSIN::TRANStype, RandomBaseTranspose<PERM, typename BSGSIN::TRANStype> > m_cbc;
101 DSetAction<PERM> m_dsetAction;
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);
113 dset* orbRepresentatives(dset element,
const std::list<typename PERM::ptr>& generators);
117 std::vector<unsigned long> m_orbVector;
118 unsigned int m_orbVectorIndex;
119 const bool m_abortSearchIfSmallerFound;
123 template<
class BSGSIN>
125 if (element.count() == element.size())
127 if (element.count() == 0)
129 CandidatePtr c0(
new Candidate(element));
131 std::list<CandidatePtr> candList0, candList1;
132 std::list<CandidatePtr>* cand0 = &candList0;
133 std::list<CandidatePtr>* cand1 = &candList1;
135 cand0->push_back(c0);
136 dset M_i(element.size());
137 std::list<unsigned long> base;
139 S_i.reserve(m_bsgs.S.size());
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))
144 std::swap(cand0, cand1);
146 std::for_each(candList0.begin(), candList0.end(),
delete_object());
147 std::for_each(candList1.begin(), candList1.end(),
delete_object());
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;)
159 for (
unsigned int j = i; j < m_bsgs.B.size(); ++j) {
160 if (m_bsgs.U[j].size() > 1) {
166 M_i = candidates.front()->D;
167 BOOST_FOREACH(
const CandidatePtr& R, candidates) {
168 if (isLexSmaller(R->D, M_i)) {
175 unsigned int m = m_bsgs.n + 1;
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;
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());
190 for (
unsigned long j = 0; j < R->D.size(); ++j) {
191 if (R->J[j] || !R->D[j])
194 unsigned long val = orbitCache[j];
195 if (val == UNDEFINED_ORBIT) {
196 val = orbMin(j, S_i);
207 }
else if (m_R == m) {
212 PERMLIB_DEBUG(std::cout <<
" found m = " << m << std::endl;)
218 m_cbc.change(m_bsgs, base.begin(), base.end());
220 std::for_each(candidatesNext.begin(), candidatesNext.end(), delete_object());
221 candidatesNext.clear();
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;
228 BOOST_FOREACH(
const CandidatePtr& R, pass) {
229 for (
unsigned long j = 0; j < R->D.size(); ++j) {
233 PERM* perm = transversalCache[j];
234 if (perm == UNDEFINED_TRANSVERSAL) {
235 perm = m_bsgs.U[i].at(j);
237 perm->invertInplace();
239 transversalCache[j] = perm;
245 CandidatePtr c(
new Candidate(R->D, R->J));
246 m_dsetAction.apply(*perm, R->D, c->D);
248 candidatesNext.push_back(c);
252 BOOST_FOREACH(PERM* p, transversalCache) {
253 if (p != UNDEFINED_TRANSVERSAL)
259 template<
class BSGSIN>
260 inline unsigned long OrbitLexMinSearch<BSGSIN>::orbMin(
unsigned long element,
const PERMvec& generators) {
264 unsigned long minElement = element;
266 m_orb.set(element, 1);
267 m_orbVectorIndex = 0;
268 m_orbVector[m_orbVectorIndex++] = element;
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;
276 if (!m_orb[alpha_p]) {
277 m_orbVector[m_orbVectorIndex++] = alpha_p;
279 if (alpha_p < minElement)
280 minElement = alpha_p;
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());
293 for (
unsigned int j = 0; j < element.size(); ++j) {
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);
320 template<
class BSGSIN>
322 dset::size_type i = a.find_first(), j = b.find_first();
323 while (i != dset::npos && j != dset::npos) {
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