#include "peddata.h"


PedData::PedData(){

  MAXTOP = 20;
  snp_size = 0;
  sample_size = 0;
  chrom = 0;
  outfp = stdout;
  sigfp = stdout;
  sigfp_min = stdout;

  beg = -1;
  end = -1;

  outoption = 98;  // a(97) -- choice a; b(98) -- choice b
  full_thresh = 15;
}


int PedData::run(ParamSet &ps){

  fprintf(outfp, "****** Results of the Case-control Association Tests ******\n\n");

  fprintf(stderr, "Loading phenotype data ... \n");
  if(read_pheno(ps.phenofile)<0)
    return -1;

  fprintf(stderr, "Loading IBD coefficients ... \n");
  if(load_ibdcoeff(ps.ibdfile)<0)
    return -1;

  fprintf(stderr, "Loading genotype data ... \n");
  if(read_geno(ps.genofile)<0)
    return -1;

  this->beg = ps.begin;
  this->end = ps.end;
  this->outoption = ps.outoption;
  this->full_thresh = ps.full_thresh;
  
  fprintf(stderr, "Loading parameter ... \n");
  if(load_parameter(ps.parameterfile)<0)
    return -1;

  fprintf(outfp, "The prevalence value used in the MIQLS calculations is %.6f.\n\n", prevalence);
  fprintf(outfp, "The number of markers in each haplotype window is %d.\n\n", nmark);
  
  if(signal_scan()<0)
    return -1;

  return 1;
}


