#ifndef __MISSINGTRIO_H_
#define __MISSINGTRIO_H_

#include <vector>
#include <utility>
#include <map>
#include <string.h>
#include "genotypetrio.h"
#include "matrix.h"
#include "missing.h"

using namespace std;


class Missingtrio  // record all possible haplotype that is compatible with the observed genotype
{
 private:

  // enumerate all possible trio_code that have the same missing patten as specified in original
  vector<Genotypetrio *> missing_pat_trio_genotype;

  int count;
  int size;

  int ind_clus;
  int ind_1df;

 public:

  int typed_o;
  int typed_par1;
  int typed_par2;
  
  Missing * M_offspring;
  Missing * M_par1;
  Missing * M_par2;

  double **weight;  // dimension dim * missing_pat_trio_genotype.size(), row: haplotype; col: possible trio genotypes
  double *gfreq;  // dimension missing_pat_trio_genotype.size()
  double **centmatrix;  // variance
  double **derivative;  // F_p
  double **deriv1;
  double **deriv2;
  double **deriv3;
  double **derivative_clus;  // clustered haplotypes
  double **deriv1_clus;
  double **deriv2_clus;
  double **deriv3_clus;
  double **derivative_1df;  // clustered haplotypes
  double **deriv1_1df;
  double **deriv2_1df;
  double **deriv3_1df;
  vector<int> sites;  // record which row and column should be kept in the IQLS
  int dim;
  vector<int> nonzero;  // record the entry of missing_pat_trio_genotype if at lest one weight is nonzero
 
  /* Constructor */
  Missingtrio(Missing * M_offspring, Missing * M_par1, Missing * M_par2, int typed_o, int typed_par1, int typed_par2, const map<int, double> & hap_list, const vector<char *> & content, int size);
  Missingtrio(Missing * M_offspring, int typed_o, int typed_par1, int typed_par2, const map<int, double> & hap_list, const vector<char *> & content, int size);

  void set_vec();
  void set_weight(map<int, double> & hap_freq_table, map<pair<int, int>, double> & geno_freq_table, map<pair<int, int>, int> & hap_pair_map, int **count);
  void compute_centmatrix(map<int, double> & hap_freq_table, double *freq, map<pair<int, int>, double> & geno_freq_table);
  void compute_derivative(map<int, double> & hap_freq_table, map<pair<int, int>, double> & geno_freq_table, map<pair<int, int>, int> & hap_pair_map, int **count, double **rst, double ***hst);
  void compute_deriv1to3(map<int, double> & hap_freq_table, map<pair<int, int>, double> & geno_freq_table, map<pair<int, int>, int> & hap_pair_map, int **count, double **rst, double ***drst);
  void compute_derivative_cluster(map<int, double> & hap_freq_table, double **mul_clus, int col_clus);
  void compute_deriv1to3_cluster(map<int, double> & hap_freq_table, double **mul_clus, int col_clus);
  void compute_derivative_1df(map<int, double> & hap_freq_table, double **mul_1df, int col_1df);
  void compute_deriv1to3_1df(map<int, double> & hap_freq_table, double **mul_1df, int col_1df);
  void keep(const map<int, double> & hap_list, const vector<char *> & content);  // determin which row should be kept
//  vector<int> *indpendent_loc(const vector<int> & missing_site, const map<int, double> & hap_list, const vector<char *> & content, int last);
  void show_gstat();
  void expected_hap_num(map<int, double> & hap_freq_table, double *freq, map<pair<int, int>, double> & geno_freq_table, map<pair<int, int>, int> & hap_pair_map, int **count);
  void cond_par1_one(double **result, map<int, double> & hap_freq_table, map<pair<int, int>, double> & geno_freq_table, map<int, int> & hap_map);
  void cond_par1_two(double **result, map<int, double> & hap_freq_table, map<pair<int, int>, double> & geno_freq_table, map<pair<int, int>, int> & hap_pair_map);
  void cond_par2_one(double **result, map<int, double> & hap_freq_table, map<pair<int, int>, double> & geno_freq_table, map<int, int> & hap_map);
  void cond_par2_two(double **result, map<int, double> & hap_freq_table, map<pair<int, int>, double> & geno_freq_table, map<pair<int, int>, int> & hap_pair_map);
  void cond_off_one(double **result, map<int, double> & hap_freq_table, map<pair<int, int>, double> & geno_freq_table, map<int, int> & hap_map);
  void cond_off_two(double **result, map<int, double> & hap_freq_table, map<pair<int, int>, double> & geno_freq_table, map<pair<int, int>, int> & hap_pair_map);
  void cond_parents_one(double ***result, map<int, double> & hap_freq_table, map<pair<int, int>, double> & geno_freq_table, map<int, int> & hap_map);
  void cond_parents_two(double ***result, map<int, double> & hap_freq_table, map<pair<int, int>, double> & geno_freq_table, map<pair<int, int>, int> & hap_pair_map);
  
  const vector<Genotypetrio *> & get_missing_pat_trio_genotype()
  {
    return missing_pat_trio_genotype;
  }
  int get_num()
  {
    return count;
  }
  void increase_num()
  {
    count++;
  }

  ~Missingtrio();  // delete


 private:

  vector<int*> *complete_missing_haplotype(int *hcode, int index=0);

};

int is_in(vector<int> & s, int m);

#endif