SeedHandlerMfe.h

Go to the documentation of this file.

 
#ifndef INTARNA_SEEDHANDLERMFE_H_
#define INTARNA_SEEDHANDLERMFE_H_
 
#include "IntaRNA/InteractionEnergy.h"
#include "IntaRNA/SeedConstraint.h"
#include "IntaRNA/SeedHandler.h"
 
#include <boost/multi_array.hpp>
 
#include "IntaRNA/Matrix.h"
 
namespace IntaRNA {
 
class SeedHandlerMfe : public SeedHandler
{
public:
 
    virtual bool guaranteesStackOnlySeeds() const override;
 
 
    typedef boost::multi_array<E_type,5> SeedRecMatrix;
 
    typedef boost::array<SeedRecMatrix::index, 5> SeedIndex;
 
    typedef Matrix< std::pair<E_type, size_t> > SeedMatrix;
 
 
public:
 
    SeedHandlerMfe(
            const InteractionEnergy & energy
            , const SeedConstraint & seedConstraint
            );
 
    virtual ~SeedHandlerMfe() override;
 
    virtual
    size_t
    fillSeed(const size_t i1, const size_t j1, const size_t i2, const size_t j2) override;
 
    virtual
    void
    traceBackSeed( Interaction & interaction, const size_t i1, const size_t i2) const override;
 
 
    virtual
    E_type
    getSeedE( const size_t i1, const size_t i2 ) const override;
 
    virtual
    bool
    isSeedBound( const size_t i1, const size_t i2 ) const override;
 
    virtual
    size_t
    getSeedLength1( const size_t i1, const size_t i2 ) const override;
 
    virtual
    size_t
    getSeedLength2( const size_t i1, const size_t i2 ) const override;
 
 
 
protected:
 
    SeedRecMatrix seedE_rec;
 
    SeedMatrix seed;
 
    size_t offset1;
 
    size_t offset2;
 
    E_type
    getSeedE( const size_t i1, const size_t i2, const size_t bpInbetween, const size_t u1, const size_t u2 ) const;
 
    void
    setSeedE( const size_t i1, const size_t i2, const size_t bpInbetween, const size_t u1, const size_t u2, const E_type E );
 
    size_t
    encodeSeedLength( const size_t l1, const size_t l2 ) const;
 
    size_t
    decodeSeedLength1( const size_t code ) const;
 
    size_t
    decodeSeedLength2( const size_t code ) const;
 
