main page namespaces classes files
ELL home

src/ell/rna/S_RNAfe_SingleM_TB.cc

Go to the documentation of this file.

#include "ell/rna/S_RNAfe_SingleM_TB.hh"
#include <cmath>    // for getFitness()
#include <biu/assertbiu.hh>
#include "biu/RandomNumberFactory.hh"


namespace ell
{

    // the specific ID string of this class
    const std::string S_RNAfe_SingleM_TB::ID = std::string("ell::S_RNAfe_SingleM_TB");


    const size_t 
        S_RNAfe_SingleM_TB::RandomNeighborList::ItState::SWITCH_MODE_THRESHOLD = 20;
        
// S_RNAfe_SingleM_TB

    S_RNAfe_SingleM_TB::S_RNAfe_SingleM_TB(
                const std::string& rnaSeqStr, 
                const std::string& rnaStructBracketNotStr)
     :  rnaFreeEnergy(rnaSeqStr, rnaStructBracketNotStr)
    {
        assertbiu(rnaFreeEnergy.isValid(),
         "S_RNAfe_SingleM_TB has to be a valid state.");
         
    }

    S_RNAfe_SingleM_TB::S_RNAfe_SingleM_TB(const RNAFreeEnergy_TB& rna)
     :  rnaFreeEnergy(rna)
    {
        assertbiu(rnaFreeEnergy.isValid(),
         "S_RNAfe_SingleM_TB has to be a valid state.");
         
    }
    
    S_RNAfe_SingleM_TB::S_RNAfe_SingleM_TB(
                const S_RNAfe_SingleM_TB& rnaFESMS)
     : rnaFreeEnergy(rnaFESMS.rnaFreeEnergy)
    {
    }

    S_RNAfe_SingleM_TB::~S_RNAfe_SingleM_TB() 
    {
    }

    /* Access to a State subclass specific ID string to identify 
     * instances of this class.
     * @return the subclass specific ID string
     */
    const std::string& 
    S_RNAfe_SingleM_TB::getID( void ) const {
        return S_RNAfe_SingleM_TB::ID;
    }

    
    S_RNAfe_SingleM_TB::SingleMove 
    S_RNAfe_SingleM_TB::firstValidSingleMove() const {
        SingleMove m(   biu::RNAStructure_TB::INVALID_INDEX
                        , biu::RNAStructure_TB::INVALID_INDEX );
        if (rnaFreeEnergy.getNextSingleMove(m.first, m.second)) {
            return m;
        } else {
            return SingleMove(0,biu::RNAStructure_TB::INVALID_INDEX);
        }
    }

    S_RNAfe_SingleM_TB::SingleMove 
    S_RNAfe_SingleM_TB::nextValidSingleMove( 
            const SingleMove lastSingleMove) const 
    {
        SingleMove ret(lastSingleMove);
        if (rnaFreeEnergy.getNextSingleMove( ret.first, ret.second ))
        {
            return ret;
        } else {
            return SingleMove(0,biu::RNAStructure_TB::INVALID_INDEX);
        }
//      // if lastSingleMove.first is an opening bond, this bond was already
//      // modified in the last move
//      // closing bonds at lastSingleMove.first are not for interest
//      if (rnaFreeEnergy.getCorrespondingBase(lastSingleMove.first) != biu::RNAStructure_TB::INVALID_INDEX ) {
//          return searchForNextValidSingleMove(lastSingleMove.first+1,
//              lastSingleMove.first+2);
//      }
//      else {
//          return searchForNextValidSingleMove(lastSingleMove.first,
//              lastSingleMove.second+1);
//      }
    }

//  S_RNAfe_SingleM_TB::SingleMove 
//  S_RNAfe_SingleM_TB:: searchForNextValidSingleMove(  size_t pos1, 
//                              size_t pos2) const 
//  {
//      assertbiu(pos1 < pos2,
//       "RNA single move search pos 1 has to be less than pos 2");
//      
//
//  if (pos2>=rnaFreeEnergy.getLength() ) {
//      pos1++;
//      pos2 = pos1+1;
//  }
//
//      
//      // distance between pos1 and pos2 has to be at least
//      //  MIN_LOOP_LENGTH
//      while ( (pos1+rnaFreeEnergy.getMinLoopLength())
//               < rnaFreeEnergy.getLength()) 
//      {
//          if ( rnaFreeEnergy.getCorrespondingBase(pos1) == biu::RNAStructure_TB::INVALID_INDEX )
//          {
//              for(pos2 = rnaFreeEnergy.nextCompatibleSingleMovePos(pos2-1);
//                  pos2 < rnaFreeEnergy.getLength();
//                  pos2 = rnaFreeEnergy.nextCompatibleSingleMovePos(pos2))
//              {
//                  if ((pos2 - pos1) > rnaFreeEnergy.getMinLoopLength() &&
//                       rnaFreeEnergy.isAllowedBasePair(pos1, pos2)) {
//                          // found bases which are able to pair
//                          // -> exit loop
//                      return SingleMove(pos1,pos2);
//                      break;
//                  }
//              }
//          }
//          else if ( rnaFreeEnergy.getCorrespondingBase(pos1) > pos1 )
//          {
//                  // existing bond found -> exit loop
//              return SingleMove(pos1,rnaFreeEnergy.getCorrespondingBase(pos1));
//          }
//              // loop wasn't exited since no next valid single move was found
//          pos1++;
//          pos2 = pos1 + 1;
//      }
//      // no next valid single move found
//      return SingleMove(0,biu::RNAStructure_TB::INVALID_INDEX);
//  }

