#ifndef __IQLAESTIMATOR_H_
#define __IQLAESTIMATOR_H_

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


class IQLAEstimator{

 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;

  vector<Family *> & fam_list;
  vector<Genotype *> & complete_geno_list;
  vector<Genotype *> & geno_list;  // unique observed genotype
  vector<Missing *> & missing_list;  // store the unique missing pattern, unit is one genotype
  map<int, int> geno_missing_map;  // record which missing patten does each genotype belong

  map<int, int> pair_missing_v;
  map<pair<int, int>, int> pair_missing_IBD1;
  map<pair<int, int>, int> pair_missing_IBD2;

  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 test;
  double *Btest;
  double pvalue;
  double *Bpvalue;
  double minpvalue;
  int df;
  int col_1df;
  int test_perform;

  int flag_pd;

  map<pair<int, int>, int> hap_pair_map;
  double **rst;  // record derivative for each pair of haplotypes
  double *freq;
  double **L;
  double ***V;
  double ***mat_cov_1;
  double ***mat_cov_2;

  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;
 
  IQLAEstimator(vector<Family *> & fam_list, vector<Genotype *> & complete_geno_list, vector<Genotype *> & geno_list, vector<Missing *> & missing_list, map<int, double> & hap_list, int size);
  void init_freq_table();
  void reset_freq_table(map<int, double> & hap_freq);
  void show_incom_geno_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);
  
  void single();
  void single_Score_test();

  ~IQLAEstimator();

 private:

  void compute_count();
  void geno_missing_patten();
  void pair_missing_pattern();
};
    
#endif