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

  TwoSubs::~TwoSubs (void)
  {}

  void
  TwoSubs::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 p0;
  }

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

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

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

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

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

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

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

}

---