    bool    
    S_RNAfe_SingleM_TB::operator== (const State& state2) const {
          // check if convertable
        const S_RNAfe_SingleM_TB * tmp 
            = dynamic_cast< const S_RNAfe_SingleM_TB*>(&state2);
        return (tmp==NULL? false : rnaFreeEnergy == tmp->rnaFreeEnergy);
    }

    bool
    S_RNAfe_SingleM_TB::operator!= (const State& state2) const {
        return !(this->operator==(state2));
    }

    bool
    S_RNAfe_SingleM_TB::operator< (const State& rna2) const
    {
        if (this == &rna2) {
            return false;
        }
        assertbiu(dynamic_cast<const S_RNAfe_SingleM_TB*>(&rna2)!=NULL
                    , "rna2 is no S_RNAfe_SingleM_TB object");
        const S_RNAfe_SingleM_TB* s2 = static_cast<const S_RNAfe_SingleM_TB*>(&rna2);
        assertbiu(rnaFreeEnergy.getLength() == s2->rnaFreeEnergy.getLength()
                    , "the two RNA molecules differ in length");
        
        return (this->getEnergy() < s2->getEnergy())
                || 
                (this->getEnergy() == s2->getEnergy()
                && rnaFreeEnergy.getStructureStringRef().compare(
                        s2->rnaFreeEnergy.getStructureStringRef() ) < 0 ) ;
    }
    
    double  
    S_RNAfe_SingleM_TB::getEnergy(void) const {
        return rnaFreeEnergy.getEnergy();
    }

    unsigned int    
    S_RNAfe_SingleM_TB::getMinimalDistance( const State& _state2) const 
    {
        assertbiu(dynamic_cast<const S_RNAfe_SingleM_TB*>(&_state2) != NULL
                , "state to calculate distance to is no S_RNAfe_SingleM_TB");
        const S_RNAfe_SingleM_TB* state2 = 
         static_cast<const S_RNAfe_SingleM_TB*>(&_state2);

        const std::string& state1Struct = this->rnaFreeEnergy.getStructureStringRef();
        const std::string& state2Struct = state2->rnaFreeEnergy.getStructureStringRef();

        assertbiu(state2Struct.size() == state1Struct.size(),
        "Length of both RNA has to be equal to compute the distance.");
            
        unsigned int distance = 0;
        
        const char BOND_OPEN = this->rnaFreeEnergy.getStructureAlphabet()->getString(biu::RNAStructure_TB::STRUCT_BND_OPEN).at(0); 
        
        for (size_t i = 0; i < state1Struct.size(); i++){
            // bond maintained or bond opened and replaced by another one
            if (    BOND_OPEN == state1Struct[i] 
                    && BOND_OPEN == state2Struct[i] ) 
            {
                // check if bond opened and replaced by another one
                if (    this->rnaFreeEnergy.getCorrespondingBase(i) 
                        != state2->rnaFreeEnergy.getCorrespondingBase(i)) 
                {
                    distance += 2;
                }
            }
            // one bond opened or closed
            else if (   BOND_OPEN == state1Struct[i]  
                        || BOND_OPEN == state2Struct[i] ) 
            {
                distance += 1;
            }
        }
        return distance;
    }

