PCR in C++

main.cpp

#include "pcr.hpp"

#include <iostream>

int main(int argc, char * argv [])
{
    DNA dna({DNA::Base::A, DNA::Base::C, DNA::Base::G, DNA::Base::T, DNA::Base::C, DNA::Base::T});
    
    DNA::container_t ret = pcr(dna, 1, 3, 10);
    
    std::size_t correct = 0, incorrect = 0;
    
    for (auto itr = ret.begin(), end = ret.end();
         itr != end;
         ++itr)
    {
        if (itr->size() == 3) ++correct;
        
        else ++incorrect;
    }
    
    std::cout << "Correct: " << correct << std::endl;
    
    std::cout << "Incorrect: " << incorrect << std::endl;
}

pcr.cpp

#include "pcr.hpp"

#include <stdexcept>


DNA::DNA() = default;

DNA::DNA(const sequence_t& left)
{
    setLeftStrand(left);
    setRightStrand(getLeftComplementary(left.size() - 1));
}

DNA::DNA(const sequence_t& left,
         const sequence_t& right)
{
    setLeftStrand(left);
    setRightStrand(right);
}

DNA::DNA(const strand_t& left)
: left_(left),
  right_(getLeftComplementary(left_.size() - 1))
{ }

DNA::DNA(const strand_t& left,
         const strand_t& right)
: left_(left),
right_(right)
{ }

void DNA::setLeftStrand(const strand_t& strand)
{
    left_ = strand;
}

void DNA::setLeftStrand(const std::vector<Base>& bases)
{
    for (locus_t i = 0, end = bases.size(); i < end; ++i)
    {
        left_.push_back({bases[i], i});
    }
}

DNA::strand_t DNA::getLeftStrand()
{
    return left_;
}

const DNA::strand_t DNA::getLeftStrand() const
{
    return left_;
}

DNA::left_iterator DNA::getLeftFivePrime()
{
    return left_.begin();
}

DNA::const_left_iterator DNA::getLeftFivePrime() const
{
    return left_.begin();
}

DNA::left_iterator DNA::getLeftThreePrime()
{
    return left_.end();
}

DNA::const_left_iterator DNA::getLeftThreePrime() const
{
    return left_.end();
}

void DNA::appendToLeftStrand(DNA::Base base)
{
    left_.push_back({base, (getLeftThreePrime() - 1)->locus});
}

DNA::strand_t DNA::getLeftComplementary(locus_t primer) const
{
    strand_t complementary;
    
    const_left_iterator itr = std::find_if(getLeftFivePrime(),
                                           getLeftThreePrime(),
                                           [&] (const Nucleotide& n)
                                           { return n.locus == primer; });
    
    if (itr == getLeftThreePrime())
    {
        throw std::invalid_argument("Could not find primer locus in left strand!");
    }
    
    for(const_left_iterator begin = getLeftFivePrime();
        itr >= begin;
        --itr)
    {
        Base base;
        
        switch(itr->base)
        {
            case DNA::Base::A: base = DNA::Base::T; break;
                
            case DNA::Base::C: base = DNA::Base::G; break;
                
            case DNA::Base::G: base = DNA::Base::C; break;
                
            case DNA::Base::T: base = DNA::Base::A; break;
        }
        
        complementary.push_front({base, itr->locus});
    }
    
    return complementary;
}

void DNA::setRightStrand(const strand_t &strand)
{
    right_ = strand;
}

void DNA::setRightStrand(const std::vector<Base>& bases)
{
    for (locus_t i = 0, end = bases.size(); i < end; ++i)
    {
        right_.push_back({bases[i], i});
    }
}

DNA::strand_t DNA::getRightStrand()
{
    return right_;
}

const DNA::strand_t DNA::getRightStrand() const
{
    return right_;
}

DNA::right_iterator DNA::getRightFivePrime()
{
    return right_.rbegin();
}

DNA::const_right_iterator DNA::getRightFivePrime() const
{
    return right_.rbegin();
}

DNA::right_iterator DNA::getRightThreePrime()
{
    return right_.rend();
}

DNA::const_right_iterator DNA::getRightThreePrime() const
{
    return right_.rend();
}

void DNA::appendToRightStrand(DNA::Base base)
{
    right_.push_front({base, (getRightFivePrime() - 1)->locus});
}

