main page namespaces classes files
ELL home

src/ell/rna/RNAFreeEnergy.cc

Go to the documentation of this file.

#include <biu/assertbiu.hh>
#include "ell/rna/RNAFreeEnergy.hh"

//#include <ViennaRNA/fold.h>
extern "C" {
    #include <ViennaRNA/fold_vars.h>
}
extern "C"  float energy_of_struct(const char *string, const char *structure);
extern "C"  void update_fold_params ();


#include <cstring>  //For the string functions
#include <limits.h>

namespace ell
{

    const biu::Alphabet RNAFreeEnergy::COMMON_ALPHABET("ACGU", 1);
    const biu::AllowedBasePairs RNAFreeEnergy::COMMON_BPAIRS(
                                &COMMON_ALPHABET, "AU,CG,GU");
    const double RNAFreeEnergy::ENERGY_INF((double)UINT_MAX);
// RNAFreeEnergy

    RNAFreeEnergy::RNAFreeEnergy(const std::string& rnaSeqStr,
                                 const std::string& rnaStructBracketDotStr)
     :  biu::RNAStructure(rnaSeqStr, rnaStructBracketDotStr, &COMMON_BPAIRS),
        energy(ENERGY_INF)
    {
          // ViennaRNA package parameter for dangling ends
        dangles = 2;
        update_fold_params();
          // inform that molecule was initialised
        moleculeChanged();
    }
    
    RNAFreeEnergy::RNAFreeEnergy(biu::Sequence* rnaSeq,
                                 const biu::Structure* const rnaStructBracketDot,
                                 const bool seqIsShared)
     :  biu::RNAStructure(rnaSeq, rnaStructBracketDot, &COMMON_BPAIRS, seqIsShared),
        energy(ENERGY_INF)
    {
          // ViennaRNA package parameter for dangling ends
        dangles = 2;
        update_fold_params();
          // inform that molecule was initialised
        moleculeChanged();
    }
    
    RNAFreeEnergy::RNAFreeEnergy(const RNAFreeEnergy& rnaFreeEnerg)
     : biu::RNAStructure(rnaFreeEnerg), energy(rnaFreeEnerg.energy)
    {
    }

    RNAFreeEnergy::~RNAFreeEnergy() {
    }

    size_t  
    RNAFreeEnergy::nextCompatibleSingleMovePos(const size_t pos) const {
        for (   size_t p = pos+1; 
                p < rnaStructBracketDot->size() && (*rnaStructBracketDot)[p] != STRUCT_BND_CLOSE; 
                (p = rnaBonds[p])++) 
        {
            if ((*rnaStructBracketDot)[p] == STRUCT_UNBOUND) {
                return p;
            }
        }
        return biu::RNAStructure::INVALID_INDEX;
    }

    void        
    RNAFreeEnergy::applySingleMoveInPlace(const SingleMove& move) {
        // open bond between move.first and move.second
        if (rnaBonds[move.first] == move.second) {
            (*rnaStructBracketDot)[move.first] = STRUCT_UNBOUND;
            (*rnaStructBracketDot)[move.second] = STRUCT_UNBOUND;
            rnaBonds[move.first] = biu::RNAStructure::INVALID_INDEX;
        }
        // close bond between pos1 and pos2
        else {
            assertbiu(isValidSingleMove(move),
                "can't close bond with single move, because it is not valid");
            (*rnaStructBracketDot)[move.first] = STRUCT_BND_OPEN;
            (*rnaStructBracketDot)[move.second] = STRUCT_BND_CLOSE;
            rnaBonds[move.first] = move.second;
        }
          // inform that structure has changed
        moleculeChanged();
    }

    void 
    RNAFreeEnergy::setStructure(const biu::Structure& str) {
        assertbiu(str.size() == rnaSeq->size(),
                "given structure is not of sequence size");
          // set structure
        RNAStructure::setStructure(str);
          // update possible bond settings
        initBonds();
          // inform that molecule has changed and e.g. energy has to be 
          // recalculated
        moleculeChanged();
    }

    bool        
    RNAFreeEnergy::isValidSingleMove(const SingleMove& move)
             const {
                
        assertbiu(move.first < move.second && move.second < rnaStructBracketDot->size(),
                "first move pos has to be smaller than second");
                
        // bond between move.first and move.second
        if (rnaBonds[move.first] == move.second) {
            return true;
        }

        // check if both unpaired and not destroying nestedness
        if (    (*rnaStructBracketDot)[move.first] == STRUCT_UNBOUND &&
                (*rnaStructBracketDot)[move.second] == STRUCT_UNBOUND &&
                    // distance between move.first and move.second has to be 
                    // at least biu::RNAStructure::MIN_LOOP_LENGTH
                (move.second-move.first) > biu::RNAStructure::MIN_LOOP_LENGTH &&
                  // check if bases are compatible
                isAllowedBasePair(move.first, move.second) ) 
        {
            size_t i = move.first; 
            while (i < move.second ) {
                if((*rnaStructBracketDot)[i] == STRUCT_BND_CLOSE) {
                    return false;
                } else if((*rnaStructBracketDot)[i] == STRUCT_BND_OPEN) {
                    i = rnaBonds[i]; 
                    i++;
                } else {
                    i++; 
                }
            }
            return i==move.second;
        }
        
        // else no valid structure
        return false;
    }

    double  
    RNAFreeEnergy::getEnergy() const {
        if (energy == ENERGY_INF) {
              // generate input for Vienna package function
//          char * seqStr = new char[rnaSeq->size()+1];
//          strcpy( seqStr, COMMON_ALPHABET.getString(*rnaSeq).c_str()); 
//          char * structStr = new char[rnaStructBracketDot->size()+1];
//          strcpy( structStr, STRUCT_ALPH.getString(*rnaStructBracketDot).c_str()); 
              // calculate energy via Vienna package
//          energy = energy_of_struct( seqStr, structStr );
            energy = energy_of_struct( 
                        COMMON_ALPHABET.getString(*rnaSeq).c_str()
                        , STRUCT_ALPH.getString(*rnaStructBracketDot).c_str() );
              // garbage collection
//          delete seqStr;
//          delete structStr;
        }
        return energy;
    }


    RNAFreeEnergy&  
    RNAFreeEnergy::operator= (const RNAFreeEnergy& rnaStruct2) {
        
        if (this != &rnaStruct2) {
              // call assignment operator of super class
            biu::RNAStructure::operator=(rnaStruct2);
            
              // apply assignments for this class
            energy = rnaStruct2.energy;
        }
        return *this;
    }

    void
    RNAFreeEnergy::moleculeChanged() {
          // force recalculation of energy
        energy = ENERGY_INF;
    }

    
} // namespace ell