#ifndef __FUNCTION_H_
#define __FUNCTION_H_

#include <stdlib.h>
#include <string.h>
#include <stdio.h>
#include <math.h>
#include <utility>
#include <map>

using namespace std;

#define        Z_EPSILON      0.000001       /* accuracy of critz approximation */
#define        Z_MAX          6.0            /* maximum meaningful z value */
#define        CHI_EPSILON     0.000001    /* accuracy of critchi approximation */
#define        CHI_MAX     99999.0         /* maximum chi square value */

#define        LOG_SQRT_PI     0.5723649429247000870717135 /* log (sqrt (pi)) */
#define        I_SQRT_PI       0.5641895835477562869480795 /* 1 / sqrt (pi) */
#define        BIGX           20.0         /* max value to represent exp (x) */
#define        ex(x)             (((x) < -BIGX) ? 0.0 : exp (x))


/* Seven functions to calculate expectation of genotypes between two outbred individuals */

double fun1(double coef2, double coef1, double coef0, double p);
double fun2(double coef1, double coef0, double p, double q);
double fun3(double coef0, double p, double q);
double fun4(double coef0, double p, double q1, double q2);
double fun7(double coef2, double coef1, double coef0, double p, double q);
double fun8(double coef1, double coef0, double p, double q1, double q2);
double fun9(double coef0, double q1, double q2, double q3, double q4);

double expfortwo(double IBD_2, double IBD_1, double IBD_0, const pair<int, int> & hp_1, const pair<int, int> & hp_2, map<int, double> & hap_freq_table);
double transProb(const pair<int, int> & hp_0, const pair<int, int> & hp_1, const pair<int, int> & hp_2, map<int, double> & hap_freq_table, map<pair<int, int>, double> & geno_freq_table, int index);
double transProbonep(const pair<int, int> & hp_0, const pair<int, int> & hp_1, map<int, double> & hap_freq_table, map<pair<int, int>, double> & geno_freq_table, int index);
double uncleaunt(const pair<int, int> & hp_0, const pair<int, int> & hp_1, const pair<int, int> & hp_2, map<int, double> & hap_freq_table, map<pair<int, int>, double> & geno_freq_table);

double poz(double z);
double critz(double p);
double pochisq(double x, int df);
double critchi(double p, int df);

double fun1(double p);
double fun2(double p, double q);
double fun7(double p, double q);
double fun8(double p, double q1, double q2);

double expfortwo(const pair<int, int> & hp, map<int, double> & hap_freq_table);
double expfortwo(const pair<int, int> & hp_1, const pair<int, int> & hp_2, map<int, double> & hap_freq_table);

#endif