#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_cdf.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_blas.h>


#include <string>
#include <iostream>
#include <iomanip>

using namespace std;

class GibbsSampler  {

  public:

    GibbsSampler()  {    }
    ~GibbsSampler(void)  
    {
        if ( Tinv != NULL )  
        {
           delete Tinv;
        }
    }
    void printTime(char* message);
    void printMatrix(gsl_matrix *theMatrix, int numRow, int numCol);
    void printMatrix(double *theMatrix, int numRow, int numCol);
    void printMatrixDiag(gsl_matrix *theMatrix, int numRow);
    void runmcmc(gsl_matrix *S, int *P, double threshold, int nindiv,int niter, int nburnin, gsl_vector* outvec, bool notRaoBlack);


    double convertHeritLiabToObserv(double h2L, double threshold);
    void calcPsuedoPheno(int nindiv, gsl_vector* liability, gsl_vector* pseudoPheno);
    void calcPsuedoPhenoObserv(gsl_matrix* observCov,
            int numFamilyMembers, gsl_vector* pmlIndiv, 
            gsl_vector* pseudoPhenoIndiv, bool testing);
    
    double sampleTruncNormal(int pheno, double mu, double sigma, double threshold);
    double expTruncNormal(int pheno, double mu, double sigma, double threshold);


  private:
        double *Tinv;   
        
        void invertT(gsl_matrix *S, int nindiv);
        void invertT2(gsl_matrix *S, int nindiv);
        void invertT3(gsl_matrix *S, int nindiv);
        void invertT4(gsl_matrix *S, int nindiv);

};