HelixHandlerUnpaired.h

Go to the documentation of this file.

 
#ifndef INTARNA_HELIXHANDLERUNPAIRED_H_
#define INTARNA_HELIXHANDLERUNPAIRED_H_
 
#include "IntaRNA/InteractionEnergy.h"
#include "IntaRNA/HelixConstraint.h"
#include "IntaRNA/HelixHandler.h"
 
#include <boost/multi_array.hpp>
 
#include "IntaRNA/Matrix.h"
 
#include <tuple>
 
namespace IntaRNA {
 
class HelixHandlerUnpaired : public HelixHandler {
 
public:
 
    typedef boost::multi_array<std::pair<E_type, size_t>,3> HelixRecMatrix;
 
    typedef boost::array<HelixRecMatrix::index, 3> HelixIndex;
 
    typedef Matrix< std::tuple<E_type, size_t, size_t> > HelixMatrix;
    typedef Matrix< std::pair<E_type, size_t> > HelixSeedMatrix;
 
 
public:
 
    HelixHandlerUnpaired(
            const InteractionEnergy & energy
            , const HelixConstraint & helixConstraint
            , SeedHandler * const seedHandler = NULL
    );
 
    virtual ~HelixHandlerUnpaired();
 
    virtual
    const HelixConstraint&
    getConstraint() const;
 
    virtual
    const InteractionEnergy&
    getInteractionEnergy() const;
 
    virtual
    size_t
    fillHelix( const size_t i1, const size_t j1, const size_t i2, const size_t j2 );
 
    virtual
    size_t
    fillHelixSeed(const size_t i1, const size_t j1, const size_t i2, const size_t j2);
 
    virtual
    void
    traceBackHelix( Interaction & interaction, const size_t i1, const size_t i2);
 
    virtual
    void
    traceBackHelixSeed( Interaction & interaction, const size_t i1, const size_t i2);
 
    virtual
    E_type
    getHelixE( const size_t i1, const size_t i2 ) const;
 
    virtual
    E_type
    getHelixSeedE( const size_t i1, const size_t i2 ) const;
 
    virtual
    size_t
    getHelixLength1( const size_t i1, const size_t i2 ) const;
 
    virtual
    size_t
    getHelixLength2( const size_t i1, const size_t i2 ) const;
 
    virtual
    size_t
    getHelixLength1( const size_t i1, const size_t i2, const size_t bp ) const;
 
    virtual
    size_t
    getHelixLength2( const size_t i1, const size_t i2, const size_t bp ) const;
 
    virtual
    size_t
    getHelixSeedLength1( const size_t i1, const size_t i2 ) const;
 
    virtual
    size_t
    getHelixSeedLength2( const size_t i1, const size_t i2 ) const;
 
    void setSeedHandler(SeedHandler & seedHandler);
 
protected:
 
    size_t
    getBestBP( const size_t i1, const size_t i2) const;
 
    E_type
    getHelixE( const size_t i1, const size_t i2, const size_t bp);
 
 
    void
    setHelixPair( const size_t i1, const size_t i2, const size_t bp
                , const E_type E, const size_t length );
 
 
    size_t
    encodeHelixLength( const size_t l1, const size_t l2 ) const;
 
    size_t
    decodeHelixLength1( const size_t code ) const;
 
 
    size_t
    decodeHelixLength2( const size_t code ) const;
 
    size_t
    encodeHelixSeedLength( const size_t l1, const size_t l2 ) const;
 
    size_t
    decodeHelixSeedLength1( const size_t code ) const;
 
    size_t
    decodeHelixSeedLength2( const size_t code ) const;
 
    void
    traceBackHelix( Interaction & interaction
            , const size_t i1, const size_t i2, const size_t bp);
 
protected:
 
    const InteractionEnergy & energy;
 
    const HelixConstraint & helixConstraint;
 
    HelixRecMatrix helixE_rec;
 
    HelixMatrix helix;
 
    HelixSeedMatrix helixSeed;
 
    size_t offset1;
 
    size_t offset2;
 