int PedData::signal_scan(){

  fprintf(stderr, "Starting IQLS ... \n");

  int st=0, en=0;  // record the order of start and end SNPs
  if(beg == -1 && end == -1){
    st = 0;
    en = snp_size-1;  // by default
  }
  else{
    if(beg<0 || end<0 || end<beg){
      fprintf(stderr, "The region is misspecified ... \n");
      return -1;
    }
    else{
      if(beg<=marker_info[0].position)
        st = 0;
      if(end>=marker_info[snp_size-1].position)
        en = snp_size-1;
      if(beg>marker_info[0].position){
        for(int i=1;i<snp_size;i++){
          if(marker_info[i-1].position<beg && marker_info[i].position>=beg){
            st = i;
            break;
          }
        }
      }
      if(end<marker_info[snp_size-1].position){
        for(int i=snp_size-2;i>0;i--){
          if(marker_info[i].position<=end && marker_info[i+1].position>end){
            en = i;
            break;
          }
        }
      }

      if(en-st+1<nmark){
        fprintf(stderr, "The region is misspecified, there are only %d markers in the region ... \n", en-st+1);
        return -1;
      }
    }
  }

  assign_gcode();
  
  vector<int> cand_list;
  for(int i=st;i<=en;i++)
    cand_list.push_back(i);

//  vector<int *> * vec = genCombSet(cand_list,nmark,nmark-1);
  vector<int *> * vec = windSet(cand_list,nmark);

  vector<punit> marker_list;
  vector<punit> marker_list_min;
  
  for(int iter=0;iter<vec->size();iter++){
    fprintf(stderr, "markers ");
    for(int i=0;i<nmark;i++)
      fprintf(stderr, "%d ", (*vec)[iter][i]);
    fprintf(stderr, "\n");

    fprintf(outfp, "\n*****************************************************\n\n");
    fprintf(outfp, "Analysis of SNP window ");
    for(int i=0;i<nmark;i++)
      fprintf(outfp, "%s ", marker_info[(*vec)[iter][i]].rsnumber.c_str());
    fprintf(outfp, "\n\n");
    fprintf(outfp, "*****************************************************\n\n");
    
    vector<char> a0v;
    vector<char> a1v;
  
    for(int i=0;i<nmark;i++){
      a0v.push_back(marker_info[(*vec)[iter][i]].allele0);
      a1v.push_back(marker_info[(*vec)[iter][i]].allele1);
    }

    assign_gcode((*vec)[iter]);
    hap_list.clear();

//    fprintf(stderr, "Storing unique genotypes ... \n");
    unique_gcode();
//    assign_info();
    
//    fprintf(stderr, "Get EW estimator ... \n");
    EW_a_estimator();

    map<int, double>::iterator hiter = hap_freq_table.begin();
    while(hiter!=hap_freq_table.end()){
      if(hiter->second>=0.0000001)
        hap_list[hiter->first] = 0;
      hiter++;
    }

    vector<char *> hpcontent;
    hiter = hap_list.begin();
    while(hiter!=hap_list.end()){
      char *pl = new char[nmark+1];
      memset(pl,0,nmark+1);
      int h = hiter->first;
      for(int j=0;j<nmark;j++){
        char bit = h%2;
        pl[j] = bit;
        h = h>>1;
      }
      hpcontent.push_back(pl);
      hiter++;
    }
    
    for(int i=0;i<complete_geno_list.size();i++)
      delete complete_geno_list[i];
    complete_geno_list.clear();

    for(int i=0;i<geno_list.size();i++)
      delete geno_list[i];
    geno_list.clear();

//    fprintf(stderr, "Storing unique genotypes ... \n");
    unique_gcode();
//    assign_info();

//    fprintf(stderr, "Get EW estimator ... \n");
    EW_a_estimator();

    while(1){
      int flag = 1;
      map<int, double>::iterator hiter = hap_freq_table.begin();
      while(hiter!=hap_freq_table.end()){
        if(hiter->second<0.0000001){
          flag = 0;
          break;
        }
        hiter++;
      }
      
      if(flag==1)
        break;
      else{
        hap_list.clear();
        for(int i=0;i<hpcontent.size();i++)
	      delete[] hpcontent[i];
	    hpcontent.clear();
	  
        hiter = hap_freq_table.begin();
        while(hiter!=hap_freq_table.end()){
          if(hiter->second>=0.0000001)
            hap_list[hiter->first] = 0;
          hiter++;
        }

        hiter = hap_list.begin();
        while(hiter!=hap_list.end()){
          char *pl = new char[nmark+1];
	      memset(pl,0,nmark+1);
          int h = hiter->first;
          for(int j=0;j<nmark;j++){
            char bit = h%2;
            pl[j] = bit;
            h = h>>1;
          }
	      hpcontent.push_back(pl);
          hiter++;
        }

        for(int i=0;i<complete_geno_list.size();i++)
          delete complete_geno_list[i];
        complete_geno_list.clear();
  
        for(int i=0;i<geno_list.size();i++)
          delete geno_list[i];
        geno_list.clear();

//        fprintf(stderr, "Storing unique genotypes ... \n");
        unique_gcode();
//        assign_info();

//        fprintf(stderr, "Get EW estimator ... \n");
        EW_a_estimator();   
      }         
    }

//    fprintf(stderr, "Storing unique missing patten ... \n");
    unique_mis_patten(hpcontent);

//    fprintf(stderr, "Get IQL_a estimator ... \n");
    IQL_a_estimator();
    
    while(1){
      int flag = 1;
      map<int, double>::iterator hiter = hap_freq_table.begin();
      while(hiter!=hap_freq_table.end()){
        if(hiter->second<0.00005){
          flag = 0;
          break;
        }
        hiter++;
      }
      
      if(flag==1)
        break;
      else{
        hap_list.clear();
        for(int i=0;i<hpcontent.size();i++)
	      delete[] hpcontent[i];
	    hpcontent.clear();
	  
        hiter = hap_freq_table.begin();
        while(hiter!=hap_freq_table.end()){
          if(hiter->second>=0.00005)
            hap_list[hiter->first] = 0;
          hiter++;
        }

        hiter = hap_list.begin();
        while(hiter!=hap_list.end()){
          char *pl = new char[nmark+1];
	      memset(pl,0,nmark+1);
          int h = hiter->first;
          for(int j=0;j<nmark;j++){
            char bit = h%2;
            pl[j] = bit;
            h = h>>1;
          }
	      hpcontent.push_back(pl);
          hiter++;
        }

        for(int i=0;i<complete_geno_list.size();i++)
          delete complete_geno_list[i];
        complete_geno_list.clear();
  
        for(int i=0;i<geno_list.size();i++)
          delete geno_list[i];
        geno_list.clear();
        
        for(int i=0;i<missing_list.size();i++)
          delete missing_list[i];
        missing_list.clear();

//        fprintf(stderr, "Storing unique genotypes ... \n");
        unique_gcode();
//        assign_info();

//        fprintf(stderr, "Get EW estimator ... \n");
        EW_a_estimator();
        
//        fprintf(stderr, "Storing unique missing patten ... \n");
        unique_mis_patten(hpcontent);
      
//        fprintf(stderr, "Get IQL_a estimator ... \n");
        IQL_a_estimator();   
      }         
    }
        
//    fprintf(stderr, "Get IQL_a estimator ... \n");
    if(outoption==97)
      IQL_a_estimator(a0v,a1v);
    else if(outoption==98){
//      fprintf(stderr, "Storing unique trio genotypes ... \n");
      unique_trio_gcode();

//      fprintf(stderr, "Storing unique trio missing patten ... \n");
      unique_mis_trio_patten(hpcontent);
      
//      fprintf(stderr, "Get IQL_b estimator ... \n");
      IQL_b_estimator();

      while(1){
        int flag = 1;
        map<int, double>::iterator hiter = hap_freq_table.begin();
        while(hiter!=hap_freq_table.end()){
          if(hiter->second<0.00005){
            flag = 0;
            break;
          }
          hiter++;
        }
      
        if(flag==1)
          break;
        else{
          hap_list.clear();
          for(int i=0;i<hpcontent.size();i++)
	        delete[] hpcontent[i];
	      hpcontent.clear();
	  
          hiter = hap_freq_table.begin();
          while(hiter!=hap_freq_table.end()){
            if(hiter->second>=0.00005)
              hap_list[hiter->first] = 0;
            hiter++;
          }

          hiter = hap_list.begin();
          while(hiter!=hap_list.end()){
            char *pl = new char[nmark+1];
	        memset(pl,0,nmark+1);
            int h = hiter->first;
            for(int j=0;j<nmark;j++){
              char bit = h%2;
              pl[j] = bit;
              h = h>>1;
            }
	        hpcontent.push_back(pl);
            hiter++;
          }

          for(int i=0;i<complete_geno_list.size();i++)
            delete complete_geno_list[i];
          complete_geno_list.clear();
  
          for(int i=0;i<geno_list.size();i++)
            delete geno_list[i];
          geno_list.clear();
        
          for(int i=0;i<missing_list.size();i++)
            delete missing_list[i];
          missing_list.clear();
          
          for(int i=0;i<geno_trio_list.size();i++)
            delete geno_trio_list[i];
          geno_trio_list.clear();

          for(int i=0;i<missing_trio_list.size();i++)
            delete missing_trio_list[i];
          missing_trio_list.clear();

//          fprintf(stderr, "Storing unique genotypes ... \n");
          unique_gcode();
//          assign_info();

//          fprintf(stderr, "Get EW estimator ... \n");
          EW_a_estimator();
        
//          fprintf(stderr, "Storing unique missing patten ... \n");
          unique_mis_patten(hpcontent);
      
//          fprintf(stderr, "Storing unique trio genotypes ... \n");
          unique_trio_gcode();

//          fprintf(stderr, "Storing unique trio missing patten ... \n");
          unique_mis_trio_patten(hpcontent);

//          fprintf(stderr, "Get IQL_b estimator ... \n");
          IQL_b_estimator();
        }         
      }      

//      fprintf(stderr, "Get IQL_b estimator ... \n");
      IQL_b_estimator(a0v,a1v);
    }
    
    for(int i=0;i<hpcontent.size();i++)
      delete[] hpcontent[i];
    hpcontent.clear();

    fprintf(outfp, "\n");

//    fprintf(stderr, "Clean heap ... \n\n");
    mem_clean();
    
    punit su;
    su.place = iter;
    su.p = p;
    marker_list.push_back(su);
    
    punit sumin;
    sumin.place = iter;
    sumin.p = p_min;
    marker_list_min.push_back(sumin);              
  }
  
  int top = MAXTOP;
  if(MAXTOP >= marker_list.size())
    top = marker_list.size();

  if(outoption==97)
    fprintf(sigfp, "Top %d results for the full-df MIQLS_a test (list of the SNP haplotype windows having the %d smallest p-values among all SNP haplotype windows that were tested) \n\n", top,top);
  else if(outoption==98)
    fprintf(sigfp, "Top %d results for the full-df MIQLS_b test (list of the SNP haplotype windows having the %d smallest p-values among all SNP haplotype windows that were tested) \n\n", top,top);
  fprintf(sigfp, "**************************************************\n\n");

//  printf("size of p is %d\n", marker_list.size());
  std::sort(marker_list.begin(),marker_list.end(),inc_sort);

  for(int j=0;j<top;j++){
    int iter = marker_list[j].place;
    fprintf(sigfp, "SNP window ");
    for(int i=0;i<nmark;i++)
      fprintf(sigfp, "%s ", marker_info[(*vec)[iter][i]].rsnumber.c_str());
    fprintf(sigfp, "has p-value %g\n", marker_list[j].p);
  }
  
  if(outoption==97)
    fprintf(sigfp_min, "Top %d results for the 1-df MIQLS_a test.  Consider the p-values from all the 1-df tests for SNP haplotype window i, and let the minimum of these be p_i.  Below we report the SNP haplotype windows having the %d smallest values of p_i. \n\n", top,top);
  else if(outoption==98)  
    fprintf(sigfp_min, "Top %d results for the 1-df MIQLS_b test.  Consider the p-values from all the 1-df tests for SNP haplotype window i, and let the minimum of these be p_i.  Below we report the SNP haplotype windows having the %d smallest values of p_i. \n\n", top,top);
  fprintf(sigfp_min, "**************************************************\n\n");

//  printf("size of p is %d\n", marker_list.size());
  std::sort(marker_list_min.begin(),marker_list_min.end(),inc_sort);

  for(int j=0;j<top;j++){
    int iter = marker_list_min[j].place;
    fprintf(sigfp_min, "SNP window ");
    for(int i=0;i<nmark;i++)
      fprintf(sigfp_min, "%s ", marker_info[(*vec)[iter][i]].rsnumber.c_str());
    fprintf(sigfp_min, "has min. 1-df p-value %g\n", marker_list_min[j].p);                 
  }

  return 1;
}


