InteractionEnergyBasePair.h

Go to the documentation of this file.

 
#ifndef INTARNA_INTERACTIONENERGYBASEPAIR_H_
#define INTARNA_INTERACTIONENERGYBASEPAIR_H_
 
#include "IntaRNA/general.h"
#include "IntaRNA/RnaSequence.h"
#include "IntaRNA/Accessibility.h"
#include "IntaRNA/NussinovHandler.h"
#include "IntaRNA/InteractionEnergy.h"
 
 
namespace IntaRNA {
 
 
class InteractionEnergyBasePair: public InteractionEnergy {
 
public:
 
    InteractionEnergyBasePair( const Accessibility & accS1
                    , const ReverseAccessibility & accS2
                    , const size_t maxInternalLoopSize1 = 16
                    , const size_t maxInternalLoopSize2 = 16
                    , const bool initES = false
                    , const Z_type RT = Z_type(1.0)
                    , const E_type bpEnergy = Ekcal_2_E(-1.0)
                    , const size_t minLoopLength = 3
                    , const E_type energyAdd = Ekcal_2_E(0.0)
                    , const bool energyWithDangles = true
                    , const bool internalLoopGU = true
                );
 
    virtual ~InteractionEnergyBasePair();
 
 
    virtual
    E_type
    getES1( const size_t i1, const size_t j1 ) const;
 
    virtual
    E_type
    getES2( const size_t i2, const size_t j2 ) const;
 
    virtual
    E_type
    getE_multiUnpaired( const size_t numUnpaired ) const;
 
 
    virtual
    E_type
    getE_multiHelix( const size_t j1, const size_t j2 ) const;
 
    virtual
    E_type
    getE_multiClosing() const;
 
 
    virtual
    E_type
    getE_init() const;
 
 
    virtual
    E_type
    getE_interLeft( const size_t i1, const size_t j1, const size_t i2, const size_t j2 ) const;
 
    virtual
    E_type
    getE_danglingLeft( const size_t i1, const size_t i2 ) const;
 
 
    virtual
    E_type
    getE_danglingRight( const size_t j1, const size_t j2 ) const;
 
    virtual
    E_type
    getE_endLeft( const size_t i1, const size_t i2 ) const;
 
    virtual
    E_type
    getE_endRight( const size_t j1, const size_t j2 ) const;
 
    virtual
    Z_type
    getRT() const;
 
    E_type
    getE_basePair() const;
 
    virtual
    E_type
    getEall1() const;
 
    virtual
    E_type
    getEall2() const;
 
 
private:
 
    const E_type basePairEnergy;
    const Z_type RT;
    const Z_type basePairWeight;
    const size_t minLoopLength;
    mutable E_type Eall1;
    mutable E_type Eall2;
 
    NussinovHandler::E2dMatrix logQ1;
    NussinovHandler::E2dMatrix logQ2;
 
    /***
     * Compute the ES for a given sequence and store it in the given lookup table.
     * @param seq RNASequence
     * @param logQ The resulting lookuptable
     */
    void computeES(const RnaSequence &seq, NussinovHandler::E2dMatrix &logQ);
 