    SeedHandler * seedHandler;
};
 
 
 
inline
HelixHandlerUnpaired::HelixHandlerUnpaired(
        const InteractionEnergy & energy
        , const HelixConstraint & helixConstraint
        , SeedHandler * const seedHandler
)
        :
        energy(energy)
        , helixConstraint(helixConstraint)
        , seedHandler(seedHandler)
        , helixE_rec( HelixIndex({{ 0, 0, 0 }}))
        , helix()
        , helixSeed()
        , offset1(0)
        , offset2(0)
{
    if (seedHandler != NULL) {
        setSeedHandler(*seedHandler);
    }
}
 
 
inline
HelixHandlerUnpaired::~HelixHandlerUnpaired()
{
}
 
 
inline
const InteractionEnergy&
HelixHandlerUnpaired::
getInteractionEnergy() const
{
    return energy;
}
 
 
inline
const HelixConstraint&
HelixHandlerUnpaired::
getConstraint() const
{
    return helixConstraint;
}
 
inline
size_t
HelixHandlerUnpaired::
getBestBP(const size_t i1, const size_t i2) const
{
    return std::get<2>(helix(i1-offset1, i2-offset2));
}
 
 
inline
void
HelixHandlerUnpaired::traceBackHelix(Interaction &interaction
        , const size_t i1
        , const size_t i2
)
{
#if INTARNA_IN_DEBUG_MODE
    if ( i1 < offset1 ) throw std::runtime_error("HelixHandlerUnpaired::traceBackHelix(i1="+toString(i1)+") is out of range (>"+toString(offset1)+")");
    if ( i1-offset1 >= helix.size1() ) throw std::runtime_error("HelixHandlerUnpaired::traceBackHelix(i1="+toString(i1)+") is out of range (<"+toString(helix.size1()+offset1)+")");
    if ( i2 < offset2 ) throw std::runtime_error("HelixHandlerUnpaired::traceBackHelix(i2="+toString(i2)+") is out of range (>"+toString(offset2)+")");
    if ( i2-offset2 >= helix.size2() ) throw std::runtime_error("HelixHandlerUnpaired::traceBackHelix(i2="+toString(i2)+") is out of range (<"+toString(helix.size2()+offset2)+")");
    if ( E_isINF( getHelixE(i1,i2) ) ) throw std::runtime_error("HelixHandlerUnpaired::traceBackHelix(i1="+toString(i1)+",i2="+toString(i2)+") no helix known (E_INF)");
    if ( i1+getHelixLength1(i1,i2)-1-offset1 >= helix.size1() ) throw std::runtime_error("HelixHandlerUnpaired::traceBackHelix(i1="+toString(i1)+") helix length ("+toString(getHelixLength1(i1,i2))+") exceeds of range (<"+toString(helix.size1()+offset1)+")");
    if ( i2+getHelixLength2(i1,i2)-1-offset2 >= helix.size2() ) throw std::runtime_error("HelixHandlerUnpaired::traceBackHelix(i2="+toString(i2)+") helix length ("+toString(getHelixLength2(i1,i2))+") exceeds of range (<"+toString(helix.size2()+offset2)+")");
#endif
    // get number of base pairs allowed within the helix
    const size_t bestBP = getBestBP(i1,i2);
 
    // trace back the according helix
    traceBackHelix( interaction, i1-offset1, i2-offset2, bestBP);
}
 
 
inline
E_type
HelixHandlerUnpaired::
getHelixE(const size_t i1, const size_t i2) const
{
    return std::get<0>(helix(i1-offset1, i2-offset2));
}
 
 
 
inline
E_type
HelixHandlerUnpaired::
getHelixSeedE(const size_t i1, const size_t i2) const
{
    return helixSeed(i1-offset1, i2-offset2).first;
}
 
 
inline
size_t
HelixHandlerUnpaired::
getHelixLength1(const size_t i1, const size_t i2) const
{
    return decodeHelixLength1(std::get<1>(helix(i1-offset1, i2-offset2)));
}
 
 
inline
size_t
HelixHandlerUnpaired::
getHelixLength2(const size_t i1, const size_t i2) const
{
    return decodeHelixLength2(std::get<1>(helix(i1-offset1, i2-offset2)));
}
 
 
inline
size_t
HelixHandlerUnpaired::
getHelixLength1(const size_t i1, const size_t i2, const size_t bp) const {
// if no base pair is given return 0 in order to simplify helixHandlerSeed computation.
    if (bp <= 1) {
        return 0;
    } else {
        return decodeHelixLength1(helixE_rec(HelixIndex(
                {{(HelixRecMatrix::index) i1, (HelixRecMatrix::index) i2, (HelixRecMatrix::index) bp}})).second);
    }
}
 
inline
size_t
HelixHandlerUnpaired::
getHelixLength2(const size_t i1, const size_t i2, const size_t bp) const
{
    // if no base pair is given return 0 in order to simplify helixHandlerSeed computation.
    if (bp <= 1) {
        return 0;
    } else {
        return decodeHelixLength2(helixE_rec(HelixIndex(
                {{(HelixRecMatrix::index) i1, (HelixRecMatrix::index) i2, (HelixRecMatrix::index) bp}})).second);
    }
}
 
 
inline
size_t
HelixHandlerUnpaired::
getHelixSeedLength1(const size_t i1, const size_t i2) const
{
    return decodeHelixSeedLength1(helixSeed(i1-offset1, i2-offset2).second);
}
 
 
inline
size_t
HelixHandlerUnpaired::
getHelixSeedLength2(const size_t i1, const size_t i2) const
{
    return decodeHelixSeedLength2(helixSeed(i1-offset1, i2-offset2).second);
}
 
 
inline
E_type
HelixHandlerUnpaired::
getHelixE(const size_t i1, const size_t i2, const size_t bp)
{
    // if no base pair is given return 0 in order to simplify helixHandlerSeed computation.
    if (bp <= 1) {
        return 0;
    } else {
        return helixE_rec(HelixIndex({{(HelixRecMatrix::index) i1, (HelixRecMatrix::index) i2, (HelixRecMatrix::index) bp}})).first;
    }
}
 
 
inline
void
HelixHandlerUnpaired::
setHelixPair(const size_t i1, const size_t i2, const size_t bp, const E_type E, const size_t length)
{
    helixE_rec( HelixIndex({{
            (HelixRecMatrix::index) i1
            , (HelixRecMatrix::index) i2
            , (HelixRecMatrix::index) bp}}) ) = std::make_pair(E, E_isINF(E) ? 0 : length);
}
 
 
inline
size_t
HelixHandlerUnpaired::
encodeHelixLength(const size_t l1, const size_t l2) const
{
    return l1 + l2 * (helixConstraint.getMaxLength1()+1);
}
 
 
inline
size_t
HelixHandlerUnpaired::
decodeHelixLength1(const size_t code) const
{
    return code % (helixConstraint.getMaxLength1()+1);
}
 
 
inline
size_t
HelixHandlerUnpaired::
decodeHelixLength2(const size_t code) const
{
    return code / (helixConstraint.getMaxLength1()+1);
}
 
 
inline
size_t
HelixHandlerUnpaired::
encodeHelixSeedLength(const size_t l1, const size_t l2) const
{
    return l1 + l2 * (helixConstraint.getMaxLength1()+ seedHandler->getConstraint().getMaxLength1()+1);
}
 
 
inline
size_t
HelixHandlerUnpaired::
decodeHelixSeedLength1(const size_t code) const
{
    return code % (helixConstraint.getMaxLength1()+ seedHandler->getConstraint().getMaxLength1()+1);
}
 
 
inline
size_t
HelixHandlerUnpaired::
decodeHelixSeedLength2(const size_t code) const
{
    return code / (helixConstraint.getMaxLength1()+ seedHandler->getConstraint().getMaxLength1()+1);
}
 
 
inline
void
HelixHandlerUnpaired::
setSeedHandler(SeedHandler & seedHandler) {
    this->seedHandler = &seedHandler;
}
 
 
 
 
} // namespace
 
#endif //INTARNA_HELIXHANDLERUNPAIRED_H_