int PedData::read_pheno(char *pheno_file){

  int fam_id, ind_id, father_id, mother_id, sex, phenotype;
  int fam_old = -1;
  int n_aff=0, n_unaff=0, n_unknown=0;

  ifstream phenodata(pheno_file);
  istringstream ins;
  string line;
  int lc=0;  // line number in the file

  vector<int *> fam_member;
  while(getline(phenodata,line)){
    lc++;
    ins.clear();
    ins.str(line);

    ins>>ws;
    if(ins.eof())
      continue;

    if(ins >> fam_id >> ind_id >> father_id >> mother_id >> sex >> phenotype){
      if(phenotype < 0 || phenotype > 2)
        fprintf(stderr, "individual %d from family %d does not have a phenotype value of 0, 1, or 2 ...skipped\n\n", ind_id,fam_id);
      else
        individual_list[ind_id] = -1;

      if(fam_id!=fam_old){
        if(fam_old!=-1){
          fam_list.push_back(new Family(fam_member,fam_old));
          fam_member.clear();
        }

        if(fam_map.find(fam_id)!=fam_map.end()){
          fprintf(stderr, "family %d is not well ordered\n\n", fam_id);
          return 0;
        }
        else{
          fam_map[fam_id] = fam_list.size();
          fam_old=fam_id;
        }
      }

      int *a = new int[5];
      memset(a,0,5*sizeof(int));
      a[0] = ind_id;
      a[1] = father_id;
      a[2] = mother_id;
      a[3] = sex;
      a[4] = phenotype;
      fam_member.push_back(a);
    }
  }

  fam_list.push_back(new Family(fam_member,fam_old));
  fam_member.clear();

  for(int i=0;i<fam_list.size();i++){
    n_aff+=fam_list[i]->Affected;
    n_unaff+=fam_list[i]->UnAffected;
    n_unknown+=fam_list[i]->Unknown;
  }
  sample_size=n_aff+n_unaff+n_unknown;

  fprintf(outfp, "There are %d individuals from %d independent families in the phenotype data file.\n", sample_size,fam_list.size());
  fprintf(outfp, "Of these, %d are affected, %d are unaffected and %d are of unknown phenotype.\n\n", n_aff,n_unaff,n_unknown);

  return 1;
}


int PedData::load_ibdcoeff(char *ibdco_file){

  int ind_1, ind_2;
  double coef8, coef7, coef6, coef5, coef4, coef3, coef2, coef1, coef0;
  int loc;

  map<int, int> ind_fam;
  for(int i=0;i<fam_list.size();i++)
    for(int j=0;j<fam_list[i]->N;j++)
      ind_fam[fam_list[i]->fam_member[j][0]] = i;

  ifstream ibdcodata(ibdco_file);
  istringstream ins;
  string line;
  int lc=0;  // line number in the file

  while(getline(ibdcodata,line)){
    lc++;
    ins.clear();
    ins.str(line);

    ins>>ws;
    if(ins.eof())
      continue;

    if(ins >> ind_1 >> ind_2 >> coef8 >> coef7 >> coef6>> coef5 >> coef4 >> coef3 >> coef2 >> coef1 >> coef0){
      if(ind_1==ind_2){
        if(ind_fam.find(ind_1)==ind_fam.end())
          fprintf(stderr, "Skip line %d in ibd coefficients file \"%s\"\n", lc,ibdco_file);
        else{
          loc = ind_fam[ind_1];
          if(coef8>0 || coef7>0 || coef6>0 || coef5>0 || coef4>0 || coef3>0 || coef1>0 || coef0>0 || fabs(coef2-1.0)>1e-8){
            fprintf(stderr, "Skip line %d in ibd coefficients file \"%s\"\n", lc,ibdco_file);
            if(individual_list.find(ind_1)!=individual_list.end())
              individual_list.erase(ind_1);
          }
          else
            fam_list[loc]->kinship_coeff[pair<int, int>(ind_1,ind_2)] = 0;
        }
      }
      else{
        if(ind_fam.find(ind_1)==ind_fam.end() || ind_fam.find(ind_2)==ind_fam.end())
          fprintf(stderr, "Skip line %d in ibd coefficients file \"%s\"\n", lc,ibdco_file);
        else{
          loc = ind_fam[ind_1];
          if(ind_fam[ind_2]!=loc)
            fprintf(stderr, "Skip line %d in ibd coefficients file \"%s\"\n", lc,ibdco_file);
          else{
            if(coef8>0 || coef7>0 || coef6>0 || coef5>0 || coef4>0 || coef3>0)
              fprintf(stderr, "Skip line %d in ibd coefficients file \"%s\"\n", lc,ibdco_file);
            else{
              fam_list[loc]->IBD2_coeff[pair<int, int>(ind_1,ind_2)] = coef2;
              fam_list[loc]->IBD1_coeff[pair<int, int>(ind_1,ind_2)] = coef1;
              fam_list[loc]->IBD0_coeff[pair<int, int>(ind_1,ind_2)] = coef0;
              fam_list[loc]->kinship_coeff[pair<int, int>(ind_1,ind_2)] = coef2/2.0+coef1/4.0;
              fam_list[loc]->IBD2_coeff[pair<int, int>(ind_2,ind_1)] = coef2;
              fam_list[loc]->IBD1_coeff[pair<int, int>(ind_2,ind_1)] = coef1;
              fam_list[loc]->IBD0_coeff[pair<int, int>(ind_2,ind_1)] = coef0;
              fam_list[loc]->kinship_coeff[pair<int, int>(ind_2,ind_1)] = coef2/2.0+coef1/4.0;
            }
          }
        }
      }
    }
  }

  ind_fam.clear();

  return 1;
}