DNA::strand_t DNA::getRightComplementary(locus_t primer) const
{
    strand_t complementary;
    
    const_right_iterator itr = std::find_if(getRightFivePrime(),
                                            getRightThreePrime(),
                                            [&] (const Nucleotide& n)
                                            { return n.locus == primer; });
    
    if (itr == getRightThreePrime())
    {
        throw std::invalid_argument("Could not find primer locus in right strand!");
    }
    
    for(const_right_iterator begin = getRightFivePrime();
        itr >= begin;
        --itr)
    {
        Base base;
        
        switch(itr->base)
        {
            case DNA::Base::A: base = DNA::Base::T; break;
                
            case DNA::Base::C: base = DNA::Base::G; break;
                
            case DNA::Base::G: base = DNA::Base::C; break;
                
            case DNA::Base::T: base = DNA::Base::A; break;
        }
        
        complementary.push_back({base, itr->locus});
    }
    
    return complementary;
}

DNA::locus_t DNA::leftToRightLocus(locus_t leftLocus) const
{
    return right_.size() - leftLocus - 1;
}

DNA::locus_t DNA::rightToLeftLocus(locus_t rightLocus) const
{
    return left_.size() - rightLocus - 1;
}

DNA::container_t pcr(const DNA& dna,
                     DNA::locus_t first,
                     DNA::locus_t last,
                     std::size_t n)
{
    if (! n) return {dna.getLeftStrand(), dna.getRightStrand()};
    
    else
    {
        DNA left {dna.getLeftStrand(), dna.getLeftComplementary(last)};
        
        DNA::container_t&& leftReturned = pcr(left, first, last, n - 1);
        
        
        DNA right {dna.getRightComplementary(first), dna.getRightStrand()};
        
        DNA::container_t&& rightReturned = pcr(right, first, last, n - 1);
        
        
        leftReturned.insert(leftReturned.end(), rightReturned.begin(), rightReturned.end());
        
        return leftReturned;
    }
}

pcr.hpp

#ifndef pcr_hpp
#define pcr_hpp

#include <vector>
#include <deque>

class DNA
{
    
public:
    
    enum class Base { A, C, G, T };
    
    typedef std::size_t locus_t;
    
    struct Nucleotide { Base base; locus_t locus; };
    
    typedef std::vector<Base> sequence_t;
    
    typedef std::deque<Nucleotide> strand_t;
    
    typedef std::vector<strand_t> container_t;
    

    typedef strand_t::iterator left_iterator;
    
    typedef strand_t::const_iterator const_left_iterator;
    

    typedef strand_t::reverse_iterator right_iterator;
    
    typedef strand_t::const_reverse_iterator const_right_iterator;
    
    
    DNA();
    
    DNA(const sequence_t& left);
    
    DNA(const sequence_t& left,
        const sequence_t& right);
    
    DNA(const strand_t& left);
    
    DNA(const strand_t& left,
        const strand_t& right);
    
    
    void setLeftStrand(const strand_t& strand);
    
    void setLeftStrand(const sequence_t& bases);
    
    
    strand_t getLeftStrand();
    
    const strand_t getLeftStrand() const;
    
    
    left_iterator getLeftFivePrime();
    
    left_iterator getLeftThreePrime();
    
    
    const_left_iterator getLeftFivePrime() const;
    
    const_left_iterator getLeftThreePrime() const;
    
    
    void appendToLeftStrand(Base base);
    
    strand_t getLeftComplementary(locus_t primer = 0) const;
    
    
    void setRightStrand(const strand_t& strand);
    
    void setRightStrand(const sequence_t& bases);
    
    
    strand_t getRightStrand();
    
    const strand_t getRightStrand() const;
    
    
    right_iterator getRightFivePrime();
    
    right_iterator getRightThreePrime();
    
    
    const_right_iterator getRightFivePrime() const;
    
    const_right_iterator getRightThreePrime() const;
    
    
    void appendToRightStrand(Base base);
    
    strand_t getRightComplementary(locus_t primer = 0) const;
    
    
    locus_t rightToLeftLocus(locus_t rightLocus) const;
    
    locus_t leftToRightLocus(locus_t leftLocus) const;
    
private:
    
    strand_t left_;
    
    strand_t right_;
};

DNA::container_t pcr(const DNA& dna,
                     DNA::locus_t first,
                     DNA::locus_t last,
                     std::size_t n);

#endif /* pcr_hpp */