    S_RNAfe_SingleM_TB* 
    S_RNAfe_SingleM_TB::clone(State* toFill) const {
        if (toFill == NULL) {
            return new S_RNAfe_SingleM_TB(*this);
        }
          // cast toFill
        assertbiu( dynamic_cast<S_RNAfe_SingleM_TB*>(toFill) != NULL
                    , "given State is no S_RNAfe_SingleM_TB instance");
        S_RNAfe_SingleM_TB* s = static_cast<S_RNAfe_SingleM_TB*>(toFill);
        
          // copy data
        s->rnaFreeEnergy = this->rnaFreeEnergy; 
        
        return s;
    }
    
    
    State*
    S_RNAfe_SingleM_TB::fromString(const std::string& stringRep) const {
          // parse structure
        std::string backboneRep = ".()";
        size_t pos = stringRep.find_first_not_of(backboneRep);
        const std::string rnaStructBracketDotStr = stringRep.substr(0, pos-1);
        assertbiu(rnaStructBracketDotStr.size() != 0, "string does not contain a dot-bracket structure encoding");
        assertbiu(rnaStructBracketDotStr.size() == this->rnaFreeEnergy.getLength(), "string is too short for this RNA molecule");
          // create new copy of this
        S_RNAfe_SingleM_TB* ret = this->clone();
          // set structure
        ret->rnaFreeEnergy.setStructure(
                ret->rnaFreeEnergy.getStructureAlphabet()->getSequence(
                        rnaStructBracketDotStr));
        return ret;
    }
    
    State*
    S_RNAfe_SingleM_TB::fromString(const std::string& stringRep, const double energy) const {
          // create state from string
        State* ret = this->fromString(stringRep);
          // set energy
        static_cast<S_RNAfe_SingleM_TB*>(ret)->rnaFreeEnergy.setEnergy(energy);
          // return updated state
        return ret;
    }
    
    std::string 
    S_RNAfe_SingleM_TB::toString() const{
        return rnaFreeEnergy.getStructureString();
    }
    
    std::string& 
    S_RNAfe_SingleM_TB::toString( std::string & toFill ) const {
          // nothing better so far
        toFill = rnaFreeEnergy.getStructureStringRef();
        return toFill;
    }
    
    std::string 
    S_RNAfe_SingleM_TB::getSequence() const {
        return rnaFreeEnergy.getSequenceString();
    }
    
    std::string 
    S_RNAfe_SingleM_TB::getStructure() const {
        return rnaFreeEnergy.getStructureString();
    }

    State::NeighborListPtr 
    S_RNAfe_SingleM_TB::getNeighborList() const {
        return NeighborListPtr(new NeighborList(this));
    }

    State::NeighborListPtr
    S_RNAfe_SingleM_TB::getRandomNeighborList() const 
    {
        return NeighborListPtr(new RandomNeighborList(this));
    }

    // NeighborList
    S_RNAfe_SingleM_TB::NeighborList::NeighborList(
                const S_RNAfe_SingleM_TB* _origin) 
     : origin(_origin) 
    {
    }

    S_RNAfe_SingleM_TB::NeighborList::~NeighborList() {
    }

    State* 
    S_RNAfe_SingleM_TB::NeighborList::first(
              State::NeighborList::ItState** itstate) const {
        if (*itstate != NULL) delete *itstate;

        ItState* myitstate = new ItState();
        *itstate = myitstate;
        
        S_RNAfe_SingleM_TB* rnaState = NULL;

        SingleMove firstMove = origin->firstValidSingleMove();
        

        if (firstMove != SingleMove(0,biu::RNAStructure_TB::INVALID_INDEX)) {
            rnaState = new S_RNAfe_SingleM_TB(*origin);
            rnaState->rnaFreeEnergy.applySingleMoveInPlace(firstMove);
            myitstate->lastMove = firstMove;
        }
        return rnaState;
    }