int PedData::read_geno(char *geno_file){

  string marker, alias;
  int chr, pos;
  char orientation, allele0, allele1;
  string geno;
  vector<string> genotype_line;

  int ind_num, count=0;
  vector<int> record;

  ifstream genodata(geno_file);
  istringstream ins;
  string line;
  int lc=0;  // line number in the file

  if(getline(genodata,line)){
    // get the first line
    ins.clear();
    ins.str(line);

    ins>>ws;
    if(ins >> alias >> alias >> alias >> alias >> alias >> alias){  // rsnumber; chromosome; position; orientation; alleleA; alleleB
      while(ins >> ind_num){
        if(individual_list.find(ind_num)!=individual_list.end()){
          individual_list[ind_num] = count;
          count++;
          record.push_back(1);
        }
        else  // individual is not in phenotype file
          record.push_back(0);
      }

      people_size = count;
      fprintf(outfp, "There are %d individuals listed in the genotype file, %d of whom are also in the phenotype file (excluding inbred individuals).\n\n", record.size(),people_size);
    }
  }

  while(getline(genodata,line)){
    lc++;
    ins.clear();
    ins.str(line);

    ins>>ws;
    if(ins.eof())
      continue;

    if(ins >> marker >> chr >> pos >> orientation >> allele0 >> allele1){
      if(chrom == 0)
        chrom = chr;

      MarkerInfo mi;    
      mi.rsnumber = marker;
      mi.position = pos;

      if(allele0>=97)  // convert to capital letter
	    allele0-=32;
      if(allele1>=97)
	    allele1-=32;
      
      // direct assignment, no assumption about allele0 < allele1
      mi.allele0 = allele0;
      mi.allele1 = allele1;

      marker_info.push_back(mi);
      snp_size++;

      count = 0;
      while(ins >> geno){
        if(record[count]==1)
          genotype_line.push_back(geno);
        count++;
      }

      if(count!=record.size())
        printf("error in genofile\n");

      if(people_size != genotype_line.size()){
        fprintf(stderr, "Error: Unequal number of persons recorded in genotyo data file\n");
        return -1;
      }

      string *cl = new string[genotype_line.size()];
      for(int i=0;i<genotype_line.size();i++)
        cl[i] = genotype_line[i];
      genotype_line.clear();
      raw_content.push_back(cl);
    }
    else
      continue;
  }

  record.clear();

  return convert_coding();
}


int PedData::convert_coding(){

  for(int i=0;i<snp_size;i++){
    char ta0 = marker_info[i].allele0;
    char ta1 = marker_info[i].allele1;

    map<string, int> geno_code;
    geno_code.clear();
    geno_code["NN"] = -1;

    if(ta0=='A' && ta1=='T'){
      geno_code["AA"] = 0;
      geno_code["AT"] = geno_code["TA"] = 1;
      geno_code["TT"] = 2;
    }
	else if(ta0=='A' && ta1=='C'){
      geno_code["AA"] = 0;
      geno_code["AC"] = geno_code["CA"] = 1;
      geno_code["CC"] = 2;
    }
    else if(ta0=='A' && ta1=='G'){
      geno_code["AA"] = 0;
      geno_code["AG"] = geno_code["GA"] = 1;
      geno_code["GG"] = 2;
    }
    else if(ta0=='T' && ta1=='A'){
      geno_code["TT"] = 0;
      geno_code["TA"] = geno_code["AT"] = 1;
      geno_code["AA"] = 2;
    }
    else if(ta0=='T' && ta1=='C'){
      geno_code["TT"] = 0;
      geno_code["TC"] = geno_code["CT"] = 1;
      geno_code["CC"] = 2;
    }
    else if(ta0=='T' && ta1=='G'){
      geno_code["TT"] = 0;
      geno_code["TG"] = geno_code["GT"] = 1;
      geno_code["GG"] = 2;
    }
    else if(ta0=='C' && ta1=='A'){
      geno_code["CC"] = 0;
      geno_code["CA"] = geno_code["AC"] = 1;
      geno_code["AA"] = 2;
    }
    else if(ta0=='C' && ta1=='T'){
      geno_code["CC"] = 0;
      geno_code["CT"] = geno_code["TC"] = 1;
      geno_code["TT"] = 2;
    }
    else if(ta0=='C' && ta1=='G'){
      geno_code["CC"] = 0;
      geno_code["CG"] = geno_code["GC"] = 1;
      geno_code["GG"] = 2;
    }
    else if(ta0=='G' && ta1=='A'){
      geno_code["GG"] = 0;
      geno_code["GA"] = geno_code["AG"] = 1;
      geno_code["AA"] = 2;
    }
    else if(ta0=='G' && ta1=='T'){
      geno_code["GG"] = 0;
      geno_code["GT"] = geno_code["TG"] = 1;
      geno_code["TT"] = 2;
    }
    else if(ta0=='G' && ta1=='C'){
      geno_code["GG"] = 0;
      geno_code["GC"] = geno_code["CG"] = 1;
      geno_code["CC"] = 2;
    }

    int *data = new int[people_size];
    for(int j=0;j<people_size;j++)
      data[j] = geno_code[raw_content[i][j]];
    delete[] raw_content[i];
    content.push_back(data);
  }
  
  raw_content.clear();

  return 1;
}


int PedData::load_parameter(char *parameter_file){

  int loci;
  double kp;

  ifstream parameterdata(parameter_file);
  istringstream ins;
  string line;

  getline(parameterdata,line);
  ins.clear();
  ins.str(line);

  ins>>ws;
  if(ins >> kp)
    prevalence = kp;

  getline(parameterdata,line);
  ins.clear();
  ins.str(line);

  ins>>ws;
  if(ins >> loci)
    nmark = loci;

  return 1;
}


void PedData::assign_gcode(){

  for(int i=0;i<fam_list.size();i++)
    fam_list[i]->set_gcode(nmark);
}


