SeedHandlerNoBulge.h

Go to the documentation of this file.

 
#ifndef INTARNA_SEEDHANDLERNOBULGE_H_
#define INTARNA_SEEDHANDLERNOBULGE_H_
 
#include "IntaRNA/InteractionEnergy.h"
#include "IntaRNA/SeedConstraint.h"
#include "IntaRNA/SeedHandler.h"
 
#include <boost/unordered_map.hpp>
 
 
namespace IntaRNA {
 
class SeedHandlerNoBulge : public SeedHandler
{
 
protected:
 
    typedef std::vector<E_type> StackingEnergyList;
 
    typedef boost::unordered_map< Interaction::BasePair, E_type, Interaction::BasePair::Hash > SeedHash;
 
    SeedHash seedForLeftEnd;
 
 
public:
 
    SeedHandlerNoBulge(
            const InteractionEnergy & energy
            , const SeedConstraint & seedConstraint
            );
 
    virtual ~SeedHandlerNoBulge();
 
    virtual
    size_t
    fillSeed(const size_t i1, const size_t j1, const size_t i2, const size_t j2);
 
    virtual
    void
    traceBackSeed( Interaction & interaction, const size_t i1, const size_t i2) const;
 
 
    virtual
    E_type
    getSeedE( const size_t i1, const size_t i2 ) const;
 
    virtual
    bool
    isSeedBound( const size_t i1, const size_t i2 ) const;
 
    virtual
    size_t
    getSeedLength1( const size_t i1, const size_t i2 ) const;
 
    virtual
    size_t
    getSeedLength2( const size_t i1, const size_t i2 ) const;
 
    virtual
    bool
    updateToNextSeed( size_t & i1, size_t & i2
            , const size_t i1min = 0, const size_t i1max = RnaSequence::lastPos
            , const size_t i2min = 0, const size_t i2max = RnaSequence::lastPos
            ) const;
 
 
protected:
 
    void
    storeSeed( const size_t j1, const size_t j2, const StackingEnergyList & bpE );
 
};
 
 
 
 
inline
SeedHandlerNoBulge::SeedHandlerNoBulge(
        const InteractionEnergy & energy
        , const SeedConstraint & seedConstraint
        )
    :
        SeedHandler(energy,seedConstraint)
        , seedForLeftEnd()
{
#if INTARNA_IN_DEBUG_MODE
    if ( ! seedConstraint.getExplicitSeeds().empty()) {
        LOG(WARNING) <<"explicit seeds definitions not supported by noBulge-seed handler (and thus ignored)";
    }
    if ( seedConstraint.getMaxUnpaired1() > 0 ) { throw std::runtime_error("SeedHandlerNoBulge() : seedConstraint.maxUnpaired1 > 0"); }
    if ( seedConstraint.getMaxUnpaired2() > 0 ) { throw std::runtime_error("SeedHandlerNoBulge() : seedConstraint.maxUnpaired2 > 0"); }
#endif
}
 
 
inline
SeedHandlerNoBulge::~SeedHandlerNoBulge()
{
    seedForLeftEnd.clear();
}
 
 
inline
void
SeedHandlerNoBulge::
traceBackSeed( Interaction & interaction
        , const size_t i1
        , const size_t i2
        ) const
{
#if INTARNA_IN_DEBUG_MODE
    if ( !( isSeedBound(i1,i2) ) ) throw std::runtime_error("SeedHandlerNoBulge::traceBackSeed(i1="+toString(i1)+",i2="+toString(i2)+") no seed known");
#endif
 
    // get number of base pairs within the seed
    const size_t seedBps = getConstraint().getBasePairs();
 
    // add seed base pairs
    // but exclude left- and right-most bp
    for (size_t k=1; (k+1)<seedBps; k++) {
        interaction.basePairs.push_back( energy.getBasePair(i1+k,i2+k) );
    }
 
}
 
 
inline
E_type
SeedHandlerNoBulge::
getSeedE( const size_t i1, const size_t i2 ) const
{
    // search for seed entry in hash
    SeedHash::const_iterator seedData = seedForLeftEnd.find( SeedHash::key_type(i1,i2) );
    // return according energy value
    return seedData == seedForLeftEnd.end() ? E_INF : seedData->second;
}
 
 
inline
bool
SeedHandlerNoBulge::
isSeedBound( const size_t i1, const size_t i2 ) const
{
    // search for seed entry in hash
    return seedForLeftEnd.find( SeedHash::key_type(i1,i2) ) != seedForLeftEnd.end();
}
 
 
inline
size_t
SeedHandlerNoBulge::
getSeedLength1( const size_t i1, const size_t i2 ) const
{
#if INTARNA_IN_DEBUG_MODE
    if ( !( isSeedBound(i1,i2) ) ) throw std::runtime_error("SeedHandlerNoBulge::getSeedLength1(i1="+toString(i1)+",i2="+toString(i2)+") no seed known");
#endif
    // check if
    return getConstraint().getBasePairs();
}
 
 
inline
size_t
SeedHandlerNoBulge::
getSeedLength2( const size_t i1, const size_t i2 ) const
{
#if INTARNA_IN_DEBUG_MODE
    if ( !( isSeedBound(i1,i2) ) ) throw std::runtime_error("SeedHandlerNoBulge::getSeedLength2(i1="+toString(i1)+",i2="+toString(i2)+") no seed known");
#endif
    // check if
    return getConstraint().getBasePairs();
}
 
 
inline
void
SeedHandlerNoBulge::
storeSeed( const size_t j1, const size_t j2, const StackingEnergyList & bpE )
{
    const size_t seedBP = seedConstraint.getBasePairs();
    E_type seedEhybrid = 0;
 
    // check if EDs of full seed are within boundaries
    if ( energy.getED1(j1+1-seedBP, j1) < seedConstraint.getMaxED()
        && energy.getED2(j2+1-seedBP, j2) < seedConstraint.getMaxED() )
    {
        // check seed boundaries if needed
        if ( ! seedConstraint.isGUendAllowed()
            && (energy.isGU(j1,j2) || energy.isGU(j1+1-seedBP,j2+1-seedBP)) )
        {
            return;
        }
 
        // compute seed energy
        for(auto e=bpE.begin(); e!=bpE.end(); e++) { seedEhybrid += *e; } // (left) stacking energies
        // check hybridization energy bound (incl E_init)
        if ( (seedEhybrid+energy.getE_init()) >= seedConstraint.getMaxEhybrid()) {
            return;
        }
 
        E_type seedEfull = energy.getE(j1+1-seedBP, j1, j2+1-seedBP, j2, seedEhybrid) + energy.getE_init(); // ED values etc.
 
        // store seed hybridization energy if overall E is below to threshold
        if( seedEfull < seedConstraint.getMaxE() ) {
            seedForLeftEnd[ SeedHash::key_type(j1+1-seedBP,j2+1-seedBP) ] = seedEhybrid ;
        }
    }
 
 
}
 
 
} // namespace
 
#endif /* SEEDHANDLERNOBULGE_H_ */