    State* 
    S_RNAfe_SingleM_TB::NeighborList::next(
              State::NeighborList::ItState* itstate, State* elem) const {
        assertbiu(elem != NULL, "elem is not allowed to be NULL");
        assertbiu(dynamic_cast<S_RNAfe_SingleM_TB*>(elem) != NULL
                , "elem is no instance of S_RNAfe_SingleM_TB");

        ItState* myitstate = static_cast<ItState*>(itstate);

        // find next move
        SingleMove nextMove = origin->nextValidSingleMove(myitstate->lastMove);
        // no more moves possible
        if (nextMove == SingleMove(0,biu::RNAStructure_TB::INVALID_INDEX)) {
            return NULL;
        }
        else {
            S_RNAfe_SingleM_TB* rnaState 
                    = static_cast<S_RNAfe_SingleM_TB*>(elem);
            // undo last move
            rnaState->rnaFreeEnergy.applySingleMoveInPlace(myitstate->lastMove);
            // perform next move
            rnaState->rnaFreeEnergy.applySingleMoveInPlace(nextMove);
            myitstate->lastMove = nextMove;
            return rnaState;
        }
    }

    // RandomNeighborList
    S_RNAfe_SingleM_TB::RandomNeighborList::RandomNeighborList(
            const S_RNAfe_SingleM_TB* _origin)
     :  origin(_origin)
    {
    }

    S_RNAfe_SingleM_TB::RandomNeighborList::~RandomNeighborList() {
    }

    void 
    S_RNAfe_SingleM_TB::RandomNeighborList::switchMode(
            ItState* itstate) const 
    {

        ItState* myitstate = static_cast<ItState*>(itstate);

        SingleMove nextMove = origin->firstValidSingleMove();

        // loop while valid moves exist
        while (nextMove != SingleMove(0,biu::RNAStructure_TB::INVALID_INDEX)) {
            // add all valid and unchosen moves to unchosenValidMoves
            if (myitstate->chosenMoves.find(nextMove) == 
                 myitstate->chosenMoves.end()) {
                myitstate->unchosenValidMoves.push_back(nextMove);
            }
            nextMove = origin->nextValidSingleMove(nextMove);
        }
    }

    State* 
    S_RNAfe_SingleM_TB::RandomNeighborList::first(
              State::NeighborList::ItState** itstate) const 
    {
        if (*itstate != NULL) delete *itstate;

        *itstate = new ItState();
        S_RNAfe_SingleM_TB* rnaState = 
         new S_RNAfe_SingleM_TB(*origin);

        return next(*itstate, rnaState);
    }

    State* 
    S_RNAfe_SingleM_TB::RandomNeighborList::next(
              State::NeighborList::ItState* itstate, State* elem) const 
    {
        assertbiu(elem != NULL, "Elem is not allowed to be NULL");
        assertbiu(dynamic_cast<S_RNAfe_SingleM_TB*>(elem) != NULL
                , "elem is no instance of S_RNAfe_SingleM_TB");

        ItState* myitstate = static_cast<ItState*>(itstate);
        S_RNAfe_SingleM_TB* rnaState = static_cast<S_RNAfe_SingleM_TB*>(elem);

        if (myitstate->switchModeValue > 0) {
            // positions the move will be applied to
            size_t pos1, pos2;
            SingleMove nextMove;
            bool foundValidMove = false;

            while (!foundValidMove) {
                pos1 = biu::RNF::getRN() % origin->rnaFreeEnergy.getLength();
                pos2 = biu::RNF::getRN() % origin->rnaFreeEnergy.getLength();
                if (pos1 == pos2) {
                    continue;
                }
                nextMove = (pos1 < pos2) ? SingleMove(pos1, pos2) : SingleMove(pos2, pos1);
        
                if (origin->rnaFreeEnergy.isValidSingleMove(nextMove)) {
                    if (!(myitstate->chosenMoves.insert(nextMove).second)) {
                        // move has already been chosen --> try next
                        myitstate->switchModeValue--;
                        if (myitstate->switchModeValue == 0) {
                            switchMode(myitstate);
                            break;
                        }
                    }
                        // move hasn't been chosen yet and is valid
                        // set rnaState to the state of origin of the neighbors
                    rnaState->rnaFreeEnergy = origin->rnaFreeEnergy;
                    rnaState->rnaFreeEnergy.applySingleMoveInPlace(
                     nextMove);
                    foundValidMove = true;
                }   
            }
        } 
        if (myitstate->switchModeValue == 0) {
            size_t n = myitstate->unchosenValidMoves.size();
            if (n == 0) {
                return NULL;
            }
            else {
                size_t i = biu::RNF::getRN() % n;
                    // set rnaState to the state of origin of the neighbors
                rnaState->rnaFreeEnergy = origin->rnaFreeEnergy;
                rnaState->rnaFreeEnergy.applySingleMoveInPlace(
                 myitstate->unchosenValidMoves[i]);
                myitstate->unchosenValidMoves[i] =
                 myitstate->unchosenValidMoves[n-1];
                myitstate->unchosenValidMoves.resize(n-1);
            }
        }
        return rnaState;
    }
    