void PedData::assign_gcode(int *loci){

  int *gd = new int[nmark];
  int count=0;
  int count_aff=0;
  int count_unaff=0;
  int count_unknown=0;

  for(int i=0;i<fam_list.size();i++){
    fam_list[i]->set_gcode();
    for(int j=0;j<fam_list[i]->N;j++){
      int ind_num = fam_list[i]->fam_member[j][0];
      if(individual_list.find(ind_num)!=individual_list.end()){
        int col = individual_list[ind_num];
        if(col>=0){
          for(int k=0;k<nmark;k++)
            gd[k] = content[loci[k]][col];
          fam_list[i]->set_gcode(gd,j);

          for(int k=0;k<nmark;k++){
            if(gd[k]!=-1){
              count++;
              if(fam_list[i]->fam_member[j][4]==2)
                count_aff++;
              else if(fam_list[i]->fam_member[j][4]==1)
                count_unaff++;
              else if(fam_list[i]->fam_member[j][4]==0)
                count_unknown++;

              break;
            }
          }
        }
      }
    }
    fam_list[i]->fill_typed();
  }
  delete[] gd;

  fprintf(outfp, "There are %d individuals with genotype data for at least one SNP in the current haplotype window.\n", count);
  fprintf(outfp, "Of these, %d are affected, %d are unaffected and %d are of unknown phenotype. \n\n", count_aff,count_unaff,count_unknown);
}


void PedData::assign_info(){

  for(int i=0;i<fam_list.size();i++)
    fam_list[i]->fill_info_gcode(geno_list);
}


void PedData::unique_gcode(){

  complete_geno_list.clear();
  int count=0;
  int *gd = new int[nmark];
  for(int k=0;k<nmark;k++)
    gd[k] = -1;  // genotypes are missing at all markers
  Genotype *geno = new Genotype(gd,hap_list,nmark);
  vector<int*> *complete_geno_code = geno->get_geno_code();  // complete-data genotypes
  count = complete_geno_code->size();
//  printf("There are %d complete genotypes\n", count);

  for(int i=0;i<count;i++){
    Genotype *complete_geno = new Genotype((*complete_geno_code)[i],hap_list,nmark);
    const vector<pair<int, int> > hv = complete_geno->get_hap_vec();
    if(hv.size()>0){
      complete_geno->compute_gname();
      complete_geno_list.push_back(complete_geno);
    }
    else
      delete complete_geno;
  }
  delete geno;
//  printf("There are %d complete genotypes stored\n", complete_geno_list.size());

  geno_list.clear();
  map<int, int> gmap;

  for(int i=0;i<fam_list.size();i++)
    for(int j=0;j<fam_list[i]->N;j++){
      int code=0;
      for(int k=0;k<nmark;k++){
        int gc = fam_list[i]->gcode[j][k];
        gd[k] = gc;
        code = 10*code+gc;
      }

      if(gmap.find(code)==gmap.end()){
        gmap[code] = geno_list.size();
        geno_list.push_back(new Genotype(gd,hap_list,nmark));
      }

      fam_list[i]->set_cat(gmap[code],j);
      geno_list[gmap[code]]->increase_num();
    }
  delete[] gd;
/*
  for(int i=0;i<complete_geno_list.size();i++)
    complete_geno_list[i]->show_gstat();
  printf("\n");

  for(int i=0;i<geno_list.size();i++)
    geno_list[i]->show_gstat();
*/
}


void PedData::EW_a_estimator(){

  HFEstimator est(fam_list,geno_list,hap_list,nmark);
  
  int count = 1;
  int flag;
  while(1){
    map<int, double> of = est.get_hap_freq();
    est.EM_Iteration(0);
    map<int, double> nf = est.get_hap_freq();

    double diff = 0;
    flag = 1;
    map<int, double>::iterator iter=nf.begin();
    while(iter!=nf.end()){
      diff += (nf[iter->first]-of[iter->first])*(nf[iter->first]-of[iter->first]);
      if(nf[iter->first]<0){
        flag = 0;
        break;                    
      }
      iter++;                                
    }

    if(diff<1e-8 || count>100 || flag==0)
      break;

    count++;
  }

  hap_freq_table = est.get_hap_freq();

//  printf("It takes %d iterations.\n", count);
//  est.show_hap_freq();
/*
  printf("EW estimator:\n");
  map<int, double>::iterator iter = hap_freq_table.begin();
  while(iter!=hap_freq_table.end()){
    printf("%d \t %.10f\n", iter->first,iter->second);
    iter++;
  }
  printf("\n");
*/
}


void PedData::CDW_a_estimator(){

  // compute complete-data weight
  for(int i=0;i<fam_list.size();i++)
    fam_list[i]->compute_weight();

  // sum up
  for(int i=0;i<fam_list.size();i++)
    for(int j=0;j<fam_list[i]->N;j++){
      int cat=fam_list[i]->cat[j];
      geno_list[cat]->increase_count(fam_list[i]->weight[j]);
    }

  HFEstimator est(fam_list,geno_list,hap_list,nmark);
  est.reset_freq_table(hap_freq_table);

  int count = 1;
  int flag;
  while(1){
    map<int, double> of = est.get_hap_freq();
    est.EM_Iteration(1);
    map<int, double> nf = est.get_hap_freq();

    double diff = 0;
    flag = 1;
    map<int, double>::iterator iter=nf.begin();
    while(iter!=nf.end()){
      diff += (nf[iter->first]-of[iter->first])*(nf[iter->first]-of[iter->first]);
      if(nf[iter->first]<0){
        flag = 0;
        break;                    
      }
      iter++;                                
    }

    if(diff<1e-6 || count>20 || flag==0)
      break;

    count++;
  }

  if(flag==1)
    hap_freq_table = est.get_hap_freq();

//  printf("It takes %d iterations.\n", count);
//  est.show_hap_freq();
}


void PedData::unique_mis_patten(const vector<char *> & hpcontent){

  missing_list.clear();
  map<int, int> mis_map;

  for(int i=0;i<fam_list.size();i++)
    for(int j=0;j<fam_list[i]->N;j++){
      int code=0;
      vector<int> missing_site;
      for(int k=0;k<nmark;k++){
        int gc = fam_list[i]->gcode[j][k];
        if(gc==-1){
          code = 10*code+gc;
          missing_site.push_back(k);
        }
        else
          code*=10;
      }

      if(mis_map.find(code)==mis_map.end()){
        mis_map[code] = missing_list.size();
        missing_list.push_back(new Missing(missing_site,hap_list,hpcontent,nmark));
      }

      fam_list[i]->set_mis_cat(mis_map[code],j);
      missing_list[mis_map[code]]->increase_num();
      missing_site.clear();
    }
/*
  for(int i=0;i<missing_list.size();i++)
    printf("%d %d, ", missing_list[i]->dim,missing_list[i]->mapping_function.size());
  printf("\n");
*/
}


