models.h

#ifndef _MODELS_H_
#define _MODELS_H_
/*****************************************************************************/
#include <cmath>

#include "lnmultinomial.h"
#include "data_types/allele.h"
#include "typedef.h"
#include "quartet.h"
#include "locus.h"

class models{
private:
lnmultinomial lnMM_, lnMm_, lnmm_;      
                        // The '4' specifies the number of categories of the distribution.
                        // Since these represent the distribution of the four nucleotides 
                        // A, C, G and T, we use 4 categories.
lnmultinomial lnMMP_, lnMmP_, lnmmP_;       

float_t E0_, E1_, E2_;          
public:

    models(void);
    ~models(void);

    float_t loglikelihood(const Locus &, const Allele &);

    float_t goodness_of_fit (Locus &, const Allele &, std::vector <float_t> &, const float_t &);

    void init_gof(const count_t *, const Allele &);

    float_t get_gof(const count_t *, const Allele &);

    float_t genotypelikelihood(const quartet_t &, const Allele &);

    inline float_t lnL(const float_t &logMM, const float_t &logMm, const float &logmm, const count_t *count){
        /*posterior = prior x likelihood */
        E0_=logMM+lnMMP_.lnprob(count);
        E1_=logMm+lnMmP_.lnprob(count);
        E2_=logmm+lnmmP_.lnprob(count);

        if (E0_>E2_) std::swap(E2_, E0_);
        if (E1_>E2_) std::swap(E2_, E1_);
                    
        //We make E2 the largest (i.e. least negative) value of the three. E0 and E1 will often times be 
        //extreamly small, and can even be negative infinity. So we are basically just ensuring that we 
        //don't throw away a lot of precission by returning log(exp(E2) ). 

        return log(exp(E0_-E2_)+exp(E1_-E2_)+1.)+E2_;
    }

    models& operator=(const models& rhs);  
};

float_t *mmModel(Allele *);
float_t *MmModel(Allele *);
float_t *MMModel(Allele *);
float_t *mmModelP(Allele *);
float_t *MmModelP(Allele *);
float_t *MMModelP(Allele *);
#endif