TwoSubs.cc

/*

  Copyright (C) 2003-2009 Kevin Thornton, krthornt[]@[]uci.edu

  Remove the brackets to email me.

  This file is part of libsequence.

  libsequence is free software: you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation, either version 3 of the License, or
  (at your option) any later version.

  libsequence is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.

  You should have received a copy of the GNU General Public License
  long with libsequence.  If not, see <http://www.gnu.org/licenses/>.

*/

#include <cmath>
#include <cfloat>
#include <Sequence/SingleSub.hpp>
#include <Sequence/RedundancyCom95.hpp>
#include <Sequence/WeightingSchemes.hpp>
#include <Sequence/PathwayHelper.hpp>
#include <Sequence/TwoSubs.hpp>

namespace Sequence
{
  struct TwoSubs::TwoSubsImpl
  {
    double p0, p2S, p2V, p4, q0, q2S, q2V,q4;
    double p0_b1, p2S_b1, p2V_b1, p4_b1, q0_b1, q2S_b1, q2V_b1, q4_b1;
    double p0_b2, p2S_b2, p2V_b2, p4_b2, q0_b2, q2S_b2, q2V_b2, q4_b2;
    double p0_b3, p2S_b3, p2V_b3, p4_b3, q0_b3, q2S_b3, q2V_b3, q4_b3;
    double p0_b4, p2S_b4, p2V_b4, p4_b4, q0_b4, q2S_b4, q2V_b4, q4_b4;
    void Calculate (const RedundancyCom95 & sitesObj, const std::string &codon1,
                    const std::string &int_1, const std::string &int_2,
                    const std::string &codon2, const double w_path1,
                    const double w_path2);
    TwoSubsImpl(void) : p0(0.), p2S(0.), p2V(0.), p4(0.), q0(0.), q2S(0.), q2V(0.),q4(0.),
                        p0_b1(0.), p2S_b1(0.), p2V_b1(0.), p4_b1(0.), q0_b1(0.), q2S_b1(0.), q2V_b1(0.), q4_b1(0.),
                        p0_b2(0.), p2S_b2(0.), p2V_b2(0.), p4_b2(0.), q0_b2(0.), q2S_b2(0.), q2V_b2(0.), q4_b2(0.),
                        p0_b3(0.), p2S_b3(0.), p2V_b3(0.), p4_b3(0.), q0_b3(0.), q2S_b3(0.), q2V_b3(0.), q4_b3(0.),
                        p0_b4(0.), p2S_b4(0.), p2V_b4(0.), p4_b4(0.), q0_b4(0.), q2S_b4(0.), q2V_b4(0.), q4_b4(0.)
    {
    }
  };
  
  TwoSubs::TwoSubs() : impl(std::unique_ptr<TwoSubsImpl>(new TwoSubsImpl()))
  {
  }

  void TwoSubs::operator() (const RedundancyCom95 & sitesObj,
                            const std::string & codon1, const std::string & codon2,
                            const Sequence::WeightingScheme2 *weights2)
  {
    impl->p0= impl->p2S= impl->p2V= impl->p4= impl->q0= impl->q2S= impl->q2V=impl->q4=0.;
    impl->p0_b1= impl->p2S_b1= impl->p2V_b1= impl->p4_b1= impl->q0_b1= impl->q2S_b1= impl->q2V_b1= impl->q4_b1=0.;
    impl->p0_b2= impl->p2S_b2= impl->p2V_b2= impl->p4_b2= impl->q0_b2= impl->q2S_b2= impl->q2V_b2= impl->q4_b2=0.;
    impl->p0_b3= impl->p2S_b3= impl->p2V_b3= impl->p4_b3= impl->q0_b3= impl->q2S_b3= impl->q2V_b3= impl->q4_b3=0.;
    impl->p0_b4= impl->p2S_b4= impl->p2V_b4= impl->p4_b4= impl->q0_b4= impl->q2S_b4= impl->q2V_b4= impl->q4_b4=0.;
    auto intermediates = Intermediates2(codon1,codon2);
    auto weights = weights2->operator()(codon1,codon2,sitesObj.gencode());
    impl->Calculate (sitesObj, codon1, intermediates[0], codon2, intermediates[1], weights[0], weights[1]);
  }

