PolyFunctional.hpp

Go to the documentation of this file.

/*

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/>.

*/

#ifndef __POLY_FUNCTIONAL_HPP__
#define __POLY_FUNCTIONAL_HPP__

#include <Sequence/stateCounter.hpp>
#include <Sequence/PolyTableManip.hpp>
#include <Sequence/typedefs.hpp>
#include <vector>
#include <string>
#include <utility>

namespace Sequence
{
  struct countStates
  {
    template<typename charItr>
    inline stateCounter operator()( charItr beg,
                                    charItr end,
                                    const bool & haveOutgroup = false,
                                    const size_t & outgroup = 0 )
    {
      stateCounter c;
      if(beg>=end || beg+outgroup>=end) return c;
      for( charItr itr=beg ; 
           itr != end ; 
           ++itr )
        {
          if ( (!haveOutgroup) 
               //casting to unsigned is safe since itr is always >= beg
               || ( haveOutgroup && ( unsigned(itr-beg) != outgroup ) ) )
            {
              c(*itr);
            }
        }
      return c;
    }
  };
  
  struct countDerivedStates
  {
    template<typename charItr>
    inline stateCounter operator()( charItr beg,
                                    charItr end,
                                    const bool & haveOutgroup = true,
                                    const size_t & outgroup = 0 )
    {
      stateCounter c;
      if(beg>=end || beg+outgroup>=end || haveOutgroup == false) return c;
      for( charItr itr=beg ; 
           itr != end ; 
           ++itr )
        {
          //casting to unsigned is safe since itr is always >= beg
          if ( (unsigned(itr-beg) != outgroup) 
               && (*itr != *(beg+outgroup))  )
            {
              c(*itr);
            }
        }
      return c;
    }
  };

  struct ssh : public std::binary_function< double,polymorphicSite,double >

  {
    inline double operator()( double & _ssh,
                              const polymorphicSite & site,
                              const bool & haveOutgroup = false,
                              const unsigned & outgroup = 0 ) const
      {
        stateCounter c = countStates()(site.second.begin(),
                                       site.second.end(),
                                       haveOutgroup,outgroup);

        unsigned nsam = site.second.length() - 
          unsigned( haveOutgroup==true ? 1 : 0 ) - c.n;
        double hom=0.;
        if (c.gap==0 && nsam > 1)
          {
            hom += (c.a > 0) ? double(c.a) * double(c.a-1) : 0.;
            hom += (c.g > 0) ? double(c.g) * double(c.g-1) : 0.;
            hom += (c.c > 0) ? double(c.c) * double(c.c-1) : 0.;
            hom += (c.t > 0) ? double(c.t) * double(c.t-1) : 0.;
            hom += (c.zero > 0) ? double(c.zero) * double(c.zero-1) : 0.;
            hom += (c.one > 0) ? double(c.one) * double(c.one-1) : 0.;
          }
        return _ssh += (hom > 0.) ? 1.- (hom/( double(nsam)*double(nsam-1) )) : 0.;
      }
  };

  template<typename counter>
  struct nmuts 
  {
    typedef unsigned result_type;
    inline unsigned operator()(unsigned & nm,
                               const polymorphicSite & site,
                               const bool & haveOutgroup = false,
                               const unsigned & outgroup = 0 ) const
    {
      stateCounter c = counter()(site.second.begin(),
                                 site.second.end(),
                                 haveOutgroup,outgroup);
      unsigned n = 0;
      n += (c.a > 0) ? 1 : 0;
      n += (c.g > 0) ? 1 : 0;
      n += (c.c > 0) ? 1 : 0;
      n += (c.t > 0) ? 1 : 0;
      n += (c.zero > 0) ? 1 : 0;
      n += (c.one > 0) ? 1 : 0;
      return nm += (n>=2) ? n-1 : 0;
    }
  };
}
#endif

---