    void
    traceBackSeed( Interaction & interaction
            , const size_t i1, const size_t i2, const size_t bp
            , const size_t u1, const size_t u2 ) const;
 
};
 
 
 
 
inline
SeedHandlerMfe::SeedHandlerMfe(
        const InteractionEnergy & energy
        , const SeedConstraint & seedConstraint
        )
    :
        SeedHandler(energy,seedConstraint)
        , seedE_rec( SeedIndex({{ 0,0,0,0,0 }}))
        , seed()
        , offset1(0)
        , offset2(0)
{
#if INTARNA_IN_DEBUG_MODE
    if ( ! seedConstraint.getExplicitSeeds().empty()) {
        LOG(WARNING) <<"explicit seeds definitions not supported by mfe-seed handler (and thus ignored)";
    }
#endif
}
 
 
inline
SeedHandlerMfe::~SeedHandlerMfe()
{
}
 
 
inline
void
SeedHandlerMfe::
traceBackSeed( Interaction & interaction
        , const size_t i1
        , const size_t i2
        ) const
{
#if INTARNA_IN_DEBUG_MODE
    if ( i1 < offset1 ) throw std::runtime_error("SeedHandlerMfe::traceBackSeed(i1="+toString(i1)+") is out of range (>"+toString(offset1)+")");
    if ( i1-offset1 >= seed.size1() ) throw std::runtime_error("SeedHandlerMfe::traceBackSeed(i1="+toString(i1)+") is out of range (<"+toString(seed.size1()+offset1)+")");
    if ( i2 < offset2 ) throw std::runtime_error("SeedHandlerMfe::traceBackSeed(i2="+toString(i2)+") is out of range (>"+toString(offset2)+")");
    if ( i2-offset2 >= seed.size2() ) throw std::runtime_error("SeedHandlerMfe::traceBackSeed(i2="+toString(i2)+") is out of range (<"+toString(seed.size2()+offset2)+")");
    if ( !( isSeedBound(i1,i2) ) ) throw std::runtime_error("SeedHandlerMfe::traceBackSeed(i1="+toString(i1)+",i2="+toString(i2)+") no seed known (E_INF)");
    if ( i1+getSeedLength1(i1,i2)-1-offset1 >= seed.size1() ) throw std::runtime_error("SeedHandlerMfe::traceBackSeed(i1="+toString(i1)+") seed length ("+toString(getSeedLength1(i1,i2))+") exceeds of range (<"+toString(seed.size1()+offset1)+")");
    if ( i2+getSeedLength2(i1,i2)-1-offset2 >= seed.size2() ) throw std::runtime_error("SeedHandlerMfe::traceBackSeed(i2="+toString(i2)+") seed length ("+toString(getSeedLength2(i1,i2))+") exceeds of range (<"+toString(seed.size2()+offset2)+")");
#endif
 
    // get number of base pairs within the seed
    const size_t seedBps = getConstraint().getBasePairs();
 
    // trace back the according seed
    traceBackSeed( interaction, i1, i2
            , seedBps-2
            , getSeedLength1(i1,i2)-seedBps
            , getSeedLength2(i1,i2)-seedBps );
}
 
 
inline
E_type
SeedHandlerMfe::
getSeedE( const size_t i1, const size_t i2 ) const
{
    return seed(i1-offset1,i2-offset2).first;
}
 
 
inline
bool
SeedHandlerMfe::
isSeedBound( const size_t i1, const size_t i2 ) const
{
    return E_isNotINF( getSeedE(i1,i2) );
}
 
 
inline
size_t
SeedHandlerMfe::
getSeedLength1( const size_t i1, const size_t i2 ) const
{
    return decodeSeedLength1(seed(i1-offset1,i2-offset2).second);
}
 
 
inline
size_t
SeedHandlerMfe::
getSeedLength2( const size_t i1, const size_t i2 ) const
{
    return decodeSeedLength2(seed(i1-offset1,i2-offset2).second);
}
 
 
inline
E_type
SeedHandlerMfe::
getSeedE( const size_t i1, const size_t i2, const size_t bpInbetween, const size_t u1, const size_t u2 ) const
{
#if INTARNA_IN_DEBUG_MODE
    if ( i1 < offset1 ) throw std::runtime_error("SeedHandlerMfe::getSeedE(i1="+toString(i1)+") is out of range (>"+toString(offset1)+")");
    if ( i1-offset1 >= seed.size1() ) throw std::runtime_error("SeedHandlerMfe::getSeedE(i1="+toString(i1)+") is out of range (<"+toString(seed.size1()+offset1)+")");
    if ( i2 < offset2 ) throw std::runtime_error("SeedHandlerMfe::getSeedE(i2="+toString(i2)+") is out of range (>"+toString(offset2)+")");
    if ( i2-offset2 >= seed.size2() ) throw std::runtime_error("SeedHandlerMfe::getSeedE(i2="+toString(i2)+") is out of range (<"+toString(seed.size2()+offset2)+")");
#endif
 
    return seedE_rec( SeedIndex({{
          (SeedRecMatrix::index) (i1-offset1)
        , (SeedRecMatrix::index) (i2-offset2)
        , (SeedRecMatrix::index) bpInbetween
        , (SeedRecMatrix::index) u1
        , (SeedRecMatrix::index) u2 }}) );
}
 
 
inline
void
SeedHandlerMfe::
setSeedE( const size_t i1, const size_t i2, const size_t bpInbetween, const size_t u1, const size_t u2, const E_type E )
{
#if INTARNA_IN_DEBUG_MODE
    if ( i1 < offset1 ) throw std::runtime_error("SeedHandlerMfe::setSeedE(i1="+toString(i1)+") is out of range (>"+toString(offset1)+")");
    if ( i1-offset1 >= seed.size1() ) throw std::runtime_error("SeedHandlerMfe::setSeedE(i1="+toString(i1)+") is out of range (<"+toString(seed.size1()+offset1)+")");
    if ( i2 < offset2 ) throw std::runtime_error("SeedHandlerMfe::setSeedE(i2="+toString(i2)+") is out of range (>"+toString(offset2)+")");
    if ( i2-offset2 >= seed.size2() ) throw std::runtime_error("SeedHandlerMfe::setSeedE(i2="+toString(i2)+") is out of range (<"+toString(seed.size2()+offset2)+")");
#endif
 
    seedE_rec( SeedIndex({{
          (SeedRecMatrix::index) (i1-offset1)
        , (SeedRecMatrix::index) (i2-offset2)
        , (SeedRecMatrix::index) bpInbetween
        , (SeedRecMatrix::index) u1
        , (SeedRecMatrix::index) u2 }}) ) = E;
}
 
 
inline
size_t
SeedHandlerMfe::
encodeSeedLength( const size_t l1, const size_t l2 ) const
{
    return l1 + l2*(seedConstraint.getMaxLength1()+1);
}
 
 
inline
size_t
SeedHandlerMfe::
decodeSeedLength1( const size_t code ) const
{
    return code % (seedConstraint.getMaxLength1()+1);
}
 
 
inline
size_t
SeedHandlerMfe::
decodeSeedLength2( const size_t code ) const
{
    return code / (seedConstraint.getMaxLength1()+1);
}
 
 
inline bool
SeedHandlerMfe::guaranteesStackOnlySeeds() const
{
    return getConstraint().getMaxUnpaired1() == 0 && getConstraint().getMaxUnpaired2() == 0;
}
 
} // namespace
 
#endif /* SEEDHANDLERMFE_H_ */