LinkageDisequilibrium.hpp

/*
    Copyright 2009 Stéphane De Mita, Mathieu Siol

    This file is part of the EggLib library.

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

    EggLib 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
    along with EggLib.  If not, see <http://www.gnu.org/licenses/>.
*/

#ifndef EGGLIB_LINKAGEDISEQUILIBRUM_HPP
#define EGGLIB_LINKAGEDISEQUILIBRUM_HPP


#include "BaseDiversity.hpp"
#include "EggException.hpp"


namespace egglib {

    class LinkageDisequilibrium : public BaseDiversity {

      public:

        LinkageDisequilibrium();

        virtual ~LinkageDisequilibrium();

        void load(CharMatrix& data,
                double minimumExploitableData=1.,
                unsigned int ignoreFrequency=0,
                std::string characterMapping=dnaMapping);
    
    
        unsigned int numberOfPairs() const;

        int d(unsigned int pair_index);

        double D(unsigned int pair_index);

        double Dp(unsigned int pair_index);

        double r(unsigned int pair_index);

        double r2(unsigned int pair_index);

        unsigned int site1(unsigned int pair_index);

        unsigned int site2(unsigned int pair_index);

        double correl() const;
        
        unsigned int Rmin(CharMatrix& data) const;



      protected:
      
        // adds a pair of polymorphic sites
        // assume position2>position1,
        //  sites are polymorphic with exactly 2 alleles
        void add(CharMatrix& data, unsigned int position1, unsigned int position2);

        // Constructor help
        void init();
                
        // Destructor helper
        void clear();
        
        // Resizes arrays
        void reset();
        
        // Small helper
        inline double min(double a, double b) { return (a>b)?a:b;}

        // Small helper
        inline double max(double a, double b) { return (a>b)?b:a;}

        // Small helper
        inline void check(unsigned int pos) {  if (pos>=_n) throw EggRuntimeError("tried to access an invalid index"); }

       /* Performs correlation
        *
        * This function works independently from the rest of the class.
        *
        * \param n length of data arrays.
        * \param x first data vector.
        * \param y second data vector.
        * \param r variable to receive the correlation coefficient.
        * \param a variable to receive the regression slope.
        */
        static void _correl(unsigned int n, const int* x, const double* y, double& r, double& a);

        // Distance between pairs
        int* _d;
        
        // D (classical) measure of LD
        double *_D;
        
        // D'
        double *_Dp;
        
        // r, correlation coefficient
        double *_r;
        
        // square r
        double *_r2;
        
        // Data array (not managed by the instance)
        unsigned int *_site1;

        // Data array (not managed by the instance)
        unsigned int *_site2;
        
        // Number of pairs
        unsigned int _n;

     private:
     
        LinkageDisequilibrium(const LinkageDisequilibrium&) { }

        LinkageDisequilibrium& operator=(const LinkageDisequilibrium&) {
            return *this;
        }


        class Interval {
            public:
                Interval(unsigned int, unsigned int);
                unsigned int a() const;
                unsigned int b() const;
                bool good() const;
                void set_false();
            private:
                unsigned int _a;
                unsigned int _b;
                unsigned int _good;
        };


  };
}

#endif