  TwoSubs::~TwoSubs (void)
  {}

  void
  TwoSubs::TwoSubsImpl::Calculate (const RedundancyCom95 & sitesObj, const std::string & codon1,
                          const std::string & int_1, const std::string & int_2,
                          const std::string & codon2, const double w_path1,
                          const double w_path2)
  {
    //count the changes along each branch
    //branches that go through stop codons are counted as
    //containing no changes on them.
    SingleSub Single;
    Single(sitesObj,codon1,int_1);
    p0_b1 = Single.P0();
    p2S_b1 =Single.P2S();
    p2V_b1 =Single.P2V();
    p4_b1 = Single.P4();
    q0_b1 = Single.Q0();
    q2S_b1 =Single.Q2S();
    q2V_b1 =Single.Q2V();
    q4_b1 = Single.Q4();

    Single (sitesObj, int_1, codon2);
    p0_b2 = Single.P0();
    p2S_b2 =Single.P2S();
    p2V_b2 =Single.P2V();
    p4_b2 = Single.P4();
    q0_b2 = Single.Q0();
    q2S_b2 =Single.Q2S();
    q2V_b2 =Single.Q2V();
    q4_b2 = Single.Q4();


    Single (sitesObj, codon1, int_2);
    p0_b3 = Single.P0();
    p2S_b3 =Single.P2S();
    p2V_b3 =Single.P2V();
    p4_b3 = Single.P4();
    q0_b3 = Single.Q0();
    q2S_b3 =Single.Q2S();
    q2V_b3 =Single.Q2V();
    q4_b3 = Single.Q4();

    Single (sitesObj, int_2, codon2);
    p0_b4 = Single.P0();
    p2S_b4 =Single.P2S();
    p2V_b4 =Single.P2V();
    p4_b4 = Single.P4();
    q0_b4 = Single.Q0();
    q2S_b4 =Single.Q2S();
    q2V_b4 =Single.Q2V();
    q4_b4 = Single.Q4();


    p0  = (p0_b1  + p0_b2)  * w_path1
      + (p0_b3  + p0_b4) * w_path2;
    p2S = (p2S_b1 + p2S_b2) * w_path1
      + (p2S_b3 + p2S_b4) * w_path2;
    p2V = (p2V_b1 + p2V_b2) * w_path1
      + (p2V_b3 + p2V_b4) * w_path2;
    p4  = (p4_b1  + p4_b2)  * w_path1
      + (p4_b3  + p4_b4) * w_path2;
    q0  = (q0_b1  + q0_b2)  * w_path1
      + (q0_b3  + q0_b4) * w_path2;
    q2S = (q2S_b1 + q2S_b2) * w_path1
      + (q2S_b3 + q2S_b4) * w_path2;
    q2V = (q2V_b1 + q2V_b2) * w_path1
      + (q2V_b3 + q2V_b4) * w_path2;
    q4  = (q4_b1  + q4_b2)  * w_path1
      + (q4_b3  + q4_b4) * w_path2;
  }



  double
  TwoSubs::P0 (void) const
  {
    return impl->p0;
  }

  double
  TwoSubs::P2S (void) const
  {
    return impl->p2S;
  }

  double
  TwoSubs::P2V (void) const
  {
    return impl->p2V;
  }

  double
  TwoSubs::P4 (void) const
  {
    return impl->p4;
  }

  double
  TwoSubs::Q0 (void) const
  {
    return impl->q0;
  }

  double
  TwoSubs::Q2S (void) const
  {
    return impl->q2S;
  }

  double
  TwoSubs::Q2V (void) const
  {
    return impl->q2V;
  }

  double
  TwoSubs::Q4 (void) const
  {
    return impl->q4;
  }

}