main page namespaces classes files
ELL home

src/ell/rna/S_RNAfe_SingleM.cc

Go to the documentation of this file.

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


namespace ell
{

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

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

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

    S_RNAfe_SingleM::S_RNAfe_SingleM(
                const S_RNAfe_SingleM& rnaFESMS)
     : rnaFreeEnergy(rnaFESMS.rnaFreeEnergy)
    {
    }

    S_RNAfe_SingleM::~S_RNAfe_SingleM() {
    }

    S_RNAfe_SingleM::SingleMove 
    S_RNAfe_SingleM::firstValidSingleMove() const {
        return searchForNextValidSingleMove(0,1);
    }

    S_RNAfe_SingleM::SingleMove 
    S_RNAfe_SingleM::nextValidSingleMove( 
            const SingleMove lastSingleMove) const 
    {
        // 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.getStructureRef().at(lastSingleMove.first) != biu::RNAStructure::STRUCT_UNBOUND ) {
            return searchForNextValidSingleMove(lastSingleMove.first+1,
                lastSingleMove.first+2);
        }
        else {
            return searchForNextValidSingleMove(lastSingleMove.first,
                lastSingleMove.second+1);
        }
    }

    S_RNAfe_SingleM::SingleMove 
    S_RNAfe_SingleM:: 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 == 0) {
        pos1++;
        pos2 = pos1+1;
    }
//          // no next valid single move found
//          return SingleMove(0,biu::RNAStructure::INVALID_INDEX);
            
        
        // distance between pos1 and pos2 has to be at least
        //  MIN_LOOP_LENGTH
        while ( (pos1+rnaFreeEnergy.getMinLoopLength())
                 < rnaFreeEnergy.getLength()) 
        {
            if ( rnaFreeEnergy.getStructureRef().at(pos1) == biu::RNAStructure::STRUCT_UNBOUND )
            {
                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.getStructureRef().at(pos1) == biu::RNAStructure::STRUCT_BND_OPEN )
            {
                pos2 = rnaFreeEnergy.getClosingBond(pos1);
                    // existing bond found -> exit loop
                return SingleMove(pos1,pos2);
            }
                // loop wasn't exited since no next valid single move was found
            pos1++;
            pos2 = pos1;
        }
        // no next valid single move found
        return SingleMove(0,biu::RNAStructure::INVALID_INDEX);
    }

    bool    
    S_RNAfe_SingleM::operator== (const State& state2) const {
        return rnaFreeEnergy == (dynamic_cast<
            const S_RNAfe_SingleM&>(state2)).rnaFreeEnergy;
    }

    bool
    S_RNAfe_SingleM::operator!= (const State& state2) const {
        return rnaFreeEnergy != (dynamic_cast<
            const S_RNAfe_SingleM&>(state2)).rnaFreeEnergy;
    }

    int 
    S_RNAfe_SingleM::getFitness() const {
            // rounds energy and multiplies it with 100, then negate it
        double energy = 100 * rnaFreeEnergy.getEnergy();
        if (energy >= 0)
            return - static_cast<int>(floor(energy + 0.5));
        else
            return - static_cast<int>(ceil(energy - 0.5));
    }

    double  
    S_RNAfe_SingleM::getEnergy(void) const {
        return rnaFreeEnergy.getEnergy();
    }

    unsigned int    
    S_RNAfe_SingleM::getMinimalDistance(
         const State& _state2) const {
        S_RNAfe_SingleM state2 = 
         dynamic_cast<const S_RNAfe_SingleM&>(_state2);

        const biu::Structure& state1Struct = rnaFreeEnergy.getStructureRef();
        const biu::Structure& state2Struct = state2.rnaFreeEnergy.getStructureRef();

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

    S_RNAfe_SingleM* 
    S_RNAfe_SingleM::clone() const {
        return new S_RNAfe_SingleM(*this);
    }
    
    State*
    S_RNAfe_SingleM::fromString(const std::string& stringRep) const {
        std::string backboneRep = ".()";
        size_t pos = stringRep.find_first_not_of(backboneRep);
        std::cout << pos << std::endl;
        const std::string rnaStructBracketNotStr = stringRep.substr(0, pos-1);
        const std::string rnaSeqStr = stringRep.substr(pos, stringRep.size()-pos-1);
        return new S_RNAfe_SingleM(rnaSeqStr, rnaStructBracketNotStr);
    }
    
    std::string 
    S_RNAfe_SingleM::toString()
             const{
        return rnaFreeEnergy.getStringRepresentation();
    }

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

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

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

    S_RNAfe_SingleM::NeighborList::~NeighborList() {
    }

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

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

        SingleMove firstMove = origin->firstValidSingleMove();
        

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

    State* 
    S_RNAfe_SingleM::NeighborList::next(
              State::NeighborList::ItState* itstate, State* elem) const {
        assertbiu(elem != NULL, "Elem is not allowed to be NULL");

        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::INVALID_INDEX)) {
            return NULL;
        }
        else {
            S_RNAfe_SingleM* rnaState =
             dynamic_cast<S_RNAfe_SingleM*>(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::RandomNeighborList::RandomNeighborList(
            const S_RNAfe_SingleM* _origin)
     :  origin(_origin)
    {
    }

    S_RNAfe_SingleM::RandomNeighborList::~RandomNeighborList() {
    }

    void 
    S_RNAfe_SingleM::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::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::RandomNeighborList::first(
              State::NeighborList::ItState** itstate) const 
    {
        if (*itstate != NULL) delete *itstate;

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

        return next(*itstate, rnaState);
    }

    State* 
    S_RNAfe_SingleM::RandomNeighborList::next(
              State::NeighborList::ItState* itstate, State* elem) const 
    {
        assertbiu(elem != NULL, "Elem is not allowed to be NULL");

        ItState* myitstate = static_cast<ItState*>(itstate);
        S_RNAfe_SingleM* rnaState =
         dynamic_cast<S_RNAfe_SingleM*>(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::getRandomNeighbor(State* inPlaceNeigh) const 
    {
        S_RNAfe_SingleM* neigh = NULL;
            // copy this state
        if (inPlaceNeigh != NULL) {
            neigh = dynamic_cast<S_RNAfe_SingleM*>(inPlaceNeigh);
            assertbiu(neigh != NULL, "inPlaceNeigh is no S_RNAfe_SingleM"); 
            (*neigh) = (*this); 
        } else {
            neigh = new S_RNAfe_SingleM(*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::compress(void) const {
        return rnaFreeEnergy.getStructureAlphabet()->compress(rnaFreeEnergy.getStructureRef());
    }
    
        /* 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::compress(CSequence& toFill) const {
        toFill = rnaFreeEnergy.getStructureAlphabet()->compress(rnaFreeEnergy.getStructureRef());
        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::uncompress(const CSequence& cseq, State* toFill) const {

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

        S_RNAfe_SingleM* ret = NULL;
            // use second function
        if (toFill == NULL) {
            ret = this->clone();
        } else {
            ret = dynamic_cast<S_RNAfe_SingleM*>(toFill);
        }
        assertbiu(ret != NULL, "given State toFill is no S_RNAfe_SingleM object.");
            
        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::uncompress(const CSequence& cseq) {
        
        biu::Structure str = rnaFreeEnergy.getStructureAlphabet()->decompress(cseq,rnaFreeEnergy.getStructureRef().size());
        
        if (str.size() != rnaFreeEnergy.getStructureRef().size())
            return NULL;
        
        rnaFreeEnergy.setStructure(str);
        
        return this;
    }
    

} // namespace ell