void PedData::IQL_a_estimator(){

  for(int i=0;i<complete_geno_list.size();i++)
    complete_geno_list[i]->compute_gname();

  // IQLa
  IQLAEstimator est(fam_list,complete_geno_list,geno_list,missing_list,hap_list,nmark);
  est.set_output(outfp);
  est.reset_freq_table(hap_freq_table);

  int count = 1;
  int flag;
  est.allocate();
  while(1){
    map<int, double> of = est.get_hap_freq();
    est.NR_Iteration_new();
    map<int, double> nf = est.get_hap_freq();

    double diff = 0;
    flag = 1;
    map<int, double>::iterator iter=nf.begin();
    while(iter!=nf.end()){
      diff += (nf[iter->first]-of[iter->first])*(nf[iter->first]-of[iter->first]);
      if(nf[iter->first]<0.0000001){
        flag = 0;
        break;                    
      }
      iter++;                                
    }

    if(diff<1e-5 || count>7 || flag==0)
      break;

    count++;
  }

  hap_freq_table = est.get_hap_freq();
  est.deallocate();
}


void PedData::IQL_a_estimator(vector<char> &a0v, vector<char> &a1v){
  // IQLa
  IQLAEstimator est(fam_list,complete_geno_list,geno_list,missing_list,hap_list,nmark);
  est.set_output(outfp);
  est.reset_freq_table(hap_freq_table);

  est.allocate();
  // need frequencies in affected, unaffected, and unknown group
  est.estimate_subgroup();
  est.estimate_subgroup_EW();

  freq_thresh = double(full_thresh)/double(est.aff+est.unaff+est.unknown)/2.0;
  
  // compute test weight
  double mean = -prevalence/(1-prevalence);
  for(int i=0;i<fam_list.size();i++)
    fam_list[i]->compute_test_weight(mean);

  int row = hap_freq_table.size();
  hap_clus_list.clear();
  int count = 0;
  map<int, double>::iterator iter=hap_freq_table.begin();
  while(iter!=hap_freq_table.end()){
    hap_unit hu;
    hu.name = iter->first;
    hu.order = count;
    hu.freq = iter->second;
    hap_clus_list.push_back(hu);
    count++;
    iter++;
  }

  std::sort(hap_clus_list.begin(),hap_clus_list.end(),dec_sort);  // sort based on haplotype frequencies
  if(row!=hap_clus_list.size()){
    printf("error, missing haplotypes\n");
    exit(1);
  }

  cluster = new vector<int *>();                                        
  hap_cluster();
  if(est.Score_test_cluster_new(cluster)>0){
    est.show_test(cluster,a0v,a1v);
    p = est.get_p_value();
    p_min = est.get_min_p_value();
  }
  else{
    p = 1;
    p_min = 1;   
  }
  est.deallocate();
  cluster->clear();
  delete cluster;
/*
  if(hap_clus_list[row-1].freq<freq_thresh){
    cluster = new vector<int *>();                                        
    hap_cluster();
    est.Score_test_cluster_new(cluster);
    cluster->clear();
    delete cluster;
  }
  else{
    // association test
    // IQLS_a test
    est.Score_test_new();
  }  
*/
  est.print_subgroup(a0v,a1v);
  est.print_subgroup_EW(a0v,a1v);
}


void PedData::unique_trio_gcode(){

  geno_trio_list.clear();
  vector<trio> triolist;
  oneparent = 0;

  for(int i=0;i<fam_list.size();i++){
    for(int j=0;j<fam_list[i]->N;j++){
      trio m;
      int id_1, id_2, ind_1, ind_2;
      int loc;
      int founder=0;
      int cat_1, cat_2;
      int typed_1, typed_2;
      int swi=0;

      m.cat_o = fam_list[i]->cat[j];
      m.typed_o = fam_list[i]->typed[j];
      
      ind_1 = fam_list[i]->fam_member[j][1];
      ind_2 = fam_list[i]->fam_member[j][2];
      if(ind_1==0 && ind_2==0){  // j is a founder
        m.cat_par1 = -1;
        m.typed_par1 = 0;
        m.cat_par2 = -1;
        m.typed_par2 = 0;
        founder = 1;
      }
      else{
        id_1 = fam_list[i]->member_map[ind_1];
        cat_1 = fam_list[i]->cat[id_1];
        typed_1 = fam_list[i]->typed[id_1];
        id_2 = fam_list[i]->member_map[ind_2];
        cat_2 = fam_list[i]->cat[id_2];
        typed_2 = fam_list[i]->typed[id_2];
        
        if(typed_1==typed_2){  // 00 or 11
          if(cat_1<=cat_2){
            m.cat_par1 = cat_1;
            m.typed_par1 = typed_1;
            m.cat_par2 = cat_2;
            m.typed_par2 = typed_2;               
          }
          else{
            m.cat_par1 = cat_2;
            m.typed_par1 = typed_2;
            m.cat_par2 = cat_1;
            m.typed_par2 = typed_1;
            swi = 1;   
          }
        }
        else if(typed_1<typed_2){  // 01
          m.cat_par1 = cat_2;
          m.typed_par1 = typed_2;
          m.cat_par2 = cat_1;
          m.typed_par2 = typed_1;
          swi = 1;
        }
        else{  // 10
          m.cat_par1 = cat_1;
          m.typed_par1 = typed_1;
          m.cat_par2 = cat_2;
          m.typed_par2 = typed_2;     
        }                                 
      }       

      loc = find_loc(triolist,m);
      if(loc==-1){  // not found
        fam_list[i]->set_trio_cat(geno_trio_list.size(),j);
        triolist.push_back(m);
        if(founder==1)
          geno_trio_list.push_back(new Genotypetrio(geno_list[m.cat_o],m.typed_o,0,0,nmark));
        else
          geno_trio_list.push_back(new Genotypetrio(geno_list[m.cat_o],geno_list[m.cat_par1],geno_list[m.cat_par2],m.typed_o,m.typed_par1,m.typed_par2,nmark));

        geno_trio_list[geno_trio_list.size()-1]->increase_num();

        if(oneparent==0){
          if((typed_1==1 && typed_2==0) || (typed_1==0 && typed_2==1))
            oneparent = 1;
        }
      }
      else{
        fam_list[i]->set_trio_cat(loc,j);
        geno_trio_list[loc]->increase_num();
      }
    }
  }

  triolist.clear();
}


