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 00033 #ifndef RBASE_H_ 00034 #define RBASE_H_ 00035 00036 #include <permlib/predicate/subgroup_predicate.h> 00037 00038 #include <permlib/search/base_search.h> 00039 00040 #include <permlib/search/partition/partition.h> 00041 #include <permlib/search/partition/refinement_family.h> 00042 #include <permlib/search/partition/backtrack_refinement.h> 00043 00044 #include <permlib/change/conjugating_base_change.h> 00045 #include <permlib/change/random_base_transpose.h> 00046 00047 #include <permlib/sorter/base_sorter.h> 00048 00049 #include <utility> 00050 00051 namespace permlib { 00052 namespace partition { 00053 00055 template<class BSGSIN,class TRANSRET> 00056 class RBase : public BaseSearch<BSGSIN,TRANSRET> { 00057 public: 00058 typedef typename BaseSearch<BSGSIN,TRANSRET>::PERM PERM; 00059 typedef typename BaseSearch<BSGSIN,TRANSRET>::TRANS TRANS; 00060 00062 00067 RBase(const BSGSIN& bsgs, unsigned int pruningLevelDCM, bool stopAfterFirstElement = false); 00068 00069 typedef typename Refinement<PERM>::RefinementPtr RefinementPtr; 00070 typedef typename RefinementFamily<PERM>::PartitionPtr PartitionPtr; 00071 typedef typename std::list<std::pair<PartitionPtr,RefinementPtr> >::const_iterator PartitionIt; 00072 00074 void search(BSGS<PERM,TRANSRET> &groupK); 00075 00076 using BaseSearch<BSGSIN,TRANSRET>::searchCosetRepresentative; 00077 virtual typename BaseSearch<BSGSIN,TRANSRET>::PERM::ptr searchCosetRepresentative(BSGS<PERM,TRANSRET> &groupK, BSGS<PERM,TRANSRET> &groupL); 00078 protected: 00080 Partition m_partition; 00081 Partition m_partition2; 00082 00084 00089 void construct(SubgroupPredicate<PERM>* pred, RefinementFamily<PERM>* predRefinement); 00090 00092 virtual unsigned int processNewFixPoints(const Partition& pi, unsigned int level); 00093 00094 virtual const std::vector<dom_int>& subgroupBase() const; 00095 private: 00097 std::vector<dom_int> m_subgroupBase; 00099 std::list<std::pair<PartitionPtr,RefinementPtr> > partitions; 00100 00102 unsigned int search(PartitionIt pIt, Partition &pi, const PERM& t, const PERM* t2, unsigned int level, unsigned int backtrackLevel, unsigned int& completed, BSGS<PERM,TRANSRET> &groupK, BSGS<PERM,TRANSRET> &groupL); 00103 00105 bool updateMappingPermutation(const BSGSIN& bsgs, const Partition& sigma, const Partition& pi2, PERM& t2) const; 00106 }; 00107 00108 template<class BSGSIN,class TRANSRET> 00109 RBase<BSGSIN,TRANSRET>::RBase(const BSGSIN& bsgs, unsigned int pruningLevelDCM, bool stopAfterFirstElement) 00110 : BaseSearch<BSGSIN,TRANSRET>(bsgs, pruningLevelDCM, stopAfterFirstElement), 00111 m_partition(bsgs.n), m_partition2(bsgs.n) 00112 { } 00113 00114 template<class BSGSIN,class TRANSRET> 00115 void RBase<BSGSIN,TRANSRET>::construct(SubgroupPredicate<PERM>* pred, RefinementFamily<PERM>* predRefinement) { 00116 this->m_pred.reset(pred); 00117 typedef typename boost::shared_ptr<RefinementFamily<PERM> > RefinementFamilyPtr; 00118 std::list<RefinementFamilyPtr> refinements; 00119 00120 if (!this->m_bsgs.isSymmetricGroup()) { 00121 RefinementFamilyPtr gr( new GroupRefinementFamily<PERM,TRANS>(this->m_bsgs) ); 00122 refinements.push_back(gr); 00123 } 00124 00125 if (predRefinement) { 00126 RefinementFamilyPtr predR( predRefinement ); 00127 refinements.push_back(predR); 00128 } 00129 00130 PERMLIB_DEBUG(print_iterable(this->m_bsgs.B.begin(), this->m_bsgs.B.end(), 1, "orig BSGS");) 00131 00132 Partition pi(m_partition); 00133 while (pi.cells() < this->m_bsgs.n) { 00134 PERMLIB_DEBUG(std::cout << std::endl << "PI0 = " << pi << std::endl;) 00135 bool found = false; 00136 do { 00137 found = false; 00138 unsigned int foo = 0; 00139 BOOST_FOREACH(RefinementFamilyPtr ref, refinements) { 00140 const unsigned int oldFixPointsSize = pi.fixPointsSize(); 00141 std::pair<PartitionPtr,RefinementPtr> newRef = ref->apply(pi); 00142 if (newRef.first) { 00143 partitions.push_back(newRef); 00144 if (oldFixPointsSize < pi.fixPointsSize()) { 00145 processNewFixPoints(pi, partitions.size()); 00146 } 00147 //std::cout << "BSGS " << this->m_bsgs; 00148 found = true; 00149 } 00150 ++foo; 00151 } 00152 } while(found); 00153 00154 PERMLIB_DEBUG(std::cout << std::endl << "PI1 = " << pi << std::endl;) 00155 00156 if (pi.cells() < this->m_bsgs.n) { 00157 unsigned int alpha = -1; 00158 //print_iterable(pi.fixPointsBegin(), pi.fixPointsEnd(), 1, " fix0"); 00159 //print_iterable(this->m_bsgs.B.begin(), this->m_bsgs.B.end(), 1, "bsgs0"); 00160 if (pi.fixPointsSize() < this->m_bsgs.B.size()) 00161 alpha = this->m_bsgs.B[pi.fixPointsSize()]; 00162 if (alpha >= this->m_bsgs.n) { 00163 for (unsigned int i = 0; i < this->m_bsgs.n; ++i) { 00164 if (std::find(pi.fixPointsBegin(), pi.fixPointsEnd(), i) == pi.fixPointsEnd()) { 00165 alpha = i; 00166 break; 00167 } 00168 } 00169 } 00170 BOOST_ASSERT( alpha < this->m_bsgs.n ); 00171 PERMLIB_DEBUG(std::cout << "choose alpha = " << alpha << std::endl;) 00172 RefinementPtr br(new BacktrackRefinement<PERM>(this->m_bsgs.n, alpha)); 00173 BacktrackRefinement<PERM>* ref = dynamic_cast<BacktrackRefinement<PERM> *>(br.get()); 00174 ref->initializeAndApply(pi); 00175 PartitionPtr newPi(new Partition(pi)); 00176 PERMLIB_DEBUG(std::cout << "BACKTRACK " << (ref->alpha()+1) << " in " << pi << " --> " << *newPi << std::endl;) 00177 partitions.push_back(std::make_pair(newPi, br)); 00178 00179 processNewFixPoints(pi, partitions.size()); 00180 00181 //std::cout << "BSGS " << this->m_bsgs; 00182 m_subgroupBase.push_back(ref->alpha()); 00183 } 00184 } 00185 00186 this->m_order = BaseSorterByReference::createOrder(this->m_bsgs.n, pi.fixPointsBegin(), pi.fixPointsEnd()); 00187 this->m_sorter.reset(new BaseSorterByReference(this->m_order)); 00188 for (typename std::list<std::pair<PartitionPtr,RefinementPtr> >::iterator pIt = partitions.begin(); pIt != partitions.end(); ++pIt) { 00189 (*pIt).second->sort(*this->m_sorter, 0); 00190 PERMLIB_DEBUG(std::cout << "SIGMA = " << *(*pIt).first << std::endl;) 00191 } 00192 00193 PERMLIB_DEBUG(print_iterable(this->m_order.begin(), this->m_order.end(), 0, "ORDER");) 00194 } 00195 00196 template<class BSGSIN,class TRANSRET> 00197 unsigned int RBase<BSGSIN,TRANSRET>::processNewFixPoints(const Partition& pi, unsigned int level) { 00198 const unsigned int basePos = this->m_baseChange.change(this->m_bsgs, pi.fixPointsBegin(), pi.fixPointsEnd(), true); 00199 if (this->m_bsgs2) 00200 this->m_baseChange.change(*this->m_bsgs2, pi.fixPointsBegin(), pi.fixPointsEnd(), true); 00201 //print_iterable(pi.fixPointsBegin(), pi.fixPointsEnd(), 1, " fix"); 00202 PERMLIB_DEBUG(print_iterable(this->m_bsgs.B.begin(), this->m_bsgs.B.end(), 1, "change base");) 00203 return basePos; 00204 } 00205 00206 template<class BSGSIN,class TRANSRET> 00207 void RBase<BSGSIN,TRANSRET>::search(BSGS<PERM,TRANSRET> &groupK) { 00208 BOOST_ASSERT( this->m_pred != 0 ); 00209 00210 this->setupEmptySubgroup(groupK); 00211 00212 unsigned int completed = partitions.size(); 00213 BSGS<PERM,TRANSRET> groupL(groupK); 00214 PERM identH(this->m_bsgs.n); 00215 search(partitions.begin(), m_partition2, PERM(this->m_bsgs.n), &identH, 0, 0, completed, groupK, groupL); 00216 } 00217 00218 template<class BSGSIN,class TRANSRET> 00219 typename BaseSearch<BSGSIN,TRANSRET>::PERM::ptr RBase<BSGSIN,TRANSRET>::searchCosetRepresentative(BSGS<PERM,TRANSRET> &groupK, BSGS<PERM,TRANSRET> &groupL) { 00220 BOOST_ASSERT( this->m_pred != 0 ); 00221 00222 // !!! 00223 // 00224 // TODO: check that groupK and groupL have the right base (starting with subgroupBase) 00225 // 00226 // !!! 00227 00228 unsigned int completed = partitions.size(); 00229 //BSGS<PERM,TRANS> groupL(groupK); 00230 PERM identH(this->m_bsgs.n); 00231 search(partitions.begin(), m_partition2, PERM(this->m_bsgs.n), &identH, 0, 0, completed, groupK, groupL); 00232 00233 return BaseSearch<BSGSIN,TRANSRET>::m_lastElement; 00234 } 00235 00236 00237 00238 template<class BSGSIN,class TRANSRET> 00239 unsigned int RBase<BSGSIN,TRANSRET>::search(PartitionIt pIt, Partition &pi, const PERM& t, const PERM* t2, unsigned int level, unsigned int backtrackLevel, unsigned int& completed, BSGS<PERM,TRANSRET> &groupK, BSGS<PERM,TRANSRET> &groupL) { 00240 ++this->m_statNodesVisited; 00241 00242 if (pIt == partitions.end() || this->checkLeaf(level)) { 00243 PERMLIB_DEBUG(std::cout << "LEAF: " << pi << " with t = " << t << std::endl;) 00244 return this->processLeaf(t, level, backtrackLevel, completed, groupK, groupL); 00245 } 00246 00247 const Partition& sigma = *((*pIt).first); 00248 const RefinementPtr& ref = (*pIt).second; 00249 ++pIt; 00250 00251 unsigned int s = ref->alternatives(); 00252 const bool isBacktrack = ref->type() == Backtrack; 00253 const bool isGroup = ref->type() == Group; 00254 const PERM* tForRefinement = &t; 00255 00256 if (isGroup) { 00257 GroupRefinement<PERM,TRANS>* gref = static_cast<GroupRefinement<PERM,TRANS>*>(ref.get()); 00258 if (this->m_bsgs2 && gref->bsgs() == *this->m_bsgs2) { 00259 tForRefinement = t2; 00260 } 00261 } 00262 00263 ref->sort(*this->m_sorter, &pi); 00264 typedef typename Refinement<PERM>::RefinementPtrIterator RefIt; 00265 for (RefIt rIt = ref->backtrackBegin(); rIt != ref->backtrackEnd(); ++rIt) { 00266 if (isBacktrack && s < groupK.U[backtrackLevel].size()) { 00267 PERMLIB_DEBUG(std::cout << "PRUNE the rest: s=" << s << " < " << groupK.U[backtrackLevel].size() << std::endl;) 00268 this->m_statNodesPrunedCosetMinimality += s; 00269 break; 00270 } 00271 00272 --s; 00273 RefinementPtr ref2 = *rIt; 00274 00275 const unsigned int oldFixPointsSize = pi.fixPointsSize(); 00276 PERMLIB_DEBUG(std::cout << " refinement from " << pi << std::endl;) 00277 const unsigned int strictRefinement = ref2->apply2(pi, *tForRefinement); 00278 PERMLIB_DEBUG(std::cout << " to " << pi << " with " << strictRefinement << std::endl;) 00279 PERMLIB_DEBUG(for(unsigned int jj=0; jj<level; ++jj) std::cout << " ";) 00280 PERMLIB_DEBUG(std::cout << "NODE " << sigma << " ~~~> " << pi << std::endl;) 00281 /* 00282 for (unsigned int q = 0; q < level; ++q) std::cout << " "; 00283 std::cout << " " << level << ": " << sigma << " <-> " << pi2 << " from " << pi << std::endl; 00284 for (unsigned int q = 0; q < level; ++q) std::cout << " "; 00285 std::cout << " t = " << t << std::endl; 00286 */ 00287 if (!strictRefinement) { 00288 PERMLIB_DEBUG(std::cout << "no strict refinement " << sigma << " -- " << pi << std::endl;) 00289 ++this->m_statNodesPrunedChildRestriction; 00290 continue; 00291 } 00292 if (pi.cells() != sigma.cells()) { 00293 PERMLIB_DEBUG(std::cout << "cell number mismatch " << sigma << " -- " << pi << std::endl;) 00294 ref2->undo(pi, strictRefinement); 00295 ++this->m_statNodesPrunedChildRestriction; 00296 continue; 00297 } 00298 if (pi.fixPointsSize() != sigma.fixPointsSize()) { 00299 PERMLIB_DEBUG(std::cout << "fix point number mismatch " << sigma << " -- " << pi << std::endl;) 00300 ref2->undo(pi, strictRefinement); 00301 ++this->m_statNodesPrunedChildRestriction; 00302 continue; 00303 } 00304 PERM tG(t); 00305 PERM* tH = 0; 00306 if (pi.fixPointsSize() != oldFixPointsSize) { 00307 if (!updateMappingPermutation(this->m_bsgs, sigma, pi, tG)) { 00308 PERMLIB_DEBUG(std::cout << "no t found " << sigma << " -- " << pi << "; tG = " << tG << std::endl;) 00309 ref2->undo(pi, strictRefinement); 00310 ++this->m_statNodesPrunedChildRestriction; 00311 continue; 00312 } 00313 if (this->m_bsgs2) { 00314 tH = new PERM(*t2); 00315 if (!updateMappingPermutation(*this->m_bsgs2, sigma, pi, *tH)) { 00316 PERMLIB_DEBUG(std::cout << "no t found " << sigma << " -- " << pi << "; tH = " << tH << std::endl;) 00317 ref2->undo(pi, strictRefinement); 00318 ++this->m_statNodesPrunedChildRestriction; 00319 continue; 00320 } 00321 } 00322 } 00323 if (this->m_pruningLevelDCM && isBacktrack) { 00324 if (this->pruneDCM(tG, backtrackLevel, groupK, groupL)) { 00325 ++this->m_statNodesPrunedCosetMinimality2; 00326 ref2->undo(pi, strictRefinement); 00327 continue; 00328 } 00329 } 00330 unsigned int ret = search(pIt, pi, tG, tH ? tH : t2, level+1, isBacktrack ? (backtrackLevel + 1) : backtrackLevel, completed, groupK, groupL); 00331 delete tH; 00332 PERMLIB_DEBUG(std::cout << "retract " << strictRefinement << " from " << pi << " to ";) 00333 ref2->undo(pi, strictRefinement); 00334 PERMLIB_DEBUG(std::cout << pi << std::endl;) 00335 if (BaseSearch<BSGSIN,TRANSRET>::m_stopAfterFirstElement && ret == 0) 00336 return 0; 00337 if (ret < level) 00338 return ret; 00339 } 00340 00341 completed = std::min(completed, level); 00342 return level; 00343 } 00344 00345 template<class BSGSIN,class TRANSRET> 00346 bool RBase<BSGSIN,TRANSRET>::updateMappingPermutation(const BSGSIN& bsgs, const Partition& sigma, const Partition& pi, PERM& t2) const { 00347 typedef std::vector<unsigned int>::const_iterator FixIt; 00348 std::vector<dom_int>::const_iterator bIt; 00349 unsigned int i = 0; 00350 FixIt fixSigmaIt = sigma.fixPointsBegin(); 00351 const FixIt fixSigmaEndIt = sigma.fixPointsEnd(); 00352 FixIt fixPiIt = pi.fixPointsBegin(); 00353 PERMLIB_DEBUG(print_iterable(bsgs.B.begin(), bsgs.B.end(), 1, "B ");) 00354 PERMLIB_DEBUG(print_iterable(fixSigmaIt, fixSigmaEndIt, 1, "Sigma");) 00355 for (bIt = bsgs.B.begin(); bIt != bsgs.B.end(); ++bIt, ++i) { 00356 PERMLIB_DEBUG(std::cout << " base: " << (*bIt)+1 << std::endl;) 00357 while (fixSigmaIt != fixSigmaEndIt && *fixSigmaIt != *bIt) { 00358 PERMLIB_DEBUG(std::cout << " skipping " << (*fixSigmaIt)+1 << " for " << (*bIt)+1 << std::endl;) 00359 ++fixSigmaIt; 00360 ++fixPiIt; 00361 } 00362 if (fixSigmaIt == fixSigmaEndIt) { 00363 PERMLIB_DEBUG(std::cout << " no more fix point found for " << (*bIt)+1 << std::endl;) 00364 return true; 00365 } 00366 const unsigned int alpha = *fixSigmaIt; 00367 const unsigned int beta = *fixPiIt; 00368 if (t2 / alpha != beta) { 00369 boost::scoped_ptr<PERM> u_beta(bsgs.U[i].at(t2 % beta)); 00370 if (u_beta) { 00371 //std::cout << " multiply with " << *u_beta << " for " << alpha+1 << "," << beta+1 << " // base " << bsgs.B[i] + 1<< std::endl; 00372 t2 ^= *u_beta; 00373 } else { 00374 //std::cout << "could not find a u_b with " << (t2 % beta) << " at " << i << "--" << bsgs.B[i] << " -- " << &bsgs << std::endl; 00375 return false; 00376 } 00377 } 00378 00379 ++fixSigmaIt; 00380 ++fixPiIt; 00381 } 00382 return true; 00383 } 00384 00385 template<class BSGSIN,class TRANSRET> 00386 const std::vector<dom_int>& RBase<BSGSIN,TRANSRET>::subgroupBase() const { 00387 return m_subgroupBase; 00388 } 00389 00390 } 00391 } 00392 00393 #endif // -- RBASE_H_