#ifndef __GENOTYPE_H_
#define __GENOTYPE_H_

#include <stdlib.h>
#include <vector>
#include <utility>
#include <map>
#include "function.h"

using namespace std;


class Genotype{
  // record all possible haplotype that is compatible with the observed genotype
 private:
  
  /* the input genotype code, may need to be completed by calling complete_genotype */
  int *gcode;
  
  /* all possible complete genotypes based on gcode */
  vector<int*> *geno_code;
  
  /* all possible decomposition of haplotypes */
  vector<pair<int, int> > hap_vec;

  double weighted_count;
  int count;
  int size;  // number of SNPs in the genotype/haplotype
  int case_num;
  int control_num;

  vector<double> hap_num;  // store Z-\mu, demension a
  

 public:
 
  int info;  // if hap_vec.size()=0 or hap_vec.size()=h(h+1)/2 or hap_vec.size()=total(total+1)/2, info=0, otherwise info=1
  int gname;
  double prob;
  double prob_IBD1;

  /* Constructor */
  Genotype(int *input_code, const map<int, double> & hap_list, int size);
  
  const vector<pair<int, int> > & get_hap_vec(){
    return hap_vec;
  }
  void increase_count(double weight){
    weighted_count+=weight;
  }
  int get_size(){
    return size;
  }
  void set_count(double m){
    weighted_count = m;
  }

  double get_count(){
    return weighted_count;
  }
  int get_num(){
    return count;
  }
  int get_case_num(){
    return case_num;
  }
  int get_control_num(){
    return control_num;
  }
  void increase_num(){
    count++;
  }
  void increase_case_num(){
    case_num++;
  }
  void increase_control_num(){
    control_num++;
  }
  int *get_gcode(){
    return gcode;
  }
  vector<int*> *get_geno_code(){
    return geno_code;
  }
  void record_hap_structure(const map<int, double> & hap_list);  // fill in hap_vec

  // function for IQL_a
  void expected_hap_num(double *freq, map<pair<int, int>, double> & geno_freq_table, map<pair<int, int>, int> & hap_pair_map, int **count, const vector<int> & sites);
  void expectation(map<pair<int, int>, double> & geno_freq_table);
  void expectation_IBD1(map<int, double> & hap_freq_table);
  void compute_gname();
  const vector<double> & get_hap_num(){
    return hap_num;
  }

  void show_content();  // print geno_code and hap_vec
  void show_gstat();  // print gcode with its count
  
  ~Genotype();  // delete


 private:

  // decompose genotype into all possible haplotypes
  vector<int*> * hap_decompose(int *geno, int index, int &bit);

  // if gcode contains missing data, call complete_genotype in constructor
  vector<int*> * complete_genotype(int index=0);
};


void translate_hap(FILE *fp, int code, int size, vector<char> &a0v, vector<char> &a1v);
void translate_hap(FILE *fp, int code, int size);

#endif