void PedData::unique_mis_trio_patten(const vector<char *> & hpcontent){

  missing_trio_list.clear();
  vector<trio> mis_triolist;

  for(int i=0;i<fam_list.size();i++){
    for(int j=0;j<fam_list[i]->N;j++){
      trio m;
      int id_1, id_2, ind_1, ind_2;
      int loc;
      int founder=0;
      int cat_1, cat_2;
      int typed_1, typed_2;
      int swi=0;

      m.cat_o = fam_list[i]->mis_cat[j];
      m.typed_o = fam_list[i]->typed[j];

      ind_1 = fam_list[i]->fam_member[j][1];
      ind_2 = fam_list[i]->fam_member[j][2];
      if(ind_1==0 && ind_2==0){  // j is a founder
        m.cat_par1 = -1;
        m.typed_par1 = 0;
        m.cat_par2 = -1;
        m.typed_par2 = 0;
        founder = 1;
      }
      else{
        id_1 = fam_list[i]->member_map[ind_1];
        cat_1 = fam_list[i]->mis_cat[id_1];
        typed_1 = fam_list[i]->typed[id_1];
        id_2 = fam_list[i]->member_map[ind_2];
        cat_2 = fam_list[i]->mis_cat[id_2];
        typed_2 = fam_list[i]->typed[id_2];
        
        if(typed_1==typed_2){  // 00 or 11
          if(cat_1<=cat_2){
            m.cat_par1 = cat_1;
            m.typed_par1 = typed_1;
            m.cat_par2 = cat_2;
            m.typed_par2 = typed_2;               
          }
          else{
            m.cat_par1 = cat_2;
            m.typed_par1 = typed_2;
            m.cat_par2 = cat_1;
            m.typed_par2 = typed_1;
            swi = 1;   
          }
        }
        else if(typed_1<typed_2){  // 01
          m.cat_par1 = cat_2;
          m.typed_par1 = typed_2;
          m.cat_par2 = cat_1;
          m.typed_par2 = typed_1;
          swi = 1;
        }
        else{  // 10
          m.cat_par1 = cat_1;
          m.typed_par1 = typed_1;
          m.cat_par2 = cat_2;
          m.typed_par2 = typed_2;     
        }                                 
      }
      
      loc = find_loc(mis_triolist,m);
      if(loc==-1){  // not found  
        fam_list[i]->set_mis_trio_cat(missing_trio_list.size(),j);
        mis_triolist.push_back(m);
        if(founder==1)
          missing_trio_list.push_back(new Missingtrio(missing_list[m.cat_o],m.typed_o,0,0,hap_list,hpcontent,nmark));
        else
          missing_trio_list.push_back(new Missingtrio(missing_list[m.cat_o],missing_list[m.cat_par1],missing_list[m.cat_par2],m.typed_o,m.typed_par1,m.typed_par2,hap_list,hpcontent,nmark));

        missing_trio_list[missing_trio_list.size()-1]->increase_num();
      }
      else{
        fam_list[i]->set_mis_trio_cat(loc,j);
        missing_trio_list[loc]->increase_num();
      }
      fam_list[i]->set_mis_par_order(swi,j);
    }
  }

  mis_triolist.clear();
}


void PedData::IQL_b_estimator(){
  
  IQLBEstimator est(fam_list,geno_trio_list,missing_trio_list,hap_list,nmark,oneparent);
  est.set_output(outfp);
  est.reset_freq_table(hap_freq_table);

  int count = 1;
  int flag;
  est.allocate();
  while(1){
    map<int, double> of = est.get_hap_freq();
    est.NR_Iteration_new();
    map<int, double> nf = est.get_hap_freq();

    double diff = 0;
    flag = 1;
    map<int, double>::iterator iter=nf.begin();
    while(iter!=nf.end()){
      diff += (nf[iter->first]-of[iter->first])*(nf[iter->first]-of[iter->first]);
      if(nf[iter->first]<0.0000001){
        flag = 0;
        break;                    
      }
      iter++;                                
    }

    if(diff<1e-5 || count>5 || flag==0)
      break;

    count++;
  }
  
  hap_freq_table = est.get_hap_freq();
  est.deallocate();
}


void PedData::IQL_b_estimator(vector<char> &a0v, vector<char> &a1v){
  
  IQLBEstimator est(fam_list,geno_trio_list,missing_trio_list,hap_list,nmark,oneparent);
  est.set_output(outfp);
  est.reset_freq_table(hap_freq_table);

  est.allocate();
  // need frequencies in affected, unaffected, and unknown group
  est.estimate_subgroup();
  est.estimate_subgroup_EW();

  freq_thresh = double(full_thresh)/double(est.aff+est.unaff+est.unknown)/2.0;

  // compute test weight
  double mean = -prevalence/(1-prevalence);
  for(int i=0;i<fam_list.size();i++){
    fam_list[i]->compute_test_weight(mean);
    fam_list[i]->compute_test_trio_weight(mean);
  }
/*
  // compute test trio weight
  double mean = -prevalence/(1-prevalence);
  for(int i=0;i<fam_list.size();i++)
    fam_list[i]->compute_test_trio_weight(mean);
*/
  int row = hap_freq_table.size();
  hap_clus_list.clear();
  int count = 0;
  map<int, double>::iterator iter=hap_freq_table.begin();
  while(iter!=hap_freq_table.end()){
    hap_unit hu;
    hu.name = iter->first;
    hu.order = count;
    hu.freq = iter->second;
    hap_clus_list.push_back(hu);
    count++;
    iter++;
  }

  std::sort(hap_clus_list.begin(),hap_clus_list.end(),dec_sort);  // sort based on haplotype frequencies
  if(row!=hap_clus_list.size()){
    printf("error, missing haplotypes\n");
    exit(1);
  }

  cluster = new vector<int *>();                                        
  hap_cluster();
  if(est.Score_test_cluster_new(cluster)>0){
    est.show_test(cluster,a0v,a1v);
    p = est.get_p_value();
    p_min = est.get_min_p_value();
  }
  else{
    p = 1;
    p_min = 1;   
  }
  est.deallocate();
  cluster->clear();
  delete cluster;  
/*
  if(hap_clus_list[row-1].freq<freq_thresh){
    cluster = new vector<int *>();                                        
    hap_cluster();
    est.Score_test_cluster_new(cluster);
    cluster->clear();
    delete cluster;
  }
  else{
    // association test
    // IQLS_b test
    est.Score_test_new();
  }
*/
/*
  if(hap_clus_list[row-1].freq>=freq_thresh){
    if(est.aff>0){
      iter = hap_freq_table.begin();
      double min_hap = est.aff_hap_freq_table[iter->first];
      iter++;
      while(iter!=hap_freq_table.end()){
        if(est.aff_hap_freq_table[iter->first]<min_hap)
          min_hap = est.aff_hap_freq_table[iter->first];
        iter++;
      }
      if(2*est.aff*min_hap<5)
        fprintf(outfp, "The full-df p-value might not be exact because of the small number of expected haplotype counts in affecteds\n");
    }    
    if(est.unaff>=est.unknown){
      if(est.unaff>0){
        iter = hap_freq_table.begin();
        double min_hap = est.unaff_hap_freq_table[iter->first];
        iter++;
        while(iter!=hap_freq_table.end()){
          if(est.unaff_hap_freq_table[iter->first]<min_hap)
            min_hap = est.unaff_hap_freq_table[iter->first];
          iter++;
        }
        if(2*est.unaff*min_hap<5)
          fprintf(outfp, "The full-df p-value might not be exact because of the small number of expected haplotype counts in unaffecteds\n");
      }
    }
    else{
      if(est.unknown>0){
        iter = hap_freq_table.begin();
        double min_hap = est.unknown_hap_freq_table[iter->first];
        iter++;
        while(iter!=hap_freq_table.end()){
          if(est.unknown_hap_freq_table[iter->first]<min_hap)
            min_hap = est.unknown_hap_freq_table[iter->first];
          iter++;
        }
        if(2*est.unknown*min_hap<5)
          fprintf(outfp, "The full-df p-value might not be exact because of the small number of expected haplotype counts in unknowns\n");
      }
    }
    fprintf(outfp, "\n");   
  }     
*/
  est.print_subgroup(a0v,a1v);
  est.print_subgroup_EW(a0v,a1v);
}