    State*  
    S_RNAfe_SingleM_TB::getRandomNeighbor(State* inPlaceNeigh) const 
    {
        assertbiu(dynamic_cast<S_RNAfe_SingleM_TB*>(inPlaceNeigh) != NULL
                    , "inPlaceNeigh is no S_RNAfe_SingleM_TB"); 
        S_RNAfe_SingleM_TB* neigh = NULL;
            // copy this state
        if (inPlaceNeigh != NULL) {
            neigh = static_cast<S_RNAfe_SingleM_TB*>(inPlaceNeigh);
            (*neigh) = (*this); 
        } else {
            neigh = new S_RNAfe_SingleM_TB(*this);
        }
        
        size_t pos1, pos2;
        SingleMove nextMove;

        while (true) {
            pos1 = biu::RNF::getRN() % rnaFreeEnergy.getLength();
            pos2 = biu::RNF::getRN() % rnaFreeEnergy.getLength();
            if (pos1 == pos2) {
                continue;
            }
            nextMove = (pos1 < pos2) ? SingleMove(pos1, pos2) : SingleMove(pos2, pos1);
            // move has already been chosen
            if (rnaFreeEnergy.isValidSingleMove(nextMove)) {
                    // set rnaState to the state of origin of the neighbors
                neigh->rnaFreeEnergy.applySingleMoveInPlace( nextMove );
                break;
            }   
        }
        return neigh;
    }

    

    
        /* Access to a compressed sequence representation of the state.
         * @return the compressed sequence representation
         */
    CSequence  
    S_RNAfe_SingleM_TB::compress(void) const {
        return rnaFreeEnergy.getStructureAlphabet()->compressS(rnaFreeEnergy.getStructureStringRef());
    }
    
        /* Access to a compressed sequence representation of the state.
         * @param toFill a data structure to write the compressed 
         *               representation too
         * @return the compressed sequence representation
         */
    CSequence&  
    S_RNAfe_SingleM_TB::compress(CSequence& toFill) const {
        toFill = rnaFreeEnergy.getStructureAlphabet()->compressS(rnaFreeEnergy.getStructureStringRef());
        return toFill;
    }
    
        /* Uncompresses a compressed sequencce representation into a new
         * State object.
         * @param cseq the compressed sequence representation of a state
         * @param toFill a state object to uncompress too or NULL if a new 
         *               object has to be created
         * @return new State object that is encoded in cseq or NULL in error
         *         case.
         */
    State*  
    S_RNAfe_SingleM_TB::uncompress(const CSequence& cseq, State* toFill) const {

        assertbiu(toFill == NULL || dynamic_cast<S_RNAfe_SingleM_TB*>(toFill) != NULL
                    , "toFill is no S_RNAfe_SingleM_TB"); 

        biu::Structure str = rnaFreeEnergy.getStructureAlphabet()->decompress(cseq,rnaFreeEnergy.getStructureStringRef().size());
          // error case
        if (str.size() != rnaFreeEnergy.getStructureStringRef().size())
            return NULL;

        S_RNAfe_SingleM_TB* ret = NULL;
            // use second function
        if (toFill == NULL) {
            ret = this->clone();
        } else {
            ret = static_cast<S_RNAfe_SingleM_TB*>(toFill);
            assertbiu(ret!=NULL, "given State toFill is no S_RNAfe_SingleM_TB instance");
        }
            
        ret->rnaFreeEnergy.setStructure(str);
        
        return ret;
    }
    
        /* Uncompresses a compressed sequencce representation into a this
         * State object.
         * @param cseq the compressed sequence representation of a state
         * @return this or NULL in error case
         */
    State*  
    S_RNAfe_SingleM_TB::uncompress(const CSequence& cseq) {
        
        biu::Structure str = rnaFreeEnergy.getStructureAlphabet()->decompress(cseq,rnaFreeEnergy.getStructureStringRef().size());
        
        if (str.size() != rnaFreeEnergy.getStructureStringRef().size())
            return NULL;
        
        rnaFreeEnergy.setStructure(str);
        
        return this;
    }
    

} // namespace ell