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