void PedData::hap_cluster(){
     
  int row = hap_freq_table.size();
  
  // determin the common haplotypes
  double max_freq = hap_clus_list[0].freq;
  double cut_off = freq_thresh;
  if(cut_off>max_freq)
    cut_off = max_freq;
    
  int count = 1;
  for(int i=1;i<row;i++){
    if(hap_clus_list[i].freq>=cut_off)
      count++;
    else
      break;
  }
//  printf("There are %d common haplotypes\n", count);
//  fprintf(outfp, "%d\t", count);

  for(int i=0;i<count;i++){
    int *a = new int[row-count+2];
    memset(a,-1,(row-count+2)*sizeof(int));
    a[0] = 1;  // record how many haplotypes in that cluster
    a[1] = hap_clus_list[i].order;
    cluster->push_back(a);
  }
  int *a = new int[row-count+1];
  memset(a,-1,(row-count+1)*sizeof(int));
  a[0] = 0;  // record how many haplotypes in the last cluster
  cluster->push_back(a);

  // cluster rare haplotypes
  for(int i=count;i<row;i++){
    int sign = 0;
    for(int j=0;j<count;j++){
      //if i and j differ at only one locus, combine i and j
      int mutation = 0;
      for(int k=0;k<nmark;k++){
        int mask = 1<<k;
        int s = (hap_clus_list[i].name&mask)>>k;
        int n = (hap_clus_list[j].name&mask)>>k;
        if(s!=n)
          mutation++;
      }

      if(mutation==0)
        printf("error, the two haplotypes are the same\n");
      else if(mutation==1){
        // sign i to j
        int loc = (*cluster)[j][0];
        (*cluster)[j][loc+1] = hap_clus_list[i].order;
        (*cluster)[j][0]++;
        sign = 1;
        break;
      }
    }

    if(sign==0){
      // sign i to the last cluster
      int loc = (*cluster)[count][0];
      (*cluster)[count][loc+1] = hap_clus_list[i].order;
      (*cluster)[count][0]++;
    }
  }
/*
  for(int i=0;i<=count;i++){
    printf("cluster %d: ", i+1);
    for(int j=0;j<=(*cluster)[i][0];j++)
      printf("%d, ", (*cluster)[i][j]);
    printf("\n");
  }
  printf("\n");
*/  
}


void PedData::print_data()
{
  for(int i=0;i<fam_list.size();i++)
    fam_list[i]->print();
}


vector<int *> *PedData::genCombSet(vector<int> & all, const int tsize, int index)
{
  // terminating condition 
  if(index==0)
  {
    vector<int *> *vec = new vector<int*>();
    for(int i=0;i<all.size();i++)
    {
      int *a1 = new int[tsize];
      memset(a1,0,tsize*sizeof(int));
      a1[0] = all[i];
      
      vec->push_back(a1);
    }
    
    //printf("returned\n");
    return vec;
  }

  // else recursive run
  vector<int *> *nvec = new vector<int*>();
  for(int i=0;i<all.size();i++)
  {
    vector<int> nall = all;
    nall.erase(nall.begin(),nall.begin()+i+1);
    vector<int *> *vec = genCombSet(nall,tsize,index-1);
    for(int j=0;j<vec->size();j++)
    {
      (*vec)[j][index] = all[i];
      nvec->push_back((*vec)[j]);
    }
    vec->clear();
    delete vec;
  }
  
  return nvec;
}


vector<int *> *PedData::windSet(vector<int> & all, const int tsize)
{
  vector<int *> *vec = new vector<int*>();
  int stop = all.size()-tsize+1;
  for(int i=0;i<stop;i++)
  {
    int *a1 = new int[tsize];
    memset(a1,0,tsize*sizeof(int));
    for(int j=0;j<tsize;j++)
      a1[j] = all[i+j];
      
    vec->push_back(a1);
  }
    
  return vec;
}


void PedData::mem_clean()
{
  for(int i=0;i<complete_geno_list.size();i++)
    delete complete_geno_list[i];
  complete_geno_list.clear();

  for(int i=0;i<geno_list.size();i++)
    delete geno_list[i];
  geno_list.clear();

  for(int i=0;i<missing_list.size();i++)
    delete missing_list[i];
  missing_list.clear();

  if(outoption==98){
  for(int i=0;i<geno_trio_list.size();i++)
    delete geno_trio_list[i];
  geno_trio_list.clear();

  for(int i=0;i<missing_trio_list.size();i++)
    delete missing_trio_list[i];
  missing_trio_list.clear();
  }
}


PedData::~PedData()
{
  for(int i=0;i<fam_list.size();i++)
    delete fam_list[i];
  fam_list.clear();

  fam_map.clear();
  hap_list.clear();
  hap_freq_table.clear();

  individual_list.clear();

  for(int i=0;i<snp_size;i++)
    delete[] content[i];
  content.clear();
}


int find_loc(const vector<trio> & s_vec, const trio & m)
{
  for(int i=0;i<s_vec.size();i++)
  {
    if(s_vec[i].cat_o==m.cat_o)
    {
      if(s_vec[i].cat_par1==m.cat_par1 || (s_vec[i].typed_par1==0 && m.typed_par1==0))
      {
        if(s_vec[i].cat_par2==m.cat_par2 || (s_vec[i].typed_par2==0 && m.typed_par2==0))
          return i;
      }
    }
  }
  return -1;
}


bool inc_sort(const punit& a, const punit& b)
{
  return a.p < b.p;
}


bool dec_sort(const hap_unit & a, const hap_unit & b)
{
  return a.freq > b.freq;
}