permlib  0.2.6
Library for permutation computations
 All Classes Functions Variables Typedefs Enumerations Friends
include/permlib/search/orbit_lex_min_search.h
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_