permlib
0.2.6
Library for permutation computations
|
00001 // --------------------------------------------------------------------------- 00002 // 00003 // This file is part of PermLib. 00004 // 00005 // Copyright (c) 2009-2011 Thomas Rehn <thomas@carmen76.de> 00006 // All rights reserved. 00007 // 00008 // Redistribution and use in source and binary forms, with or without 00009 // modification, are permitted provided that the following conditions 00010 // are met: 00011 // 1. Redistributions of source code must retain the above copyright 00012 // notice, this list of conditions and the following disclaimer. 00013 // 2. Redistributions in binary form must reproduce the above copyright 00014 // notice, this list of conditions and the following disclaimer in the 00015 // documentation and/or other materials provided with the distribution. 00016 // 3. The name of the author may not be used to endorse or promote products 00017 // derived from this software without specific prior written permission. 00018 // 00019 // THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR 00020 // IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES 00021 // OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. 00022 // IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT, 00023 // INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT 00024 // NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, 00025 // DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY 00026 // THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT 00027 // (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF 00028 // THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 00029 // 00030 // --------------------------------------------------------------------------- 00031 00032 #ifndef ORBIT_LEX_MIN_SEARCH_H_ 00033 #define ORBIT_LEX_MIN_SEARCH_H_ 00034 00035 #include <permlib/predicate/pointwise_stabilizer_predicate.h> 00036 #include <permlib/change/conjugating_base_change.h> 00037 #include <permlib/change/random_base_transpose.h> 00038 #include <permlib/search/dset.h> 00039 00040 #include <vector> 00041 #include <limits> 00042 #include <boost/dynamic_bitset.hpp> 00043 00044 namespace permlib { 00045 00047 00051 template<class BSGSIN> 00052 class OrbitLexMinSearch { 00053 public: 00055 00058 OrbitLexMinSearch(const BSGSIN& bsgs) 00059 : m_bsgs(bsgs), m_cbc(bsgs), m_dsetAction(bsgs.n), m_orb(m_bsgs.n), m_orbVector(m_bsgs.n, 0), m_orbVectorIndex(0) {} 00060 00062 00067 dset lexMin(const dset& element, const BSGSIN* stabilizer = NULL); 00068 00070 00075 static bool isLexSmaller(const dset& a, const dset& b); 00076 00077 private: 00078 BSGSIN m_bsgs; 00079 const BSGSIN* m_bsgsStabilizer; 00080 typedef typename BSGSIN::PERMtype PERM; 00081 typedef std::vector<boost::shared_ptr<PERM> > PERMvec; 00082 00083 struct Candidate { 00084 dset D; 00085 dset J; 00086 00087 Candidate(dset D_) : D(D_), J(D_.size()) {} 00088 Candidate(dset D_, dset J_) : D(D_), J(J_) {} 00089 00090 void print(const char* prefix) const { 00091 std::cout << prefix << ".J = " << J << " ; " << prefix << ".D = " << D << std::endl; 00092 } 00093 }; 00094 00095 typedef Candidate* CandidatePtr; 00096 00097 ConjugatingBaseChange<PERM, typename BSGSIN::TRANStype, RandomBaseTranspose<PERM, typename BSGSIN::TRANStype> > m_cbc; 00098 DSetAction<PERM> m_dsetAction; 00099 00100 bool lexMin(unsigned int i, unsigned int k, const BSGSIN* stabilizer, const std::list<CandidatePtr>& candidates, std::list<CandidatePtr>& candidatesNext, dset& M_i, std::list<unsigned long>& base, PERMvec& S_i); 00102 unsigned long orbMin(unsigned long element, const PERMvec& generators); 00103 00105 00110 dset* orbRepresentatives(dset element, const std::list<typename PERM::ptr>& generators); 00111 00112 // temporary variables for the orbMin calculation 00113 dset m_orb; 00114 std::vector<unsigned long> m_orbVector; 00115 unsigned int m_orbVectorIndex; 00116 }; 00117 00118 00119 template<class BSGSIN> 00120 inline dset OrbitLexMinSearch<BSGSIN>::lexMin(const dset& element, const BSGSIN* stabilizer) { 00121 if (element.count() == element.size()) 00122 return element; 00123 if (element.count() == 0) 00124 return element; 00125 CandidatePtr c0(new Candidate(element)); 00126 00127 std::list<CandidatePtr> candList0, candList1; 00128 std::list<CandidatePtr>* cand0 = &candList0; 00129 std::list<CandidatePtr>* cand1 = &candList1; 00130 00131 cand0->push_back(c0); 00132 dset M_i(element.size()); 00133 std::list<unsigned long> base; 00134 PERMvec S_i; 00135 S_i.reserve(m_bsgs.S.size()); 00136 00137 for (unsigned int i = 0; i < element.count(); ++i) { 00138 if (lexMin(i, element.count(), stabilizer, *cand0, *cand1, M_i, base, S_i)) 00139 break; 00140 std::swap(cand0, cand1); 00141 } 00142 std::for_each(candList0.begin(), candList0.end(), delete_object()); 00143 std::for_each(candList1.begin(), candList1.end(), delete_object()); 00144 00145 return M_i; 00146 } 00147 00148 template<class BSGSIN> 00149 inline bool OrbitLexMinSearch<BSGSIN>::lexMin(unsigned int i, unsigned int k, const BSGSIN* stabilizer, const std::list<CandidatePtr>& candidates, std::list<CandidatePtr>& candidatesNext, dset& M_i, std::list<unsigned long>& base, PERMvec& S_i) { 00150 PERMLIB_DEBUG(std::cout << "### START " << i << " with #" << candidates.size() << std::endl;) 00151 00152 // if current stabilizer in the stabilizer chain is trivial we may 00153 // choose the minimal candidate and abort the search 00154 bool allOne = true; 00155 for (unsigned int j = i; j < m_bsgs.B.size(); ++j) { 00156 if (m_bsgs.U[j].size() > 1) { 00157 allOne = false; 00158 break; 00159 } 00160 } 00161 if (allOne) { 00162 M_i = candidates.front()->D; 00163 BOOST_FOREACH(const CandidatePtr& R, candidates) { 00164 if (isLexSmaller(R->D, M_i)) { 00165 M_i = R->D; 00166 } 00167 } 00168 return true; 00169 } 00170 00171 unsigned int m = m_bsgs.n + 1; 00172 S_i.clear(); 00173 PointwiseStabilizerPredicate<PERM> stab_i(m_bsgs.B.begin(), m_bsgs.B.begin() + i); 00174 std::copy_if(m_bsgs.S.begin(), m_bsgs.S.end(), std::back_inserter(S_i), stab_i); 00175 const unsigned long UNDEFINED_ORBIT = std::numeric_limits<unsigned long>::max(); 00176 std::vector<unsigned long> orbitCache(m_bsgs.n, UNDEFINED_ORBIT); 00177 std::list<CandidatePtr> pass; 00178 00179 BOOST_FOREACH(const CandidatePtr& R, candidates) { 00180 unsigned long m_R = m; 00181 for (unsigned long j = 0; j < R->D.size(); ++j) { 00182 if (R->J[j] || !R->D[j]) 00183 continue; 00184 00185 unsigned long val = orbitCache[j]; 00186 if (val == UNDEFINED_ORBIT) { 00187 val = orbMin(j, S_i); 00188 orbitCache[j] = val; 00189 } 00190 if (m_R > val) 00191 m_R = val; 00192 } 00193 00194 if (m_R < m) { 00195 m = m_R; 00196 pass.clear(); 00197 pass.push_back(R); 00198 } else if (m_R == m) { 00199 pass.push_back(R); 00200 } 00201 } 00202 00203 PERMLIB_DEBUG(std::cout << " found m = " << m << std::endl;) 00204 M_i.set(m, 1); 00205 if (i == k-1) 00206 return true; 00207 00208 base.push_back(m); 00209 m_cbc.change(m_bsgs, base.begin(), base.end()); 00210 00211 std::for_each(candidatesNext.begin(), candidatesNext.end(), delete_object()); 00212 candidatesNext.clear(); 00213 00214 PERM* UNDEFINED_TRANSVERSAL = reinterpret_cast<PERM*>(1L); 00215 std::vector<PERM*> transversalCache(m_bsgs.n); 00216 BOOST_FOREACH(PERM*& p, transversalCache) { 00217 p = UNDEFINED_TRANSVERSAL; 00218 } 00219 BOOST_FOREACH(const CandidatePtr& R, pass) { 00220 for (unsigned long j = 0; j < R->D.size(); ++j) { 00221 if (!R->D[j]) 00222 continue; 00223 00224 PERM* perm = transversalCache[j]; 00225 if (perm == UNDEFINED_TRANSVERSAL) { 00226 perm = m_bsgs.U[i].at(j); 00227 if (perm) { 00228 perm->invertInplace(); 00229 } 00230 transversalCache[j] = perm; 00231 } 00232 00233 if (!perm) 00234 continue; 00235 00236 CandidatePtr c(new Candidate(R->D, R->J)); 00237 m_dsetAction.apply(*perm, R->D, c->D); 00238 c->J.set(m); 00239 candidatesNext.push_back(c); 00240 } 00241 } 00242 00243 BOOST_FOREACH(PERM* p, transversalCache) { 00244 if (p != UNDEFINED_TRANSVERSAL) 00245 delete p; 00246 } 00247 return false; 00248 } 00249 00250 template<class BSGSIN> 00251 inline unsigned long OrbitLexMinSearch<BSGSIN>::orbMin(unsigned long element, const PERMvec& generators) { 00252 if (element == 0) 00253 return 0; 00254 00255 unsigned long minElement = element; 00256 m_orb.reset(); 00257 m_orb.set(element, 1); 00258 m_orbVectorIndex = 0; 00259 m_orbVector[m_orbVectorIndex++] = element; 00260 00261 for (unsigned int i = 0; i < m_orbVectorIndex; ++i) { 00262 const unsigned long &alpha = m_orbVector[i]; 00263 BOOST_FOREACH(const typename PERM::ptr& p, generators) { 00264 unsigned long alpha_p = *p / alpha; 00265 if (alpha_p == 0) 00266 return 0; 00267 if (!m_orb[alpha_p]) { 00268 m_orbVector[m_orbVectorIndex++] = alpha_p; 00269 m_orb.set(alpha_p); 00270 if (alpha_p < minElement) 00271 minElement = alpha_p; 00272 } 00273 } 00274 } 00275 00276 return minElement; 00277 } 00278 00279 00280 template<class BSGSIN> 00281 inline dset* OrbitLexMinSearch<BSGSIN>::orbRepresentatives(dset element, const std::list<typename PERM::ptr>& generators) { 00282 dset* ret = new dset(element.size()); 00283 00284 for (unsigned int j = 0; j < element.size(); ++j) { 00285 if (!element[j]) 00286 continue; 00287 00288 m_orb.set(); 00289 m_orb.set(j, 0); 00290 m_orbVectorIndex = 0; 00291 m_orbVector[m_orbVectorIndex++] = j; 00292 for (unsigned int i = 0; i < m_orbVectorIndex; ++i) { 00293 const unsigned long &alpha = m_orbVector[i]; 00294 BOOST_FOREACH(const typename PERM::ptr& p, generators) { 00295 unsigned long alpha_p = *p / alpha; 00296 if (m_orb[alpha_p]) { 00297 m_orbVector[m_orbVectorIndex++] = alpha_p; 00298 m_orb.reset(alpha_p); 00299 } 00300 } 00301 } 00302 00303 element &= m_orb; 00304 ret->set(j); 00305 } 00306 00307 return ret; 00308 } 00309 00310 00311 template<class BSGSIN> 00312 inline bool OrbitLexMinSearch<BSGSIN>::isLexSmaller(const dset& a, const dset& b) { 00313 dset::size_type i = a.find_first(), j = b.find_first(); 00314 while (i != dset::npos && j != dset::npos) { 00315 if (i < j) 00316 return true; 00317 if (i > j) 00318 return false; 00319 i = a.find_next(i); 00320 j = b.find_next(j); 00321 } 00322 return false; 00323 } 00324 00325 } // ns permlib 00326 00327 #endif // ORBIT_LEX_MIN_SEARCH_H_