Generated on Tue Dec 16 13:34:01 2008 for ell-3.0.0 by doxygen 1.5.1

src/ell/rna/S_RNAfe_SingleM.cc

Go to the documentation of this file.
00001 #include "ell/rna/S_RNAfe_SingleM.hh"
00002 #include <cmath>    // for getFitness()
00003 #include <biu/assertbiu.hh>
00004 #include "biu/RandomNumberFactory.hh"
00005 
00006 
00007 namespace ell
00008 {
00009 
00010     // the specific ID string of this class
00011     const std::string S_RNAfe_SingleM::ID = std::string("ell::S_RNAfe_SingleM");
00012 
00013 
00014     const size_t 
00015         S_RNAfe_SingleM::RandomNeighborList::ItState::SWITCH_MODE_THRESHOLD = 20;
00016         
00018 // S_RNAfe_SingleM
00020 
00021     S_RNAfe_SingleM::S_RNAfe_SingleM(
00022                 const std::string& rnaSeqStr, 
00023                 const std::string& rnaStructBracketNotStr)
00024      :  rnaFreeEnergy(rnaSeqStr, rnaStructBracketNotStr)
00025     {
00026         assertbiu(rnaFreeEnergy.isValid(),
00027          "S_RNAfe_SingleM has to be a valid state.");
00028          
00029     }
00030 
00031     S_RNAfe_SingleM::S_RNAfe_SingleM(const RNAFreeEnergy& rna)
00032      :  rnaFreeEnergy(rna)
00033     {
00034         assertbiu(rnaFreeEnergy.isValid(),
00035          "S_RNAfe_SingleM has to be a valid state.");
00036          
00037     }
00038     
00039     S_RNAfe_SingleM::S_RNAfe_SingleM(
00040                 const S_RNAfe_SingleM& rnaFESMS)
00041      : rnaFreeEnergy(rnaFESMS.rnaFreeEnergy)
00042     {
00043     }
00044 
00045     S_RNAfe_SingleM::~S_RNAfe_SingleM() {
00046     }
00047 
00048     /* Access to a State subclass specific ID string to identify 
00049      * instances of this class.
00050      * @return the subclass specific ID string
00051      */
00052     const std::string& 
00053     S_RNAfe_SingleM::getID( void ) const {
00054         return S_RNAfe_SingleM::ID;
00055     }
00056 
00057     
00058     S_RNAfe_SingleM::SingleMove 
00059     S_RNAfe_SingleM::firstValidSingleMove() const {
00060         return searchForNextValidSingleMove(0,1);
00061     }
00062 
00063     S_RNAfe_SingleM::SingleMove 
00064     S_RNAfe_SingleM::nextValidSingleMove( 
00065             const SingleMove lastSingleMove) const 
00066     {
00067         // if lastSingleMove.first is an opening bond, this bond was already
00068         // modified in the last move
00069         // closing bonds at lastSingleMove.first are not for interest
00070         if (rnaFreeEnergy.getStructureRef().at(lastSingleMove.first) != biu::RNAStructure::STRUCT_UNBOUND ) {
00071             return searchForNextValidSingleMove(lastSingleMove.first+1,
00072                 lastSingleMove.first+2);
00073         }
00074         else {
00075             return searchForNextValidSingleMove(lastSingleMove.first,
00076                 lastSingleMove.second+1);
00077         }
00078     }
00079 
00080     S_RNAfe_SingleM::SingleMove 
00081     S_RNAfe_SingleM:: searchForNextValidSingleMove( size_t pos1, 
00082                                 size_t pos2) const 
00083     {
00084         assertbiu(pos1 < pos2,
00085          "RNA single move search pos 1 has to be less than pos 2");
00086         
00087 
00088     if (pos2==rnaFreeEnergy.getLength() && pos1 == 0) {
00089         pos1++;
00090         pos2 = pos1+1;
00091     }
00092 //          // no next valid single move found
00093 //          return SingleMove(0,biu::RNAStructure::INVALID_INDEX);
00094             
00095         
00096         // distance between pos1 and pos2 has to be at least
00097         //  MIN_LOOP_LENGTH
00098         while ( (pos1+rnaFreeEnergy.getMinLoopLength())
00099                  < rnaFreeEnergy.getLength()) 
00100         {
00101             if ( rnaFreeEnergy.getStructureRef().at(pos1) == biu::RNAStructure::STRUCT_UNBOUND )
00102             {
00103                 for(pos2 = rnaFreeEnergy.nextCompatibleSingleMovePos(pos2-1);
00104                     pos2 < rnaFreeEnergy.getLength();
00105                     pos2 = rnaFreeEnergy.nextCompatibleSingleMovePos(pos2))
00106                 {
00107                     if ((pos2 - pos1) > rnaFreeEnergy.getMinLoopLength() &&
00108                          rnaFreeEnergy.isAllowedBasePair(pos1, pos2)) {
00109                             // found bases which are able to pair
00110                             // -> exit loop
00111                         return SingleMove(pos1,pos2);
00112                         break;
00113                     }
00114                 }
00115             }
00116             else if ( rnaFreeEnergy.getStructureRef().at(pos1) == biu::RNAStructure::STRUCT_BND_OPEN )
00117             {
00118                 pos2 = rnaFreeEnergy.getClosingBond(pos1);
00119                     // existing bond found -> exit loop
00120                 return SingleMove(pos1,pos2);
00121             }
00122                 // loop wasn't exited since no next valid single move was found
00123             pos1++;
00124             pos2 = pos1;
00125         }
00126         // no next valid single move found
00127         return SingleMove(0,biu::RNAStructure::INVALID_INDEX);
00128     }
00129 
00130     bool    
00131     S_RNAfe_SingleM::operator== (const State& state2) const {
00132         return rnaFreeEnergy == (dynamic_cast<
00133             const S_RNAfe_SingleM&>(state2)).rnaFreeEnergy;
00134     }
00135 
00136     bool
00137     S_RNAfe_SingleM::operator!= (const State& state2) const {
00138         return rnaFreeEnergy != (dynamic_cast<
00139             const S_RNAfe_SingleM&>(state2)).rnaFreeEnergy;
00140     }
00141 
00142     bool
00143     S_RNAfe_SingleM::operator< (const State& rna2) const
00144     {
00145         if (this == &rna2) {
00146             return false;
00147         }
00148         const S_RNAfe_SingleM* s2 = dynamic_cast<const S_RNAfe_SingleM*>(&rna2);
00149         assertbiu(s2!=NULL, "rna2 is no S_RNAfe_SingleM object");
00150         assertbiu(rnaFreeEnergy.getLength() == s2->rnaFreeEnergy.getLength()
00151                     , "the two RNA molecules differ in length");
00152         
00153         return (this->getEnergy() < s2->getEnergy())
00154                 || (this->getEnergy() == s2->getEnergy()
00155                         && std::lexicographical_compare(
00156                                 this->rnaFreeEnergy.getStructureRef().begin()
00157                                 , this->rnaFreeEnergy.getStructureRef().end()
00158                                 , s2->rnaFreeEnergy.getStructureRef().begin()
00159                                 , s2->rnaFreeEnergy.getStructureRef().end()) );
00160     }
00161     
00162     double  
00163     S_RNAfe_SingleM::getEnergy(void) const {
00164         return rnaFreeEnergy.getEnergy();
00165     }
00166 
00167     unsigned int    
00168     S_RNAfe_SingleM::getMinimalDistance(
00169          const State& _state2) const {
00170         S_RNAfe_SingleM state2 = 
00171          dynamic_cast<const S_RNAfe_SingleM&>(_state2);
00172 
00173         const biu::Structure& state1Struct = rnaFreeEnergy.getStructureRef();
00174         const biu::Structure& state2Struct = state2.rnaFreeEnergy.getStructureRef();
00175 
00176         assertbiu(state2Struct.size() == state1Struct.size(),
00177         "Length of both RNA has to be equal to compute the distance.");
00178             
00179         unsigned int distance = 0;
00180         for (size_t i = 0; i < state1Struct.size(); i++){
00181             // bond maintained or bond opened and replaced by another one
00182             if (state1Struct[i] == biu::RNAStructure::STRUCT_BND_OPEN &&
00183                state2Struct[i] == biu::RNAStructure::STRUCT_BND_OPEN ) {
00184                 // bond opened and replaced by another one
00185                 if (rnaFreeEnergy.getClosingBond(i) !=
00186                      state2.rnaFreeEnergy.getClosingBond(i)) {
00187                     distance += 2;
00188                 }
00189             }
00190             // bond opened or closed
00191             else if (state1Struct[i] == biu::RNAStructure::STRUCT_BND_OPEN  ||
00192                       state2Struct[i] == biu::RNAStructure::STRUCT_BND_OPEN ) {
00193                 distance += 1;
00194             }
00195         }
00196         return distance;
00197     }
00198 
00199     S_RNAfe_SingleM* 
00200     S_RNAfe_SingleM::clone(State* toFill) const {
00201         if (toFill == NULL) {
00202             return new S_RNAfe_SingleM(*this);
00203         }
00204           // cast toFill
00205         assertbiu( dynamic_cast<S_RNAfe_SingleM*>(toFill) != NULL
00206                     , "given State is no S_RNAfe_SingleM instance");
00207         S_RNAfe_SingleM* s = static_cast<S_RNAfe_SingleM*>(toFill);
00208         
00209           // copy data
00210         s->rnaFreeEnergy = this->rnaFreeEnergy; 
00211         
00212         return s;
00213     }
00214     
00215     State*
00216     S_RNAfe_SingleM::fromString(const std::string& stringRep) const {
00217           // parse structure
00218         std::string backboneRep = ".()";
00219         size_t pos = stringRep.find_first_not_of(backboneRep);
00220         const std::string rnaStructBracketDotStr = stringRep.substr(0, pos-1);
00221         assertbiu(rnaStructBracketDotStr.size() != 0, "string does not contain a dot-bracket structure encoding");
00222         assertbiu(rnaStructBracketDotStr.size() == this->rnaFreeEnergy.getLength(), "string is too short for this RNA molecule");
00223           // create new copy of this
00224         S_RNAfe_SingleM* ret = this->clone();
00225           // set structure
00226         ret->rnaFreeEnergy.setStructure(
00227                 ret->rnaFreeEnergy.getStructureAlphabet()->getSequence(
00228                         rnaStructBracketDotStr));
00229         return ret;
00230     }
00231     
00232     std::string 
00233     S_RNAfe_SingleM::toString() const{
00234         return rnaFreeEnergy.getStructureString();
00235     }
00236     
00237     std::string& 
00238     S_RNAfe_SingleM::toString( std::string & toFill ) const {
00239           // nothing better so far
00240         toFill = rnaFreeEnergy.getStructureString();
00241         return toFill;
00242     }
00243     
00244 
00245     std::string 
00246     S_RNAfe_SingleM::getSequence() const {
00247         return rnaFreeEnergy.getSequenceString();
00248     }
00249     
00250     std::string 
00251     S_RNAfe_SingleM::getStructure() const {
00252         return rnaFreeEnergy.getStructureString();
00253     }
00254 
00255     State::NeighborListPtr 
00256     S_RNAfe_SingleM::getNeighborList() const {
00257         return NeighborListPtr(new NeighborList(this));
00258     }
00259 
00260     State::NeighborListPtr
00261     S_RNAfe_SingleM::getRandomNeighborList() const 
00262     {
00263         return NeighborListPtr(new RandomNeighborList(this));
00264     }
00265 
00266     // NeighborList
00267     S_RNAfe_SingleM::NeighborList::NeighborList(
00268                 const S_RNAfe_SingleM* _origin) 
00269      : origin(_origin) 
00270     {
00271     }
00272 
00273     S_RNAfe_SingleM::NeighborList::~NeighborList() {
00274     }
00275 
00276     State* 
00277     S_RNAfe_SingleM::NeighborList::first(
00278               State::NeighborList::ItState** itstate) const {
00279         if (*itstate != NULL) delete *itstate;
00280 
00281         ItState* myitstate = new ItState();
00282         *itstate = myitstate;
00283         
00284         S_RNAfe_SingleM* rnaState = NULL;
00285 
00286         SingleMove firstMove = origin->firstValidSingleMove();
00287         
00288 
00289         if (firstMove != SingleMove(0,biu::RNAStructure::INVALID_INDEX)) {
00290             rnaState = new S_RNAfe_SingleM(*origin);
00291             rnaState->rnaFreeEnergy.applySingleMoveInPlace(firstMove);
00292             myitstate->lastMove = firstMove;
00293         }
00294         return rnaState;
00295     }
00296 
00297     State* 
00298     S_RNAfe_SingleM::NeighborList::next(
00299               State::NeighborList::ItState* itstate, State* elem) const {
00300         assertbiu(elem != NULL, "Elem is not allowed to be NULL");
00301 
00302         ItState* myitstate = static_cast<ItState*>(itstate);
00303 
00304         // find next move
00305         SingleMove nextMove = origin->nextValidSingleMove(myitstate->lastMove);
00306         // no more moves possible
00307         if (nextMove == SingleMove(0,biu::RNAStructure::INVALID_INDEX)) {
00308             return NULL;
00309         }
00310         else {
00311             S_RNAfe_SingleM* rnaState =
00312              dynamic_cast<S_RNAfe_SingleM*>(elem);
00313             // undo last move
00314             rnaState->rnaFreeEnergy.applySingleMoveInPlace(myitstate->lastMove);
00315             // perform next move
00316             rnaState->rnaFreeEnergy.applySingleMoveInPlace(nextMove);
00317             myitstate->lastMove = nextMove;
00318             return rnaState;
00319         }
00320     }
00321 
00322     // RandomNeighborList
00323     S_RNAfe_SingleM::RandomNeighborList::RandomNeighborList(
00324             const S_RNAfe_SingleM* _origin)
00325      :  origin(_origin)
00326     {
00327     }
00328 
00329     S_RNAfe_SingleM::RandomNeighborList::~RandomNeighborList() {
00330     }
00331 
00332     void 
00333     S_RNAfe_SingleM::RandomNeighborList::switchMode(
00334             ItState* itstate) const 
00335     {
00336 
00337         ItState* myitstate = static_cast<ItState*>(itstate);
00338 
00339         SingleMove nextMove = origin->firstValidSingleMove();
00340 
00341         // loop while valid moves exist
00342         while (nextMove != SingleMove(0,biu::RNAStructure::INVALID_INDEX)) {
00343             // add all valid and unchosen moves to unchosenValidMoves
00344             if (myitstate->chosenMoves.find(nextMove) == 
00345                  myitstate->chosenMoves.end()) {
00346                 myitstate->unchosenValidMoves.push_back(nextMove);
00347             }
00348             nextMove = origin->nextValidSingleMove(nextMove);
00349         }
00350     }
00351 
00352     State* 
00353     S_RNAfe_SingleM::RandomNeighborList::first(
00354               State::NeighborList::ItState** itstate) const 
00355     {
00356         if (*itstate != NULL) delete *itstate;
00357 
00358         *itstate = new ItState();
00359         S_RNAfe_SingleM* rnaState = 
00360          new S_RNAfe_SingleM(*origin);
00361 
00362         return next(*itstate, rnaState);
00363     }
00364 
00365     State* 
00366     S_RNAfe_SingleM::RandomNeighborList::next(
00367               State::NeighborList::ItState* itstate, State* elem) const 
00368     {
00369         assertbiu(elem != NULL, "Elem is not allowed to be NULL");
00370 
00371         ItState* myitstate = static_cast<ItState*>(itstate);
00372         S_RNAfe_SingleM* rnaState =
00373          dynamic_cast<S_RNAfe_SingleM*>(elem);
00374 
00375         if (myitstate->switchModeValue > 0) {
00376             // positions the move will be applied to
00377             size_t pos1, pos2;
00378             SingleMove nextMove;
00379             bool foundValidMove = false;
00380 
00381             while (!foundValidMove) {
00382                 pos1 = biu::RNF::getRN() % origin->rnaFreeEnergy.getLength();
00383                 pos2 = biu::RNF::getRN() % origin->rnaFreeEnergy.getLength();
00384                 if (pos1 == pos2) {
00385                     continue;
00386                 }
00387                 nextMove = (pos1 < pos2) ? SingleMove(pos1, pos2) : SingleMove(pos2, pos1);
00388         
00389                 if (origin->rnaFreeEnergy.isValidSingleMove(nextMove)) {
00390                     if (!(myitstate->chosenMoves.insert(nextMove).second)) {
00391                         // move has already been chosen --> try next
00392                         myitstate->switchModeValue--;
00393                         if (myitstate->switchModeValue == 0) {
00394                             switchMode(myitstate);
00395                             break;
00396                         }
00397                     }
00398                         // move hasn't been chosen yet and is valid
00399                         // set rnaState to the state of origin of the neighbors
00400                     rnaState->rnaFreeEnergy = origin->rnaFreeEnergy;
00401                     rnaState->rnaFreeEnergy.applySingleMoveInPlace(
00402                      nextMove);
00403                     foundValidMove = true;
00404                 }   
00405             }
00406         } 
00407         if (myitstate->switchModeValue == 0) {
00408             size_t n = myitstate->unchosenValidMoves.size();
00409             if (n == 0) {
00410                 return NULL;
00411             }
00412             else {
00413                 size_t i = biu::RNF::getRN() % n;
00414                     // set rnaState to the state of origin of the neighbors
00415                 rnaState->rnaFreeEnergy = origin->rnaFreeEnergy;
00416                 rnaState->rnaFreeEnergy.applySingleMoveInPlace(
00417                  myitstate->unchosenValidMoves[i]);
00418                 myitstate->unchosenValidMoves[i] =
00419                  myitstate->unchosenValidMoves[n-1];
00420                 myitstate->unchosenValidMoves.resize(n-1);
00421             }
00422         }
00423         return rnaState;
00424     }
00425     
00426     State*  
00427     S_RNAfe_SingleM::getRandomNeighbor(State* inPlaceNeigh) const 
00428     {
00429         S_RNAfe_SingleM* neigh = NULL;
00430             // copy this state
00431         if (inPlaceNeigh != NULL) {
00432             neigh = dynamic_cast<S_RNAfe_SingleM*>(inPlaceNeigh);
00433             assertbiu(neigh != NULL, "inPlaceNeigh is no S_RNAfe_SingleM"); 
00434             (*neigh) = (*this); 
00435         } else {
00436             neigh = new S_RNAfe_SingleM(*this);
00437         }
00438         
00439         size_t pos1, pos2;
00440         SingleMove nextMove;
00441 
00442         while (true) {
00443             pos1 = biu::RNF::getRN() % rnaFreeEnergy.getLength();
00444             pos2 = biu::RNF::getRN() % rnaFreeEnergy.getLength();
00445             if (pos1 == pos2) {
00446                 continue;
00447             }
00448             nextMove = (pos1 < pos2) ? SingleMove(pos1, pos2) : SingleMove(pos2, pos1);
00449             // move has already been chosen
00450             if (rnaFreeEnergy.isValidSingleMove(nextMove)) {
00451                     // set rnaState to the state of origin of the neighbors
00452                 neigh->rnaFreeEnergy.applySingleMoveInPlace( nextMove );
00453                 break;
00454             }   
00455         }
00456         return neigh;
00457     }
00458 
00459     
00460 
00462     
00463         /* Access to a compressed sequence representation of the state.
00464          * @return the compressed sequence representation
00465          */
00466     CSequence  
00467     S_RNAfe_SingleM::compress(void) const {
00468         return rnaFreeEnergy.getStructureAlphabet()->compress(rnaFreeEnergy.getStructureRef());
00469     }
00470     
00471         /* Access to a compressed sequence representation of the state.
00472          * @param toFill a data structure to write the compressed 
00473          *               representation too
00474          * @return the compressed sequence representation
00475          */
00476     CSequence&  
00477     S_RNAfe_SingleM::compress(CSequence& toFill) const {
00478         toFill = rnaFreeEnergy.getStructureAlphabet()->compress(rnaFreeEnergy.getStructureRef());
00479         return toFill;
00480     }
00481     
00482         /* Uncompresses a compressed sequencce representation into a new
00483          * State object.
00484          * @param cseq the compressed sequence representation of a state
00485          * @param toFill a state object to uncompress too or NULL if a new 
00486          *               object has to be created
00487          * @return new State object that is encoded in cseq or NULL in error
00488          *         case.
00489          */
00490     State*  
00491     S_RNAfe_SingleM::uncompress(const CSequence& cseq, State* toFill) const {
00492 
00493         biu::Structure str = rnaFreeEnergy.getStructureAlphabet()->decompress(cseq,rnaFreeEnergy.getStructureRef().size());
00494           // error case
00495         if (str.size() != rnaFreeEnergy.getStructureRef().size())
00496             return NULL;
00497 
00498         S_RNAfe_SingleM* ret = NULL;
00499             // use second function
00500         if (toFill == NULL) {
00501             ret = this->clone();
00502         } else {
00503             ret = dynamic_cast<S_RNAfe_SingleM*>(toFill);
00504             assertbiu(ret!=NULL, "given State toFill is no S_RNAfe_SingleM instance");
00505         }
00506             
00507         ret->rnaFreeEnergy.setStructure(str);
00508         
00509         return ret;
00510     }
00511     
00512         /* Uncompresses a compressed sequencce representation into a this
00513          * State object.
00514          * @param cseq the compressed sequence representation of a state
00515          * @return this or NULL in error case
00516          */
00517     State*  
00518     S_RNAfe_SingleM::uncompress(const CSequence& cseq) {
00519         
00520         biu::Structure str = rnaFreeEnergy.getStructureAlphabet()->decompress(cseq,rnaFreeEnergy.getStructureRef().size());
00521         
00522         if (str.size() != rnaFreeEnergy.getStructureRef().size())
00523             return NULL;
00524         
00525         rnaFreeEnergy.setStructure(str);
00526         
00527         return this;
00528     }
00529     
00530 
00531 } // namespace ell