permlib  0.2.6
Library for permutation computations
 All Classes Functions Variables Typedefs Enumerations Friends
include/permlib/search/partition/group_refinement.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 GROUPREFINEMENT_H_
00034 #define GROUPREFINEMENT_H_
00035 
00036 #include <permlib/predicate/pointwise_stabilizer_predicate.h>
00037 #include <permlib/search/partition/refinement.h>
00038 
00039 namespace permlib {
00040 namespace partition {
00041 
00043 template<class PERM,class TRANS>
00044 class GroupRefinement : public Refinement<PERM> {
00045 public:
00047         explicit GroupRefinement(const BSGSCore<PERM,TRANS>& bsgs);
00048         
00049         virtual unsigned int apply(Partition& pi) const;
00050         virtual unsigned int apply2(Partition& pi, const PERM& t) const;
00051         
00052         virtual bool init(Partition& pi);
00053         
00055         const BSGSCore<PERM,TRANS>& bsgs() const;
00056 private:
00057         const BSGSCore<PERM,TRANS>& m_bsgs;
00058         
00059         std::vector<unsigned int> thetaOrbit;
00060         std::vector<int> thetaBorder;
00061         //WARNING: not thread-safe
00063         mutable Partition::vector_t m_myTheta;
00064         
00065         unsigned int apply2(Partition& pi, const PERM* t) const;
00066 };
00067 
00068 template<class PERM,class TRANS>
00069 GroupRefinement<PERM,TRANS>::GroupRefinement(const BSGSCore<PERM,TRANS>& bsgs_) 
00070         : Refinement<PERM>(bsgs_.n, Group), m_bsgs(bsgs_), thetaOrbit(m_bsgs.n), thetaBorder(m_bsgs.n, -1), m_myTheta(m_bsgs.n)
00071 {
00072 }
00073 
00074 template<class PERM,class TRANS>
00075 unsigned int GroupRefinement<PERM,TRANS>::apply(Partition& pi) const {
00076         return apply2(pi, 0);
00077 }
00078 
00079 template<class PERM,class TRANS>
00080 unsigned int GroupRefinement<PERM,TRANS>::apply2(Partition& pi, const PERM& t) const {
00081         return apply2(pi, &t);
00082 }
00083 
00084 template<class PERM,class TRANS>
00085 unsigned int GroupRefinement<PERM,TRANS>::apply2(Partition& pi, const PERM* t) const {
00086         BOOST_ASSERT( this->initialized() );
00087         
00088         Partition::vector_t& myTheta = m_myTheta;
00089         
00090         Partition::vector_t::iterator thetaIt;
00091         Partition::vector_t::iterator thetaBeginIt, thetaEndIt;
00092         Partition::vector_t::const_iterator thetaOrbitIt;
00093         std::list<int>::const_iterator cellPairIt = Refinement<PERM>::m_cellPairs.begin();
00094         unsigned int ret = false;
00095         while (cellPairIt != Refinement<PERM>::m_cellPairs.end()) {
00096                 const int thetaC = *cellPairIt;
00097                 ++cellPairIt;
00098                 if (*cellPairIt < 0) {
00099                         ++cellPairIt;
00100                         continue;
00101                 }
00102                 
00103                 int borderLo = 0;
00104                 if (thetaC > 0)
00105                         borderLo = thetaBorder[thetaC-1];
00106                 thetaBeginIt = myTheta.begin() + borderLo;
00107                 thetaEndIt   = myTheta.begin() + thetaBorder[thetaC];
00108                 
00109                 if (t) {
00110                         for (thetaIt = thetaBeginIt, thetaOrbitIt =  thetaOrbit.begin() + borderLo;
00111                                          thetaIt != thetaEndIt,  thetaOrbitIt != thetaOrbit.begin() + thetaBorder[thetaC];
00112                                          ++thetaIt, ++thetaOrbitIt) 
00113                         {
00114                                 *thetaIt = *t / *thetaOrbitIt;
00115                         }
00116                         std::sort(thetaBeginIt, thetaEndIt);
00117                 }
00118                 
00119                 for (int c = *cellPairIt; c >= 0; c = *(++cellPairIt)) {
00120                         if (t) {
00121                                 PERMLIB_DEBUG(std::cout << "apply theta " << thetaC << "," << c << " t = " << *t << " to " << pi << std::endl;)
00122                         } else {
00123                                 PERMLIB_DEBUG(std::cout << "apply theta " << thetaC << "," << c << " t = 0 to " << pi << std::endl;)
00124                         }
00125                         //print_iterable(thetaBeginIt, thetaEndIt, 1, "theta apply");
00126                         if (pi.intersect(thetaBeginIt, thetaEndIt, c))
00127                                 ++ret;
00128                 }
00129                 ++cellPairIt;
00130         }
00131         
00132         return ret;
00133 }
00134 
00135 template<class PERM,class TRANS>
00136 bool GroupRefinement<PERM,TRANS>::init(Partition& pi) {
00137         unsigned int fixSize = pi.fixPointsSize();
00138         if (fixSize > 0) {
00139                 boost::dynamic_bitset<> orbitCharacterstic(m_bsgs.n);
00140                 
00141                 std::vector<dom_int>::const_iterator Bit;
00142                 Partition::vector_t::const_iterator fixIt = pi.fixPointsBegin();
00143                 Partition::vector_t::const_iterator fixEndIt = pi.fixPointsEnd();
00144                 for (Bit = m_bsgs.B.begin(); Bit != m_bsgs.B.end(); ++Bit) {
00145                         while (fixIt != fixEndIt && *fixIt != *Bit) {
00146                                 PERMLIB_DEBUG(std::cout << "skip " << (*fixIt + 1) << " for " << (*Bit + 1) << std::endl;)
00147                                 ++fixIt;
00148                         }
00149                         if (fixIt == fixEndIt)
00150                                 break;
00151                         ++fixIt;
00152                 }
00153                 //PointwiseStabilizerPredicate<PERM> fixStab(m_bsgs.B.begin(), m_bsgs.B.begin() + std::min(fixSize, static_cast<unsigned int>(m_bsgs.B.size())));
00154 #ifdef PERMLIB_DEBUG_OUTPUT
00155                 print_iterable(m_bsgs.B.begin(), m_bsgs.B.end(), 1, " BSGS ");
00156                 print_iterable(m_bsgs.B.begin(), Bit, 1, "to fix");
00157                 print_iterable(pi.fixPointsBegin(), pi.fixPointsEnd(), 1, "   fix");
00158 #endif
00159                 PointwiseStabilizerPredicate<PERM> fixStab(m_bsgs.B.begin(), Bit);
00160                 std::list<PERM> S_fix;
00161                 BOOST_FOREACH(const typename PERM::ptr& p, m_bsgs.S) {
00162                         if (fixStab(p)) {
00163                                 //std::cout << "$ " << *p << " fixes " << std::min(fixSize, static_cast<unsigned long>(m_bsgs.B.size())) << " points" << std::endl;
00164                                 S_fix.push_back(*p);
00165                         }
00166                 }
00167                 
00168                 unsigned int thetaIndex = 0;
00169                 std::vector<int>::iterator thetaBorderIt = thetaBorder.begin();
00170                 std::vector<unsigned int>::iterator thetaIt = thetaOrbit.begin();
00171                 for (unsigned long alpha = 0; alpha < m_bsgs.n; ++alpha) {
00172                         if (orbitCharacterstic[alpha])
00173                                 continue;
00174                         orbitCharacterstic.flip(alpha);
00175                         std::vector<unsigned int>::iterator thetaOrbitBeginIt = thetaIt;
00176                         *thetaIt = alpha;
00177                         ++thetaIt;
00178                         ++thetaIndex;
00179                         std::vector<unsigned int>::iterator thetaOrbitEndIt = thetaIt;
00180                         
00181                         std::vector<unsigned int>::iterator it;
00182                         for (it = thetaOrbitBeginIt; it != thetaOrbitEndIt; ++it) {
00183                                 unsigned int beta = *it;
00184                                 BOOST_FOREACH(const PERM &p, S_fix) {
00185                                         unsigned int beta_p = p / beta;
00186                                         if (!orbitCharacterstic[beta_p]) {
00187                                                 *thetaIt = beta_p;
00188                                                 ++thetaIt;
00189                                                 ++thetaOrbitEndIt;
00190                                                 ++thetaIndex;
00191                                                 orbitCharacterstic.flip(beta_p);
00192                                         }
00193                                 }
00194                         }
00195                         *thetaBorderIt = thetaIndex;
00196                         ++thetaBorderIt;
00197                 }
00198                 
00199                 thetaIt = thetaOrbit.begin();
00200                 std::vector<unsigned int>::iterator thetaItEnd;
00201                 thetaBorderIt = thetaBorder.begin();
00202                 unsigned int thetaC = 0;
00203                 int oldBorder = 0;
00204                 while (thetaBorderIt != thetaBorder.end() && *thetaBorderIt >= 0) {
00205                         thetaItEnd = thetaOrbit.begin() + *thetaBorderIt;
00206                         std::sort(thetaIt, thetaItEnd);
00207                         
00208                         if (*thetaBorderIt - oldBorder != 1 || std::find(pi.fixPointsBegin(), pi.fixPointsEnd(), *thetaIt) == pi.fixPointsEnd()) {
00209                                 bool hasTheta = false;
00210                                 const unsigned int oldCellNumber = pi.cells();
00211                                 for (unsigned int c = 0; c < oldCellNumber; ++c) {
00212                                         if (pi.cellSize(c) == 1)
00213                                                 continue;
00214                                         
00215                                         //std::cout << "  theta pi = " << pi << std::endl;
00216                                         //print_iterable(thetaIt, thetaItEnd, 1, "theta prep");
00217                                         if (pi.intersect(thetaIt, thetaItEnd, c)) {
00218                                                 PERMLIB_DEBUG(std::cout << "prepare theta " << thetaC << "," << c << std::endl;)
00219                                                 //print_iterable(thetaIt, thetaItEnd, 1, "theta prep");
00220                                                 if (!hasTheta) {
00221                                                         Refinement<PERM>::m_cellPairs.push_back(thetaC);
00222                                                         hasTheta = true;
00223                                                 }
00224                                                 Refinement<PERM>::m_cellPairs.push_back(c);
00225                                                 //std::cout << (thetaIt - thetaOrbit.begin()) << " - " << (thetaItEnd - thetaOrbit.begin()) << std::endl;
00226                                         }
00227                                         //std::cout << "pii = " << pi << std::endl;
00228                                 }
00229                                 
00230                                 if (hasTheta)
00231                                         Refinement<PERM>::m_cellPairs.push_back(-1);
00232                         }
00233                         
00234                         oldBorder = *thetaBorderIt;
00235                         thetaIt = thetaItEnd;
00236                         ++thetaC;
00237                         ++thetaBorderIt;
00238                 }
00239                 //print_iterable(Refinement<PERM>::m_cellPairs.begin(), Refinement<PERM>::m_cellPairs.end(), 0, "cell pairs");
00240                 if (!Refinement<PERM>::m_cellPairs.empty()) {
00241                         typename Refinement<PERM>::RefinementPtr ref(new GroupRefinement<PERM,TRANS>(*this));
00242                         Refinement<PERM>::m_backtrackRefinements.push_back(ref);
00243                         return true;
00244                 }
00245         } 
00246         
00247         return false;
00248 }
00249 
00250 template<class PERM,class TRANS>
00251 const BSGSCore<PERM,TRANS>& GroupRefinement<PERM,TRANS>::bsgs() const {
00252         return m_bsgs;
00253 }
00254 
00255 }
00256 }
00257 
00258 #endif // -- GROUPREFINEMENT_H_