    E_type
    computeIntraEall( const Accessibility & acc ) const;
 
};
 
 
 
 
inline
InteractionEnergyBasePair::InteractionEnergyBasePair(
        const Accessibility & accS1
        , const ReverseAccessibility & accS2
        , const size_t maxInternalLoopSize1
        , const size_t maxInternalLoopSize2
        , const bool initES
    , const Z_type _RT
    , const E_type bpEnergy
    , const size_t minLoopLen
    , const E_type energyAdd
    , const bool energyWithDangles
    , const bool internalLoopGU
    )
 :
    InteractionEnergy(accS1, accS2
            , maxInternalLoopSize1, maxInternalLoopSize2
            , energyAdd, energyWithDangles,internalLoopGU ),
  RT(_RT),
  basePairEnergy(bpEnergy),
  minLoopLength(minLoopLen),
  basePairWeight(Z_exp(E_2_Z(-bpEnergy) / _RT)),
  Eall1(E_INF),
  Eall2(E_INF),
  logQ1(),
  logQ2()
{
    if (initES) {
     computeES(accS1.getSequence(), logQ1);
     computeES(accS2.getSequence(), logQ2);
    }
}
 
 
inline
InteractionEnergyBasePair::~InteractionEnergyBasePair()
{
}
 
 
inline
E_type
InteractionEnergyBasePair::
getES1( const size_t i1, const size_t j1 ) const
{
#if INTARNA_IN_DEBUG_MODE
    // sanity check
    if (i1>j1) throw std::runtime_error("InteractionEnergy::getES1(i1="+toString(i1)+" > j1="+toString(j1));
    if (j1>=size1()) throw std::runtime_error("InteractionEnergy::getES1() : j1="+toString(j1)+" >= size1()="+toString(size1()));
    if (logQ1.size1() != size1()) throw std::runtime_error("InteractionEnergy::getES1() : ES wasn't computed yet.");
#endif
    return logQ1(i1, j1);
}
 
 
inline
E_type
InteractionEnergyBasePair::
getES2( const size_t i2, const size_t j2 ) const
{
#if INTARNA_IN_DEBUG_MODE
    // sanity check
    if (i2>j2) throw std::runtime_error("InteractionEnergy::getES2(i2="+toString(i2)+" > j2="+toString(j2));
    if (j2>=size2()) throw std::runtime_error("InteractionEnergy::getES2() : j2="+toString(j2)+" >= size2()="+toString(size2()));
    if (logQ2.size1() != size2()) throw std::runtime_error("InteractionEnergy::getES2() : ES wasn't computed yet.");
#endif
    return logQ2(i2, j2);
}
 
 
inline
E_type
InteractionEnergyBasePair::
getE_multiUnpaired( const size_t numUnpaired ) const
{
    return E_type(0);
}
 
 
inline
E_type
InteractionEnergyBasePair::
getE_multiHelix( const size_t j1, const size_t j2 ) const
{
    return E_type(0);
}
 
 
inline
E_type
InteractionEnergyBasePair::
getE_multiClosing() const
{
    return E_type(0);
}
 
 
inline
E_type
InteractionEnergyBasePair::
getE_init() const
{
    return basePairEnergy;
}
 
 
inline
E_type
InteractionEnergyBasePair::
getE_interLeft( const size_t i1, const size_t j1, const size_t i2, const size_t j2 ) const
{
    // if valid internal loop
    if ( isValidInternalLoop(i1,j1,i2,j2) ) {
        // return negated number of gained base pairs by closing this loop = -1
        return basePairEnergy;
    } else {
        return E_INF;
    }
}
 
 
inline
E_type
InteractionEnergyBasePair::
getE_danglingLeft( const size_t i1, const size_t i2 ) const
{
    // no dangling end contribution
    return E_type(0);
}
 
 
inline
E_type
InteractionEnergyBasePair::
getE_danglingRight( const size_t j1, const size_t j2 ) const
{
    // no dangling end contribution
    return E_type(0);
}
 
 
inline
E_type
InteractionEnergyBasePair::
getE_endLeft( const size_t i1, const size_t i2 ) const {
    return E_type(0);
}
 
 
inline
E_type
InteractionEnergyBasePair::
getE_endRight( const size_t j1, const size_t j2 ) const {
    return E_type(0);
}
 
 
inline
E_type
InteractionEnergyBasePair::
getE_basePair() const {
    return basePairEnergy;
}
 
 
inline
Z_type
InteractionEnergyBasePair::
getRT() const {
    return this->RT;
}
 
 
}  // namespace IntaRNA
 
#endif /* INTERACTIONENERGYBASEPAIR_H_ */