#ifndef __IQLBESTIMATOR_H_
#define __IQLBESTIMATOR_H_

#include <map>
#include "family.h"
#include "missing.h"
#include "function.h"
#include "matrix.h"
#include "missingtrio.h"


class IQLBEstimator
{
 private:
  
  FILE *outfp;
  map<int, double> hap_freq_table;  // record haplotype frequencies
  map<pair<int, int>, double> geno_freq_table;

  int loci_num;
  int sample_size;
  int oneparent;

  vector<Family *> & fam_list;
  vector<Genotypetrio *> & geno_trio_list;  // unique observed trio genotype
  vector<Missingtrio *> & missing_trio_list;  // store the unique trio missing pattern
  map<int, int> geno_missing_trio_map;  // record which trio missing patten does each trio genotype belong

  map<int, int> par1_one;
  map<int, int> par1_two;
  map<int, int> par2_one;
  map<int, int> par2_two;
  map<int, int> off_one;
  map<int, int> off_two;
  map<int, int> parents_one;
  map<int, int> parents_two;

  int **count;  // record how many copies of h_k in pair of haplotypes (i,j), demension (a-1) * a(a+1)/2, entries are 0, 1 or 2
  double ***hst;  // derivative of one parent transition probability
  double test;
  double *Btest;
  double pvalue;
  double *Bpvalue;
  double minpvalue;
  int df;
  int col_1df;
  int test_perform;

  double **work1;
  double **work2;

  int flag_pd;

  map<int, int> hap_map;
  map<pair<int, int>, int> hap_pair_map;
  double **rst;  // record derivative for each pair of haplotypes
  double *freq;
  double *genofreq;
  double ***result_par1_one;
  double ***result_par1_two;
  double ***result_par2_one;
  double ***result_par2_two;
  double ***result_off_one;
  double ***result_off_two;
  double ****result_parents_one;
  double ****result_parents_two;
  
  double **denominator;
  double *numerator;

  int tot_dim;
  double **cholent;
  double **chol;
  double **cholaug;
  double **covMatrix;

 public:
  
  int aff;
  int unaff;
  int unknown;
  map<int, double> aff_hap_freq_table;
  map<int, double> unaff_hap_freq_table;
  map<int, double> unknown_hap_freq_table;
  map<int, double> full_hap_freq_table;
  map<int, double> aff_hap_freq_table_EW;
  map<int, double> unaff_hap_freq_table_EW;
  map<int, double> unknown_hap_freq_table_EW;
  map<int, double> full_hap_freq_table_EW;
  
  IQLBEstimator(vector<Family *> & fam_list, vector<Genotypetrio *> & geno_trio_list, vector<Missingtrio *> & missing_trio_list, map<int, double> & hap_list, int size, int oneparent);
  void init_freq_table();
  void reset_freq_table(map<int, double> & hap_freq);
  void show_hap_freq();
  void show_hap_freq(vector<char> &a0v, vector<char> &a1v);
  void show_geno_freq();
  void update_hap_freq();
  void update_geno_freq();

  const map<int, double> & get_hap_freq(){
    return hap_freq_table;
  }
  int get_loci_num(){
    return loci_num;
  }
  int get_sample_size(){
    return sample_size;
  }
  double get_p_value(){
    return pvalue;
  }
  double get_min_p_value(){
    return minpvalue;
  }
  void show_test(vector<int *> *cluster, vector<char> &a0v, vector<char> &a1v);
  void allocate();
  void NR_Iteration_new();
  void Score_test_new();
  int Score_test_cluster_new(vector<int *> *cluster);
  void deallocate();

  void set_output(FILE *fp)
  {
    outfp = fp;
  }

  void estimate_subgroup();
  void NR_Iteration_subgroup(int group);
  void print_subgroup(vector<char> &a0v, vector<char> &a1v);
  void NR_Iteration_EW();
  void estimate_subgroup_EW();
  void NR_Iteration_subgroup_EW(int group);
  void print_subgroup_EW(vector<char> &a0v, vector<char> &a1v);
  
  ~IQLBEstimator();


 private:

  void compute_count();
  void compute_hst();
  void geno_missing_trio_patten();
  void pair_missing_pattern();

  int relation(int fam, int ind_1, int ind_2, int ind_2f, int ind_2m, double value_o, double IBD_2, double IBD_1, int &diff);
  int relation(int fam, int ind_1, int ind_1f, int ind_1m, int ind_2, int ind_2f, int ind_2m, double value_o, double IBD_2, double IBD_1, int &diff);

};

int compare_equ(int *a, int *b, int loci_num);
    
#endif