#include "iqlbestimator.h"


IQLBEstimator::IQLBEstimator(vector<Family *> & fam_list, vector<Genotypetrio *> & geno_trio_list, vector<Missingtrio *> & missing_trio_list, map<int, double> & hap_list, int size, int oneparent)
:fam_list(fam_list), geno_trio_list(geno_trio_list), missing_trio_list(missing_trio_list), hap_freq_table(hap_list), loci_num(size), oneparent(oneparent)
{
  sample_size=0;

  for(int i=0;i<fam_list.size();i++)
    sample_size+=fam_list[i]->N_typed;

  init_freq_table();
  compute_count();
  if(oneparent==1)
    compute_hst();
  geno_missing_trio_patten();
  pair_missing_pattern();

  flag_pd = 1;
  pvalue = 0;
}


void IQLBEstimator::init_freq_table()
{
  int tot;
  if(hap_freq_table.size()==0)
  {
    tot = 1<<loci_num;
    for(int i=0;i<tot;i++)
      hap_freq_table[i] = 1.0/double(tot);
  }
  else
  {
    tot = hap_freq_table.size();
    map<int, double>::iterator iter=hap_freq_table.begin();
    while(iter!=hap_freq_table.end())
    {
      iter->second = 1.0/double(tot);
      iter++;
    }
  }
  
  geno_freq_table.clear();
  map<int, double>::iterator iter1=hap_freq_table.begin();
  map<int, double>::iterator iter2;
  while(iter1!=hap_freq_table.end())
  {
    iter2 = iter1;
    while(iter2!=hap_freq_table.end())
    {
      int h1 = iter1->first;
      int h2 = iter2->first;

      if(h1 == h2)
        geno_freq_table[pair<int, int>(h1,h2)] = iter1->second*iter2->second;
      else
        geno_freq_table[pair<int, int>(h1,h2)] = 2*iter1->second*iter2->second;

      iter2++;
    }
    iter1++;
  }
}


void IQLBEstimator::reset_freq_table(map<int,double> & hap_freq)
{
  hap_freq_table.clear();
  hap_freq_table = hap_freq;
  
  geno_freq_table.clear();
  map<int, double>::iterator iter1=hap_freq_table.begin();
  map<int, double>::iterator iter2;
  while(iter1!=hap_freq_table.end())
  {
    iter2 = iter1;
    while(iter2!=hap_freq_table.end())
    {
      int h1 = iter1->first;
      int h2 = iter2->first;

      if(h1 == h2)
        geno_freq_table[pair<int, int>(h1,h2)] = iter1->second*iter2->second;
      else
        geno_freq_table[pair<int, int>(h1,h2)] = 2*iter1->second*iter2->second;

      iter2++;
    }
    iter1++;
  }
}


void IQLBEstimator::show_hap_freq()
{
/*
  map<int, double>::iterator iter = hap_freq_table.begin();
  while(iter!=hap_freq_table.end())
  {
    if(iter->second<0.00005)
    {
      iter++;
      continue;
    }
    printf("%2d ", iter->first);
    printf("   %.4f\n", iter->second);
    iter++;
  }
  printf("\n");

  map<int, double>::iterator iter = hap_freq_table.begin();
  while(iter!=hap_freq_table.end())
  {
    if(iter->second<0.00005)
    {
      iter++;
      continue;
    }
    fprintf(outfp, "%.4f\t", iter->second);
    iter++;
  }
*/

  fprintf(outfp, "Tag SNP haplotype frequency estimates in full sample\n");
  fprintf(outfp, "Haplotype\tFrequency\n");
  map<int, double>::iterator iter = hap_freq_table.begin();
//  iter = hap_freq_table.begin();
  while(iter!=hap_freq_table.end()){
    if(iter->second<0.00005){
      iter++;
      continue;
    }
    translate_hap(outfp,iter->first,loci_num);
    fprintf(outfp, " \t%.4f\n", iter->second);
    iter++;
  }
  fprintf(outfp, "\n");
}


void IQLBEstimator::show_hap_freq(map<int, int> &trans_map, vector<char> &a0v, vector<char> &a1v)
{
  fprintf(outfp, "Tag SNP haplotype frequency estimates in full sample\n");
  fprintf(outfp, "Haplotype\tFrequency\n");
  map<int, double>::iterator iter = hap_freq_table.begin();
  while(iter!=hap_freq_table.end()){
    if(iter->second<0.00005){
      iter++;
      continue;
    }
    translate_hap(outfp,trans_map[iter->first],loci_num,a0v,a1v);
    if(loci_num==2)
      fprintf(outfp, "\t");
    fprintf(outfp, " \t%.4f\n", iter->second);
    iter++;
  }
  fprintf(outfp, "\n");
}


void IQLBEstimator::show_geno_freq()
{
  map<pair<int, int>, double>::iterator giter = geno_freq_table.begin();
  while(giter!=geno_freq_table.end())
  {
    if(giter->second<0.00005)
    {
      giter++;
      continue;
    }
    printf("(%2d,%2d) ", (giter->first).first,(giter->first).second);
    printf("   %.4f\n", giter->second);
    giter++;
  }
  printf("\n");
}


void IQLBEstimator::compute_count()
{
  int row = hap_freq_table.size();
  int col = geno_freq_table.size();

  count = new int*[row-1];
  for(int i=0;i<row-1;i++)
  {
    count[i] = new int[col];
    memset(count[i],0,col*sizeof(int));
  }

  map<int, double>::iterator iter=hap_freq_table.begin();
  int loc = 0;
  while(iter!=hap_freq_table.end())
  {
    hap_map[iter->first] = loc;
    loc++;
    iter++;
  }

  loc = 0;
  map<pair<int, int>, double>::iterator giter=geno_freq_table.begin();
  while(giter!=geno_freq_table.end())
  {
    int h1=(giter->first).first;
    int h2=(giter->first).second;
    int loc1 = hap_map[h1];
    int loc2 = hap_map[h2];

    if(loc1!=row-1)
      count[loc1][loc]++;
    if(loc2!=row-1)
      count[loc2][loc]++;

    loc++;
    giter++;
  }
}


void IQLBEstimator::compute_hst()
{
  int row = hap_freq_table.size();
  int col = geno_freq_table.size();

  hst = new double**[row-1];
  for(int i=0;i<row-1;i++)
  {
    hst[i] = new double*[col];
    for(int j=0;j<col;j++)
    {
      hst[i][j] = new double[col];
      memset(hst[i][j],0,col*sizeof(double));
    }
  }

  int loc, loc1, loc2;
  map<pair<int, int>, double>::iterator giter1;
  map<pair<int, int>, double>::iterator giter2;

  loc1 = 0;
  giter1=geno_freq_table.begin();
  while(giter1!=geno_freq_table.end())
  {
    int i=(giter1->first).first;
    int j=(giter1->first).second;
    loc2 = 0;
    giter2=geno_freq_table.begin();
    while(giter2!=geno_freq_table.end())
    {
      int k=(giter2->first).first;
      int l=(giter2->first).second;

      if(k==l)  // parent is homozygote
      {
        if(i==k)
        {
          loc = hap_map[j];
          if(loc!=row-1)
            hst[loc][loc1][loc2]+=1;
          else
          {
            for(int s=0;s<row-1;s++)
              hst[s][loc1][loc2]-=1;
          }
        }
        else if(j==k)
        {
          loc = hap_map[i];
          if(loc!=row-1)
            hst[loc][loc1][loc2]+=1;
          else
          {
            for(int s=0;s<row-1;s++)
              hst[s][loc1][loc2]-=1;
          }
        }
      }
      else  // parent is heterozygote
      {
        if(i==j)  // offspring is homozygote
        {
          if(i==k || i==l)
          {
            loc = hap_map[j];
            if(loc!=row-1)
              hst[loc][loc1][loc2]+=0.5;
            else
            {
              for(int s=0;s<row-1;s++)
                hst[s][loc1][loc2]-=0.5;
            }
          }
        }
        else  // offspring is heterozygote
        {
          if(i==k && j==l)
          {
            loc = hap_map[i];
            if(loc!=row-1)
              hst[loc][loc1][loc2]+=0.5;
            else
            {
              for(int s=0;s<row-1;s++)
                hst[s][loc1][loc2]-=0.5;
            }
            loc = hap_map[j];
            if(loc!=row-1)
              hst[loc][loc1][loc2]+=0.5;
            else
            {
              for(int s=0;s<row-1;s++)
                hst[s][loc1][loc2]-=0.5;
            }
          }
          else if(i==k || i==l)
          {
            loc = hap_map[j];
            if(loc!=row-1)
              hst[loc][loc1][loc2]+=0.5;
            else
            {
              for(int s=0;s<row-1;s++)
                hst[s][loc1][loc2]-=0.5;
            }
          }
          else if(j==k || j==l)
          {
            loc = hap_map[i];
            if(loc!=row-1)
              hst[loc][loc1][loc2]+=0.5;
            else
            {
              for(int s=0;s<row-1;s++)
                hst[s][loc1][loc2]-=0.5;
            }
          }
        }
      }
      loc2++;
      giter2++;
    }
    loc1++;
    giter1++;
  }
}


void IQLBEstimator::geno_missing_trio_patten()
{
  int typed_o, typed_par1, typed_par2;

  for(int i=0;i<geno_trio_list.size();i++)
  {
    typed_o = geno_trio_list[i]->typed_o;
    typed_par1 = geno_trio_list[i]->typed_par1;
    typed_par2 = geno_trio_list[i]->typed_par2;

    if(typed_par1==1 || typed_par2==1)
    {
      int *a_o=(geno_trio_list[i]->offspring)->get_gcode();
      int *a_f=(geno_trio_list[i]->par1)->get_gcode();
      int *a_m=(geno_trio_list[i]->par2)->get_gcode();

      for(int j=0;j<missing_trio_list.size();j++)
      {
        int flag;

        if(typed_o!=missing_trio_list[j]->typed_o || typed_par1!=missing_trio_list[j]->typed_par1 || typed_par2!=missing_trio_list[j]->typed_par2)
        {
          flag = 0;
          continue;
        }

        int *b_o=(missing_trio_list[j]->M_offspring)->get_gcode();
        int *b_f=(missing_trio_list[j]->M_par1)->get_gcode();
        int *b_m=(missing_trio_list[j]->M_par2)->get_gcode();

        flag = compare_equ(a_o,b_o,loci_num);
        if(flag==0)
          continue;

        if((compare_equ(a_f,b_f,loci_num)==1 && compare_equ(a_m,b_m,loci_num)==1) || (compare_equ(a_f,b_m,loci_num)==1 && compare_equ(a_m,b_f,loci_num)==1))
        {
          geno_missing_trio_map[i] = j;
          break;
        }
      }
    }
    else if(typed_par1==0 && typed_par2==0)
    {
      int *a_o=(geno_trio_list[i]->offspring)->get_gcode();

      for(int j=0;j<missing_trio_list.size();j++)
      {
        int flag=1;
        if(typed_o!=missing_trio_list[j]->typed_o || typed_par1!=missing_trio_list[j]->typed_par1 || typed_par2!=missing_trio_list[j]->typed_par2)
        {
          flag = 0;
          continue;
        }

        int *b_o=(missing_trio_list[j]->M_offspring)->get_gcode();

        flag = compare_equ(a_o,b_o,loci_num);
        if(flag==0)
          continue;

        geno_missing_trio_map[i] = j;
        break;
      }
    }
  }
/*
  if(geno_missing_trio_map.size()!=geno_trio_list.size())
  {
    printf("geno_missing_trio_map assignment is wrong. Please check.\n");
    exit(1);
  }
*/
}

  
void IQLBEstimator::show_test(double num)
{
/*
  printf("the M-type IQLS_a test stat is %.4f\n", test);
  fprintf(outfp, "%.4f\t", test);
*/
  fprintf(outfp, "ATRIUM test\n");
  if(flag_pd==0){
    test = num;
    pvalue = pochisq(num,1);
  }
  fprintf(outfp, "ATRIUM statistic value = %.6f\t pvalue = %g\n", test,pvalue);
}


void IQLBEstimator::update_hap_freq()
{
  map<int, double>::iterator iter = hap_freq_table.begin();
  while(iter!=hap_freq_table.end())
  {
    iter->second = 0;
    iter++;
  }
    
  map<pair<int, int>, double>::iterator giter = geno_freq_table.begin();
  
  while(giter!= geno_freq_table.end())
  {
    pair<int, int> conf = giter->first;
    if(conf.first==conf.second)
      hap_freq_table[conf.first] += geno_freq_table[conf];
    else
    {
      hap_freq_table[conf.first] += .5*geno_freq_table[conf];
      hap_freq_table[conf.second] += .5*geno_freq_table[conf];
    }
    giter++;
  }
}


void IQLBEstimator::update_geno_freq()
{
  map<pair<int, int>, double>::iterator giter = geno_freq_table.begin();
  while(giter!=geno_freq_table.end())
  {
    giter->second = 0;
    giter++;
  }

  giter = geno_freq_table.begin();
  while(giter!=geno_freq_table.end())
  {
    pair<int, int> conf = giter->first;
    int k1 = conf.first;
    int k2 = conf.second;
    if(k1==k2)
      giter->second = hap_freq_table[k1]*hap_freq_table[k2];
    else
      giter->second = 2*hap_freq_table[k1]*hap_freq_table[k2];
    
    giter++;
  }
}


IQLBEstimator::~IQLBEstimator()
{
  int row = hap_freq_table.size();
  int col = geno_freq_table.size();

  for(int k=0;k<row-1;k++)
    delete[] count[k];
  delete[] count;

  if(oneparent==1)
  {
    for(int k=0;k<row-1;k++)
    {
      for(int l=0;l<col;l++)
        delete[] hst[k][l];
      delete[] hst[k];
    }
    delete[] hst;
  }

  hap_freq_table.clear();
  geno_freq_table.clear();
  geno_missing_trio_map.clear();
  hap_map.clear();
}


// find all possible pair of missing pattern
void IQLBEstimator::pair_missing_pattern(){

  int ind_1, ind_2, ind_1f, ind_1m, ind_2f, ind_2m;
  int par_1, par_2;  // how many parents is typed
  int dim_1, dim_2;
  int trio_cat_1, trio_cat_2;
  int mis_trio_cat_1, mis_trio_cat_2;
  double IBD_2, IBD_1, kin;
  map<pair<int, int>, double>::iterator iterk;
  map<pair<int, int>, double>::iterator iter1;
  map<pair<int, int>, double>::iterator iter2;
          
  for(int fam=0;fam<fam_list.size();fam++){
    if(fam_list[fam]->N_typed>0){
      for(int i=0;i<fam_list[fam]->N;i++){
        trio_cat_1 = fam_list[fam]->trio_cat[i];
        mis_trio_cat_1 = fam_list[fam]->mis_trio_cat[i];
        dim_1 = missing_trio_list[mis_trio_cat_1]->dim;
        if(dim_1 > 0 && geno_trio_list[trio_cat_1]->info==1){
          ind_1 = fam_list[fam]->fam_member[i][0];
          par_1 = 0;
          par_1+=missing_trio_list[mis_trio_cat_1]->typed_par1;
          par_1+=missing_trio_list[mis_trio_cat_1]->typed_par2;
          ind_1f = fam_list[fam]->fam_member[i][1];
          ind_1m = fam_list[fam]->fam_member[i][2];

          for(int j=i+1;j<fam_list[fam]->N;j++){
            trio_cat_2 = fam_list[fam]->trio_cat[j];
            mis_trio_cat_2 = fam_list[fam]->mis_trio_cat[j];
            dim_2 = missing_trio_list[mis_trio_cat_2]->dim;
            if(dim_2 > 0 && geno_trio_list[trio_cat_2]->info==1){
              ind_2 = fam_list[fam]->fam_member[j][0];
              iterk = (fam_list[fam]->kinship_coeff).find(pair<int, int>(ind_1,ind_2));
              if(iterk == (fam_list[fam]->kinship_coeff).end()){
                printf("No IBD coefficient between individual %d and individual %d from family %d. Please check...\n\n", ind_1,ind_2,fam_list[fam]->fam_id);
                exit(1);
              }

              if(iterk->second>1e-8){ // i is related to j
                par_2 = 0;
                par_2+=missing_trio_list[mis_trio_cat_2]->typed_par1;
                par_2+=missing_trio_list[mis_trio_cat_2]->typed_par2;
                ind_2f = fam_list[fam]->fam_member[j][1];
                ind_2m = fam_list[fam]->fam_member[j][2];

                kin = iterk->second;
                iter2 = (fam_list[fam]->IBD2_coeff).find(pair<int, int>(ind_1,ind_2));
                iter1 = (fam_list[fam]->IBD1_coeff).find(pair<int, int>(ind_1,ind_2));
                IBD_2 = iter2->second;
                IBD_1 = iter1->second;
                
                if(par_1==0 && par_2==0){
                  if(fam_list[fam]->IBD1_coeff[pair<int, int>(ind_1,ind_2)]>1e-8){
                    if(off_one.find(mis_trio_cat_1)==off_one.end())
                      off_one[mis_trio_cat_1] = 0;
                    if(off_one.find(mis_trio_cat_2)==off_one.end())
                      off_one[mis_trio_cat_2] = 0;
                  }
                  if(fam_list[fam]->IBD2_coeff[pair<int, int>(ind_1,ind_2)]>1e-8){
                    if(off_two.find(mis_trio_cat_1)==off_two.end())
                      off_two[mis_trio_cat_1] = 0;
                    if(off_two.find(mis_trio_cat_2)==off_two.end())
                      off_two[mis_trio_cat_2] = 0;
                  }
                }
                else if(par_1==0 && (par_2==1 || par_2==2)){
                  int order = fam_list[fam]->mis_par_order[j];
                  int ind_2par1 = ind_2f;
                  int ind_2par2 = ind_2m;
                  if(order==1){
                    ind_2par1 = ind_2m;
                    ind_2par2 = ind_2f;
                  }  
                  int diff = 0;
                  int index = relation(fam,ind_1,ind_2,ind_2par1,ind_2par2,kin,IBD_2,IBD_1,diff);
                  if(index==4){
                    if(off_two.find(mis_trio_cat_1)==off_two.end())
                      off_two[mis_trio_cat_1] = 0;
                    if(par1_two.find(mis_trio_cat_2)==par1_two.end())
                      par1_two[mis_trio_cat_2] = 0;           
                  }
                  else if(index==5){
                    if(off_two.find(mis_trio_cat_1)==off_two.end())
                      off_two[mis_trio_cat_1] = 0;
                    if(par2_two.find(mis_trio_cat_2)==par2_two.end())
                      par2_two[mis_trio_cat_2] = 0;   
                  }
                  else if(index==1){
                    if(fam_list[fam]->IBD1_coeff[pair<int, int>(ind_1,ind_2par1)]>1e-8){
                      if(off_one.find(mis_trio_cat_1)==off_one.end())
                        off_one[mis_trio_cat_1] = 0;
                      if(par1_one.find(mis_trio_cat_2)==par1_one.end())
                        par1_one[mis_trio_cat_2] = 0;
                    }
                    if(fam_list[fam]->IBD2_coeff[pair<int, int>(ind_1,ind_2par1)]>1e-8){
                      if(off_two.find(mis_trio_cat_1)==off_two.end())
                        off_two[mis_trio_cat_1] = 0;
                      if(par1_two.find(mis_trio_cat_2)==par1_two.end())
                        par1_two[mis_trio_cat_2] = 0;
                    }
                  }
                  else if(index==2){
                    if(fam_list[fam]->IBD1_coeff[pair<int, int>(ind_1,ind_2par2)]>1e-8){
                      if(off_one.find(mis_trio_cat_1)==off_one.end())
                        off_one[mis_trio_cat_1] = 0;
                      if(par2_one.find(mis_trio_cat_2)==par2_one.end())
                        par2_one[mis_trio_cat_2] = 0;
                    }
                    if(fam_list[fam]->IBD2_coeff[pair<int, int>(ind_1,ind_2par2)]>1e-8){
                      if(off_two.find(mis_trio_cat_1)==off_two.end())
                        off_two[mis_trio_cat_1] = 0;
                      if(par2_two.find(mis_trio_cat_2)==par2_two.end())
                        par2_two[mis_trio_cat_2] = 0;
                    }   
                  }
                  else if(index==3){
                    if(fam_list[fam]->IBD1_coeff[pair<int, int>(ind_1,ind_2)]>1e-8){
                      if(off_one.find(mis_trio_cat_1)==off_one.end())
                        off_one[mis_trio_cat_1] = 0;
                      if(off_one.find(mis_trio_cat_2)==off_one.end())
                        off_one[mis_trio_cat_2] = 0;
                    }
                    if(fam_list[fam]->IBD2_coeff[pair<int, int>(ind_1,ind_2)]>1e-8){
                      if(off_two.find(mis_trio_cat_1)==off_two.end())
                        off_two[mis_trio_cat_1] = 0;
                      if(off_two.find(mis_trio_cat_2)==off_two.end())
                        off_two[mis_trio_cat_2] = 0;
                    }   
                  }
                  else if(index==6){  // can not be sibs
                    if(diff==1)
                      printf("please check relationships\n");
                    if(off_one.find(mis_trio_cat_1)==off_one.end())
                        off_one[mis_trio_cat_1] = 0;
                    if(par1_one.find(mis_trio_cat_2)==par1_one.end())
                        par1_one[mis_trio_cat_2] = 0;
                    if(par2_one.find(mis_trio_cat_2)==par2_one.end())
                        par2_one[mis_trio_cat_2] = 0;             
                  }
                }
                else if((par_1==1 || par_1==2) && par_2==0){
                  int order = fam_list[fam]->mis_par_order[i];
                  int ind_1par1 = ind_1f;
                  int ind_1par2 = ind_1m;
                  if(order==1){
                    ind_1par1 = ind_1m;
                    ind_1par2 = ind_1f;
                  }  
                  int diff = 0;
                  int index = relation(fam,ind_2,ind_1,ind_1par1,ind_1par2,kin,IBD_2,IBD_1,diff);
                  if(index==4){
                    if(off_two.find(mis_trio_cat_2)==off_two.end())
                      off_two[mis_trio_cat_2] = 0;
                    if(par1_two.find(mis_trio_cat_1)==par1_two.end())
                      par1_two[mis_trio_cat_1] = 0;           
                  }
                  else if(index==5){
                    if(off_two.find(mis_trio_cat_2)==off_two.end())
                      off_two[mis_trio_cat_2] = 0;
                    if(par2_two.find(mis_trio_cat_1)==par2_two.end())
                      par2_two[mis_trio_cat_1] = 0;   
                  }                  
                  else if(index==1){
                    if(fam_list[fam]->IBD1_coeff[pair<int, int>(ind_1par1,ind_2)]>1e-8){
                      if(off_one.find(mis_trio_cat_2)==off_one.end())
                        off_one[mis_trio_cat_2] = 0;
                      if(par1_one.find(mis_trio_cat_1)==par1_one.end())
                        par1_one[mis_trio_cat_1] = 0;
                    }
                    if(fam_list[fam]->IBD2_coeff[pair<int, int>(ind_1par1,ind_2)]>1e-8){
                      if(off_two.find(mis_trio_cat_2)==off_two.end())
                        off_two[mis_trio_cat_2] = 0;
                      if(par1_two.find(mis_trio_cat_1)==par1_two.end())
                        par1_two[mis_trio_cat_1] = 0;
                    }
                  }
                  else if(index==2){
                    if(fam_list[fam]->IBD1_coeff[pair<int, int>(ind_1par2,ind_2)]>1e-8){
                      if(off_one.find(mis_trio_cat_2)==off_one.end())
                        off_one[mis_trio_cat_2] = 0;
                      if(par2_one.find(mis_trio_cat_1)==par2_one.end())
                        par2_one[mis_trio_cat_1] = 0;
                    }
                    if(fam_list[fam]->IBD2_coeff[pair<int, int>(ind_1par2,ind_2)]>1e-8){
                      if(off_two.find(mis_trio_cat_2)==off_two.end())
                        off_two[mis_trio_cat_2] = 0;
                      if(par2_two.find(mis_trio_cat_1)==par2_two.end())
                        par2_two[mis_trio_cat_1] = 0;
                    }   
                  }
                  else if(index==3){
                    if(fam_list[fam]->IBD1_coeff[pair<int, int>(ind_1,ind_2)]>1e-8){
                      if(off_one.find(mis_trio_cat_1)==off_one.end())
                        off_one[mis_trio_cat_1] = 0;
                      if(off_one.find(mis_trio_cat_2)==off_one.end())
                        off_one[mis_trio_cat_2] = 0;
                    }
                    if(fam_list[fam]->IBD2_coeff[pair<int, int>(ind_1,ind_2)]>1e-8){
                      if(off_two.find(mis_trio_cat_1)==off_two.end())
                        off_two[mis_trio_cat_1] = 0;
                      if(off_two.find(mis_trio_cat_2)==off_two.end())
                        off_two[mis_trio_cat_2] = 0;
                    }   
                  }
                  else if(index==6){  // can not be sibs
                    if(diff==1)
                      printf("please check relationships\n");
                    if(off_one.find(mis_trio_cat_2)==off_one.end())
                        off_one[mis_trio_cat_2] = 0;
                    if(par1_one.find(mis_trio_cat_1)==par1_one.end())
                        par1_one[mis_trio_cat_1] = 0;
                    if(par2_one.find(mis_trio_cat_1)==par2_one.end())
                        par2_one[mis_trio_cat_1] = 0;             
                  }
                }                
                else if((par_1==1 || par_1==2) && (par_2==1 || par_2==2)){
                  int order_1 = fam_list[fam]->mis_par_order[i];
                  int ind_1par1 = ind_1f;
                  int ind_1par2 = ind_1m;
                  if(order_1==1){
                    ind_1par1 = ind_1m;
                    ind_1par2 = ind_1f;
                  }
                  int order_2 = fam_list[fam]->mis_par_order[j];
                  int ind_2par1 = ind_2f;
                  int ind_2par2 = ind_2m;
                  if(order_2==1){
                    ind_2par1 = ind_2m;
                    ind_2par2 = ind_2f;
                  }                    
                  int diff = 0;
                  int index = relation(fam,ind_1,ind_1par1,ind_1par2,ind_2,ind_2par1,ind_2par2,kin,IBD_2,IBD_1,diff);
                  if(index==13){
                    if(off_two.find(mis_trio_cat_1)==off_two.end())
                      off_two[mis_trio_cat_1] = 0;
                    if(par1_two.find(mis_trio_cat_2)==par1_two.end())
                      par1_two[mis_trio_cat_2] = 0;              
                  }
                  else if(index==14){
                    if(off_two.find(mis_trio_cat_1)==off_two.end())
                      off_two[mis_trio_cat_1] = 0;
                    if(par2_two.find(mis_trio_cat_2)==par2_two.end())
                      par2_two[mis_trio_cat_2] = 0;   
                  }
                  else if(index==15){
                    if(off_two.find(mis_trio_cat_2)==off_two.end())
                      off_two[mis_trio_cat_2] = 0;
                    if(par1_two.find(mis_trio_cat_1)==par1_two.end())
                      par1_two[mis_trio_cat_1] = 0;              
                  }
                  else if(index==16){
                    if(off_two.find(mis_trio_cat_2)==off_two.end())
                      off_two[mis_trio_cat_2] = 0;
                    if(par2_two.find(mis_trio_cat_1)==par2_two.end())
                      par2_two[mis_trio_cat_1] = 0;   
                  }
                  else if(index==17){
                    if(parents_two.find(mis_trio_cat_1)==parents_two.end())
                      parents_two[mis_trio_cat_1] = 0;
                    if(parents_two.find(mis_trio_cat_2)==parents_two.end())
                      parents_two[mis_trio_cat_2] = 0;     
                  }
                  else if(index==18){
                    if(ind_1par1==ind_2par1){
                      if(par1_two.find(mis_trio_cat_1)==par1_two.end())
                        par1_two[mis_trio_cat_1] = 0;
                      if(par1_two.find(mis_trio_cat_2)==par1_two.end())
                        par1_two[mis_trio_cat_2] = 0;                   
                    }
                    else if(ind_1par1==ind_2par2){
                      if(par1_two.find(mis_trio_cat_1)==par1_two.end())
                        par1_two[mis_trio_cat_1] = 0;
                      if(par2_two.find(mis_trio_cat_2)==par2_two.end())
                        par2_two[mis_trio_cat_2] = 0;   
                    }
                    else if(ind_1par2==ind_2par1){
                      if(par2_two.find(mis_trio_cat_1)==par2_two.end())
                        par2_two[mis_trio_cat_1] = 0;
                      if(par1_two.find(mis_trio_cat_2)==par1_two.end())
                        par1_two[mis_trio_cat_2] = 0;   
                    }
                    else if(ind_1par2==ind_2par2){
                      if(par2_two.find(mis_trio_cat_1)==par2_two.end())
                        par2_two[mis_trio_cat_1] = 0;
                      if(par2_two.find(mis_trio_cat_2)==par2_two.end())
                        par2_two[mis_trio_cat_2] = 0;   
                    }     
                  }
                  else if(index==1){
                    if(fam_list[fam]->IBD1_coeff[pair<int, int>(ind_1par1,ind_2par1)]>1e-8){
                      if(par1_one.find(mis_trio_cat_1)==par1_one.end())
                        par1_one[mis_trio_cat_1] = 0;
                      if(par1_one.find(mis_trio_cat_2)==par1_one.end())
                        par1_one[mis_trio_cat_2] = 0;
                    }
                    if(fam_list[fam]->IBD2_coeff[pair<int, int>(ind_1par1,ind_2par1)]>1e-8){
                      if(par1_two.find(mis_trio_cat_1)==par1_two.end())
                        par1_two[mis_trio_cat_1] = 0;
                      if(par1_two.find(mis_trio_cat_2)==par1_two.end())
                        par1_two[mis_trio_cat_2] = 0;
                    }   
                  }
                  else if(index==2){
                    if(fam_list[fam]->IBD1_coeff[pair<int, int>(ind_1par1,ind_2par2)]>1e-8){
                      if(par1_one.find(mis_trio_cat_1)==par1_one.end())
                        par1_one[mis_trio_cat_1] = 0;
                      if(par2_one.find(mis_trio_cat_2)==par2_one.end())
                        par2_one[mis_trio_cat_2] = 0;
                    }
                    if(fam_list[fam]->IBD2_coeff[pair<int, int>(ind_1par1,ind_2par2)]>1e-8){
                      if(par1_two.find(mis_trio_cat_1)==par1_two.end())
                        par1_two[mis_trio_cat_1] = 0;
                      if(par2_two.find(mis_trio_cat_2)==par2_two.end())
                        par2_two[mis_trio_cat_2] = 0;
                    }   
                  }
                  else if(index==3){
                    if(fam_list[fam]->IBD1_coeff[pair<int, int>(ind_1par2,ind_2par1)]>1e-8){
                      if(par2_one.find(mis_trio_cat_1)==par2_one.end())
                        par2_one[mis_trio_cat_1] = 0;
                      if(par1_one.find(mis_trio_cat_2)==par1_one.end())
                        par1_one[mis_trio_cat_2] = 0;
                    }
                    if(fam_list[fam]->IBD2_coeff[pair<int, int>(ind_1par2,ind_2par1)]>1e-8){
                      if(par2_two.find(mis_trio_cat_1)==par2_two.end())
                        par2_two[mis_trio_cat_1] = 0;
                      if(par1_two.find(mis_trio_cat_2)==par1_two.end())
                        par1_two[mis_trio_cat_2] = 0;
                    }   
                  }   
                  else if(index==4){
                    if(fam_list[fam]->IBD1_coeff[pair<int, int>(ind_1par2,ind_2par2)]>1e-8){
                      if(par2_one.find(mis_trio_cat_1)==par2_one.end())
                        par2_one[mis_trio_cat_1] = 0;
                      if(par2_one.find(mis_trio_cat_2)==par2_one.end())
                        par2_one[mis_trio_cat_2] = 0;
                    }
                    if(fam_list[fam]->IBD2_coeff[pair<int, int>(ind_1par2,ind_2par2)]>1e-8){
                      if(par2_two.find(mis_trio_cat_1)==par2_two.end())
                        par2_two[mis_trio_cat_1] = 0;
                      if(par2_two.find(mis_trio_cat_2)==par2_two.end())
                        par2_two[mis_trio_cat_2] = 0;
                    }   
                  }
                  else if(index==5){
                    if(fam_list[fam]->IBD1_coeff[pair<int, int>(ind_1par1,ind_2)]>1e-8){
                      if(par1_one.find(mis_trio_cat_1)==par1_one.end())
                        par1_one[mis_trio_cat_1] = 0;
                      if(off_one.find(mis_trio_cat_2)==off_one.end())
                        off_one[mis_trio_cat_2] = 0;
                    }
                    if(fam_list[fam]->IBD2_coeff[pair<int, int>(ind_1par1,ind_2)]>1e-8){
                      if(par1_two.find(mis_trio_cat_1)==par1_two.end())
                        par1_two[mis_trio_cat_1] = 0;
                      if(off_two.find(mis_trio_cat_2)==off_two.end())
                        off_two[mis_trio_cat_2] = 0;
                    }     
                  }
                  else if(index==6){
                    if(fam_list[fam]->IBD1_coeff[pair<int, int>(ind_1par2,ind_2)]>1e-8){
                      if(par2_one.find(mis_trio_cat_1)==par2_one.end())
                        par2_one[mis_trio_cat_1] = 0;
                      if(off_one.find(mis_trio_cat_2)==off_one.end())
                        off_one[mis_trio_cat_2] = 0;
                    }
                    if(fam_list[fam]->IBD2_coeff[pair<int, int>(ind_1par2,ind_2)]>1e-8){
                      if(par2_two.find(mis_trio_cat_1)==par2_two.end())
                        par2_two[mis_trio_cat_1] = 0;
                      if(off_two.find(mis_trio_cat_2)==off_two.end())
                        off_two[mis_trio_cat_2] = 0;
                    }     
                  }
                  else if(index==7){
                    if(fam_list[fam]->IBD1_coeff[pair<int, int>(ind_1,ind_2par1)]>1e-8){
                      if(off_one.find(mis_trio_cat_1)==off_one.end())
                        off_one[mis_trio_cat_1] = 0;
                      if(par1_one.find(mis_trio_cat_2)==par1_one.end())
                        par1_one[mis_trio_cat_2] = 0;                      
                    }
                    if(fam_list[fam]->IBD2_coeff[pair<int, int>(ind_1,ind_2par1)]>1e-8){
                      if(off_two.find(mis_trio_cat_1)==off_two.end())
                        off_two[mis_trio_cat_1] = 0;
                      if(par1_two.find(mis_trio_cat_2)==par1_two.end())
                        par1_two[mis_trio_cat_2] = 0;                      
                    }     
                  }
                  else if(index==8){
                    if(fam_list[fam]->IBD1_coeff[pair<int, int>(ind_1,ind_2par2)]>1e-8){
                      if(off_one.find(mis_trio_cat_1)==off_one.end())
                        off_one[mis_trio_cat_1] = 0;
                      if(par2_one.find(mis_trio_cat_2)==par2_one.end())
                        par2_one[mis_trio_cat_2] = 0;                      
                    }
                    if(fam_list[fam]->IBD2_coeff[pair<int, int>(ind_1,ind_2par2)]>1e-8){
                      if(off_two.find(mis_trio_cat_1)==off_two.end())
                        off_two[mis_trio_cat_1] = 0;
                      if(par2_two.find(mis_trio_cat_2)==par2_two.end())
                        par2_two[mis_trio_cat_2] = 0;                      
                    }     
                  }
                  else if(index==9){
                    if(diff==1){  // j and if are sibs
                      if(par1_two.find(mis_trio_cat_1)==par1_two.end())
                        par1_two[mis_trio_cat_1] = 0;
                      if(parents_one.find(mis_trio_cat_2)==parents_one.end())
                        parents_one[mis_trio_cat_2] = 0;          
                    }
                    else{
                      if(par1_one.find(mis_trio_cat_1)==par1_one.end())
                        par1_one[mis_trio_cat_1] = 0;
                      if(par1_one.find(mis_trio_cat_2)==par1_one.end())
                        par1_one[mis_trio_cat_2] = 0;
                      if(par2_one.find(mis_trio_cat_2)==par2_one.end())
                        par2_one[mis_trio_cat_2] = 0;   
                    }                          
                  }
                  else if(index==10){
                    if(diff==1){  // j and im are sibs
                      if(par2_two.find(mis_trio_cat_1)==par2_two.end())
                        par2_two[mis_trio_cat_1] = 0;
                      if(parents_one.find(mis_trio_cat_2)==parents_one.end())
                        parents_one[mis_trio_cat_2] = 0;          
                    }
                    else{
                      if(par2_one.find(mis_trio_cat_1)==par2_one.end())
                        par2_one[mis_trio_cat_1] = 0;
                      if(par1_one.find(mis_trio_cat_2)==par1_one.end())
                        par1_one[mis_trio_cat_2] = 0;
                      if(par2_one.find(mis_trio_cat_2)==par2_one.end())
                        par2_one[mis_trio_cat_2] = 0;   
                    }                          
                  }
                  else if(index==11){
                    if(diff==1){  // i and jf are sibs
                      if(par1_two.find(mis_trio_cat_2)==par1_two.end())
                        par1_two[mis_trio_cat_2] = 0;
                      if(parents_one.find(mis_trio_cat_1)==parents_one.end())
                        parents_one[mis_trio_cat_1] = 0;          
                    }
                    else{
                      if(par1_one.find(mis_trio_cat_2)==par1_one.end())
                        par1_one[mis_trio_cat_2] = 0;
                      if(par1_one.find(mis_trio_cat_1)==par1_one.end())
                        par1_one[mis_trio_cat_1] = 0;
                      if(par2_one.find(mis_trio_cat_1)==par2_one.end())
                        par2_one[mis_trio_cat_1] = 0;   
                    }                          
                  }
                  else if(index==12){
                    if(diff==1){  // i and jm are sibs
                      if(par2_two.find(mis_trio_cat_2)==par2_two.end())
                        par2_two[mis_trio_cat_2] = 0;
                      if(parents_one.find(mis_trio_cat_1)==parents_one.end())
                        parents_one[mis_trio_cat_1] = 0;          
                    }
                    else{
                      if(par2_one.find(mis_trio_cat_2)==par2_one.end())
                        par2_one[mis_trio_cat_2] = 0;
                      if(par1_one.find(mis_trio_cat_1)==par1_one.end())
                        par1_one[mis_trio_cat_1] = 0;
                      if(par2_one.find(mis_trio_cat_1)==par2_one.end())
                        par2_one[mis_trio_cat_1] = 0;   
                    }                          
                  }
                }
              }  // i, j related
            }  // j is typed
          }  // complete j
        }  // i is typed
      }  // complete i
    }
  }
}


void IQLBEstimator::allocate(){

  map<pair<int, int>, double>::iterator giter=geno_freq_table.begin();
  int loc = 0;
  while(giter!=geno_freq_table.end()){
    hap_pair_map[giter->first] = loc;
    loc++;
    giter++;
  }

  int h_size = hap_freq_table.size();
  int col = geno_freq_table.size();

  rst = new double*[h_size-1];
  for(int i=0;i<h_size-1;i++){
    rst[i] = new double[col];
    memset(rst[i],0,col*sizeof(double));
  }  // record derivative for each pair of haplotypes

  freq = new double[h_size];
  memset(freq,0,h_size*sizeof(double));
  
  genofreq = new double[col];
  memset(genofreq,0,col*sizeof(double));

  int par1_one_size = par1_one.size();
  if(par1_one_size>0){
    result_par1_one = new double **[par1_one_size];
    map<int, int>::iterator pi=par1_one.begin();
    for(int i=0;i<par1_one_size;i++){
      pi->second = i;
      int mis = pi->first;
      int dim_1 = missing_trio_list[mis]->dim;
      result_par1_one[i] = new double*[dim_1];
      for(int j=0;j<dim_1;j++){
        result_par1_one[i][j] = new double[h_size];
        memset(result_par1_one[i][j],0,h_size*sizeof(double));
      }
      pi++;
    }
  }
  
  int par1_two_size = par1_two.size();
  if(par1_two_size>0){
    result_par1_two = new double **[par1_two_size];
    map<int, int>::iterator pi=par1_two.begin();
    for(int i=0;i<par1_two_size;i++){
      pi->second = i;
      int mis = pi->first;
      int dim_1 = missing_trio_list[mis]->dim;
      result_par1_two[i] = new double*[dim_1];
      for(int j=0;j<dim_1;j++){
        result_par1_two[i][j] = new double[col];
        memset(result_par1_two[i][j],0,col*sizeof(double));
      }
      pi++;
    }
  }
  
  int par2_one_size = par2_one.size();
  if(par2_one_size>0){
    result_par2_one = new double **[par2_one_size];
    map<int, int>::iterator pi=par2_one.begin();
    for(int i=0;i<par2_one_size;i++){
      pi->second = i;
      int mis = pi->first;
      int dim_1 = missing_trio_list[mis]->dim;
      result_par2_one[i] = new double*[dim_1];
      for(int j=0;j<dim_1;j++){
        result_par2_one[i][j] = new double[h_size];
        memset(result_par2_one[i][j],0,h_size*sizeof(double));
      }
      pi++;
    }
  }
  
  int par2_two_size = par2_two.size();
  if(par2_two_size>0){
    result_par2_two = new double **[par2_two_size];
    map<int, int>::iterator pi=par2_two.begin();
    for(int i=0;i<par2_two_size;i++){
      pi->second = i;
      int mis = pi->first;
      int dim_1 = missing_trio_list[mis]->dim;
      result_par2_two[i] = new double*[dim_1];
      for(int j=0;j<dim_1;j++){
        result_par2_two[i][j] = new double[col];
        memset(result_par2_two[i][j],0,col*sizeof(double));
      }
      pi++;
    }
  }  

  int off_one_size = off_one.size();
  if(off_one_size>0){
    result_off_one = new double **[off_one_size];
    map<int, int>::iterator pi=off_one.begin();
    for(int i=0;i<off_one_size;i++){
      pi->second = i;
      int mis = pi->first;
      int dim_1 = missing_trio_list[mis]->dim;
      result_off_one[i] = new double*[dim_1];
      for(int j=0;j<dim_1;j++){
        result_off_one[i][j] = new double[h_size];
        memset(result_off_one[i][j],0,h_size*sizeof(double));
      }
      pi++;
    }
  }
  
  int off_two_size = off_two.size();
  if(off_two_size>0){
    result_off_two = new double **[off_two_size];
    map<int, int>::iterator pi=off_two.begin();
    for(int i=0;i<off_two_size;i++){
      pi->second = i;
      int mis = pi->first;
      int dim_1 = missing_trio_list[mis]->dim;
      result_off_two[i] = new double*[dim_1];
      for(int j=0;j<dim_1;j++){
        result_off_two[i][j] = new double[col];
        memset(result_off_two[i][j],0,col*sizeof(double));
      }
      pi++;
    }
  }  

  int parents_one_size = parents_one.size();
  if(parents_one_size>0){
    result_parents_one = new double ***[parents_one_size];
    map<int, int>::iterator pi=parents_one.begin();
    for(int i=0;i<parents_one_size;i++){
      pi->second = i;
      int mis = pi->first;
      int dim_1 = missing_trio_list[mis]->dim;
      result_parents_one[i] = new double**[dim_1];
      for(int j=0;j<dim_1;j++){
        result_parents_one[i][j] = new double*[h_size];
        for(int k=0;k<h_size;k++){
          result_parents_one[i][j][k] = new double[h_size];
          memset(result_parents_one[i][j][k],0,h_size*sizeof(double));
        }  
      }
      pi++;
    }
  }
  
  int parents_two_size = parents_two.size();
  if(parents_two_size>0){
    result_parents_two = new double ***[parents_two_size];
    map<int, int>::iterator pi=parents_two.begin();
    for(int i=0;i<parents_two_size;i++){
      pi->second = i;
      int mis = pi->first;
      int dim_1 = missing_trio_list[mis]->dim;
      result_parents_two[i] = new double**[dim_1];
      for(int j=0;j<dim_1;j++){
        result_parents_two[i][j] = new double*[col];
        for(int k=0;k<col;k++){
          result_parents_two[i][j][k] = new double[col];
          memset(result_parents_two[i][j][k],0,col*sizeof(double));
        }
      }
      pi++;
    }
  }

  denominator = new double*[h_size-1];
  for(int i=0;i<h_size-1;i++){
    denominator[i] = new double[h_size-1];
    memset(denominator[i],0,(h_size-1)*sizeof(double));
  }  // F_p^T Omega^(-1) F_p

  numerator = new double[h_size-1];
  memset(numerator,0,(h_size-1)*sizeof(double));  // F_p^T Omega^(-1) (Z-\mu)

  tot_dim = 0;
  for(int fam=0;fam<fam_list.size();fam++){
    if(fam_list[fam]->N_typed>0){
      int trio_cat;
      int mis_trio_cat;
      int dim=0;
      for(int i=0;i<fam_list[fam]->N;i++){
        trio_cat = fam_list[fam]->trio_cat[i];
        if(geno_trio_list[trio_cat]->info==1){
          mis_trio_cat = fam_list[fam]->mis_trio_cat[i];
          dim+=missing_trio_list[mis_trio_cat]->dim;
        }
      }

      if(dim>tot_dim)
        tot_dim = dim;
    }
  }

  cholent = new double*[tot_dim];
  for(int i=0;i<tot_dim;i++){
    cholent[i] = new double[h_size];
    memset(cholent[i],0,h_size*sizeof(double));
  }

  chol = new double*[tot_dim];
  for(int i=0;i<tot_dim;i++){
    chol[i] = new double[tot_dim];
    memset(chol[i],0,tot_dim*sizeof(double));
  }

  cholaug = new double*[tot_dim];
  for(int i=0;i<tot_dim;i++){
    cholaug[i] = new double[h_size];
    memset(cholaug[i],0,h_size*sizeof(double));
  }

  covMatrix = new double*[tot_dim];
  for(int i=0;i<tot_dim;i++){
    covMatrix[i] = new double[tot_dim];
    memset(covMatrix[i],0,tot_dim*sizeof(double));
  }
  
  work1 = new double*[h_size-1];
  for(int i=0;i<h_size-1;i++){
    work1[i] = new double[h_size-1];
    memset(work1[i],0,(h_size-1)*sizeof(double));
  }
  
  work2 = new double*[h_size-1];
  for(int i=0;i<h_size-1;i++){
    work2[i] = new double[h_size-1];
    memset(work2[i],0,(h_size-1)*sizeof(double));
  }  
}


void IQLBEstimator::deallocate(){

  hap_pair_map.clear();

  int h_size = hap_freq_table.size();
  int col = geno_freq_table.size();
  clear_matrix(rst,h_size-1,col);
  delete[] freq;
  delete[] genofreq;

  int par1_one_size = par1_one.size();
  if(par1_one_size>0){
    map<int, int>::iterator pi=par1_one.begin();
    for(int i=0;i<par1_one_size;i++){
      int dim_1 = missing_trio_list[pi->first]->dim;
      for(int j=0;j<dim_1;j++)
        delete[] result_par1_one[i][j];
      delete[] result_par1_one[i];
      pi++;
    }
    delete[] result_par1_one;
  }
  
  int par1_two_size = par1_two.size();
  if(par1_two_size>0){
    map<int, int>::iterator pi=par1_two.begin();
    for(int i=0;i<par1_two_size;i++){
      int dim_1 = missing_trio_list[pi->first]->dim;
      for(int j=0;j<dim_1;j++)
        delete[] result_par1_two[i][j];
      delete[] result_par1_two[i];
      pi++;
    }
    delete[] result_par1_two;
  }  

  int par2_one_size = par2_one.size();
  if(par2_one_size>0){
    map<int, int>::iterator pi=par2_one.begin();
    for(int i=0;i<par2_one_size;i++){
      int dim_1 = missing_trio_list[pi->first]->dim;
      for(int j=0;j<dim_1;j++)
        delete[] result_par2_one[i][j];
      delete[] result_par2_one[i];
      pi++;
    }
    delete[] result_par2_one;
  }  
  
  int par2_two_size = par2_two.size();
  if(par2_two_size>0){
    map<int, int>::iterator pi=par2_two.begin();
    for(int i=0;i<par2_two_size;i++){
      int dim_1 = missing_trio_list[pi->first]->dim;
      for(int j=0;j<dim_1;j++)
        delete[] result_par2_two[i][j];
      delete[] result_par2_two[i];
      pi++;
    }
    delete[] result_par2_two;
  }  

  int off_one_size = off_one.size();
  if(off_one_size>0){
    map<int, int>::iterator pi=off_one.begin();
    for(int i=0;i<off_one_size;i++){
      int dim_1 = missing_trio_list[pi->first]->dim;
      for(int j=0;j<dim_1;j++)
        delete[] result_off_one[i][j];
      delete[] result_off_one[i];
      pi++;
    }
    delete[] result_off_one;
  }
  
  int off_two_size = off_two.size();
  if(off_two_size>0){
    map<int, int>::iterator pi=off_two.begin();
    for(int i=0;i<off_two_size;i++){
      int dim_1 = missing_trio_list[pi->first]->dim;
      for(int j=0;j<dim_1;j++)
        delete[] result_off_two[i][j];
      delete[] result_off_two[i];
      pi++;
    }
    delete[] result_off_two;
  }  

  int parents_one_size = parents_one.size();
  if(parents_one_size>0){
    map<int, int>::iterator pi=parents_one.begin();
    for(int i=0;i<parents_one_size;i++){
      int dim_1 = missing_trio_list[pi->first]->dim;
      for(int j=0;j<dim_1;j++){
        for(int k=0;k<h_size;k++)
          delete[] result_parents_one[i][j][k];
        delete[] result_parents_one[i][j];
      }
      delete[] result_parents_one[i];
      pi++;
    }
    delete[] result_parents_one;
  }
  
  int parents_two_size = parents_two.size();
  if(parents_two_size>0){
    map<int, int>::iterator pi=parents_two.begin();
    for(int i=0;i<parents_two_size;i++){
      int dim_1 = missing_trio_list[pi->first]->dim;
      for(int j=0;j<dim_1;j++){
        for(int k=0;k<col;k++)
          delete[] result_parents_two[i][j][k];
        delete[] result_parents_two[i][j];
      }  
      delete[] result_parents_two[i];
      pi++;
    }
    delete[] result_parents_two;
  }  

  for(int i=0;i<tot_dim;i++){
    delete[] cholent[i];
    delete[] chol[i];
    delete[] cholaug[i];
    delete[] covMatrix[i];
  }
  delete[] cholent;
  delete[] chol;
  delete[] cholaug;
  delete[] covMatrix;

  clear_matrix(denominator,h_size-1,h_size-1);
  delete[] numerator;
  
  for(int i=0;i<h_size-1;i++){
    delete[] work1[i];
    delete[] work2[i];
  }
  delete[] work1;
  delete[] work2; 
}


// a single Newton-Raphson iteration
void IQLBEstimator::NR_Iteration_new(){

  int h_size = hap_freq_table.size();
  int col = geno_freq_table.size();

  for(int i=0;i<h_size-1;i++)
    memset(rst[i],0,col*sizeof(double));  // record derivative for each pair of haplotypes

  int loc = 0;
  map<int, double>::iterator iter=hap_freq_table.begin();
  while(iter!=hap_freq_table.end()){
    freq[loc] = iter->second;
    loc++;
    iter++;
  }
  
  loc = 0;
  map<pair<int, int>, double>::iterator giter=geno_freq_table.begin();
  while(giter!=geno_freq_table.end()){
    genofreq[loc] = giter->second;
    loc++;
    giter++;
  }

  int h_last=0;
  iter = hap_freq_table.end();
  iter--;
  h_last=iter->first;

  giter=geno_freq_table.begin();
  while(giter!=geno_freq_table.end()){
    pair<int, int> conf = giter->first;
    loc = 0;
    iter = hap_freq_table.begin();
    while(loc<h_size-1){
      int h=iter->first;

      if(conf.first == conf.second && conf.first == h)
        rst[loc][hap_pair_map[conf]] = 2*hap_freq_table[h];
      else if(conf.first == conf.second && conf.first == h_last)
        rst[loc][hap_pair_map[conf]] = -2*hap_freq_table[h_last];
      else if(conf.first < conf.second && conf.first == h && conf.second != h_last)
        rst[loc][hap_pair_map[conf]] = 2*hap_freq_table[conf.second];
      else if(conf.first < conf.second && conf.second == h)
        rst[loc][hap_pair_map[conf]] = 2*hap_freq_table[conf.first];
      else if(conf.first == h && conf.second == h_last)
        rst[loc][hap_pair_map[conf]] = 2*(hap_freq_table[h_last]-hap_freq_table[h]);
      else if(conf.first != h && conf.first != h_last && conf.second == h_last)
        rst[loc][hap_pair_map[conf]] = -2*hap_freq_table[conf.first];

      loc++;
      iter++;
    }
    giter++;
  }

  for(int i=0;i<geno_trio_list.size();i++){
    if(geno_trio_list[i]->get_num()>0 && geno_trio_list[i]->info==1 && missing_trio_list[geno_missing_trio_map[i]]->dim>0)
      geno_trio_list[i]->expected_hap_num(hap_freq_table,freq,geno_freq_table,hap_pair_map,count,missing_trio_list[geno_missing_trio_map[i]]->sites);
  }

  for(int i=0;i<missing_trio_list.size();i++){
    if(missing_trio_list[i]->get_num()>0 && missing_trio_list[i]->dim>0){
      missing_trio_list[i]->set_weight(hap_freq_table,geno_freq_table,hap_pair_map,count);
      missing_trio_list[i]->compute_centmatrix(hap_freq_table,freq,geno_freq_table);
      missing_trio_list[i]->compute_derivative(hap_freq_table,geno_freq_table,hap_pair_map,count,rst,hst);
      missing_trio_list[i]->expected_hap_num(hap_freq_table,freq,geno_freq_table,hap_pair_map,count);
    }
  }

  int par1_one_size = par1_one.size();
  if(par1_one_size>0){
    map<int, int>::iterator pi=par1_one.begin();
    for(int i=0;i<par1_one_size;i++){
      pi->second = i;
      int mis = pi->first;
      int dim_1 = missing_trio_list[mis]->dim;
      for(int j=0;j<dim_1;j++)
        memset(result_par1_one[i][j],0,h_size*sizeof(double));
      missing_trio_list[mis]->cond_par1_one(result_par1_one[i],hap_freq_table,geno_freq_table,hap_map);
      pi++;
    }
  }

  int par1_two_size = par1_two.size();
  if(par1_two_size>0){
    map<int, int>::iterator pi=par1_two.begin();
    for(int i=0;i<par1_two_size;i++){
      pi->second = i;
      int mis = pi->first;
      int dim_1 = missing_trio_list[mis]->dim;
      for(int j=0;j<dim_1;j++)
        memset(result_par1_two[i][j],0,col*sizeof(double));
      missing_trio_list[mis]->cond_par1_two(result_par1_two[i],hap_freq_table,geno_freq_table,hap_pair_map);  
      pi++;
    }
  }
  
  int par2_one_size = par2_one.size();
  if(par2_one_size>0){
    map<int, int>::iterator pi=par2_one.begin();
    for(int i=0;i<par2_one_size;i++){
      pi->second = i;
      int mis = pi->first;
      int dim_1 = missing_trio_list[mis]->dim;
      for(int j=0;j<dim_1;j++)
        memset(result_par2_one[i][j],0,h_size*sizeof(double));
      missing_trio_list[mis]->cond_par2_one(result_par2_one[i],hap_freq_table,geno_freq_table,hap_map);  
      pi++;
    }
  }
  
  int par2_two_size = par2_two.size();
  if(par2_two_size>0){
    map<int, int>::iterator pi=par2_two.begin();
    for(int i=0;i<par2_two_size;i++){
      pi->second = i;
      int mis = pi->first;
      int dim_1 = missing_trio_list[mis]->dim;
      for(int j=0;j<dim_1;j++)
        memset(result_par2_two[i][j],0,col*sizeof(double));
      missing_trio_list[mis]->cond_par2_two(result_par2_two[i],hap_freq_table,geno_freq_table,hap_pair_map);  
      pi++;
    }
  }  

  int off_one_size = off_one.size();
  if(off_one_size>0){
    map<int, int>::iterator pi=off_one.begin();
    for(int i=0;i<off_one_size;i++){
      pi->second = i;
      int mis = pi->first;
      int dim_1 = missing_trio_list[mis]->dim;
      for(int j=0;j<dim_1;j++)
        memset(result_off_one[i][j],0,h_size*sizeof(double));
      missing_trio_list[mis]->cond_off_one(result_off_one[i],hap_freq_table,geno_freq_table,hap_map);  
      pi++;
    }
  }
  
  int off_two_size = off_two.size();
  if(off_two_size>0){
    map<int, int>::iterator pi=off_two.begin();
    for(int i=0;i<off_two_size;i++){
      pi->second = i;
      int mis = pi->first;
      int dim_1 = missing_trio_list[mis]->dim;
      for(int j=0;j<dim_1;j++)
        memset(result_off_two[i][j],0,col*sizeof(double));
      missing_trio_list[mis]->cond_off_two(result_off_two[i],hap_freq_table,geno_freq_table,hap_pair_map);
      pi++;
    }
  }  

  int parents_one_size = parents_one.size();
  if(parents_one_size>0){
    map<int, int>::iterator pi=parents_one.begin();
    for(int i=0;i<parents_one_size;i++){
      pi->second = i;
      int mis = pi->first;
      int dim_1 = missing_trio_list[mis]->dim;
      for(int j=0;j<dim_1;j++){
        for(int k=0;k<h_size;k++)
          memset(result_parents_one[i][j][k],0,h_size*sizeof(double));  
      }
      missing_trio_list[mis]->cond_parents_one(result_parents_one[i],hap_freq_table,geno_freq_table,hap_map);
      pi++;
    }
  }
  
  int parents_two_size = parents_two.size();
  if(parents_two_size>0){
    map<int, int>::iterator pi=parents_two.begin();
    for(int i=0;i<parents_two_size;i++){
      pi->second = i;
      int mis = pi->first;
      int dim_1 = missing_trio_list[mis]->dim;
      for(int j=0;j<dim_1;j++){
        for(int k=0;k<col;k++)
          memset(result_parents_two[i][j][k],0,col*sizeof(double));
      }
      missing_trio_list[mis]->cond_parents_two(result_parents_two[i],hap_freq_table,geno_freq_table,hap_pair_map);
      pi++;
    }
  }

  for(int i=0;i<h_size-1;i++)
    memset(denominator[i],0,(h_size-1)*sizeof(double));  // F_p^T Omega^(-1) F_p

  memset(numerator,0,(h_size-1)*sizeof(double));  // F_p^T Omega^(-1) (Z-\mu)

  for(int fam=0;fam<fam_list.size();fam++){
    if(fam_list[fam]->N_typed>0){
      int trio_cat;
      int mis_trio_cat;
      int dim=0;
      for(int i=0;i<fam_list[fam]->N;i++){
        trio_cat = fam_list[fam]->trio_cat[i];
        if(geno_trio_list[trio_cat]->info==1){
          mis_trio_cat = fam_list[fam]->mis_trio_cat[i];
          dim+=missing_trio_list[mis_trio_cat]->dim;
        }
      }

      if(dim<=0)
        continue;

      for(int i=0;i<tot_dim;i++){
        memset(cholent[i],0,h_size*sizeof(double));
        memset(chol[i],0,tot_dim*sizeof(double));
        memset(cholaug[i],0,h_size*sizeof(double));
        memset(covMatrix[i],0,tot_dim*sizeof(double));
      }

      // populate covariance matrix
      int nb1 = 0;
      int nb2 = 0;
      int ind_1, ind_2, ind_1f, ind_1m, ind_2f, ind_2m;
      int par_1, par_2;  // how many parents is typed
      int dim_1, dim_2;
      int trio_cat_1, trio_cat_2;
      int mis_trio_cat_1, mis_trio_cat_2;
      double IBD_2, IBD_1, kin;
      map<pair<int, int>, double>::iterator iterk;
      map<pair<int, int>, double>::iterator iter1;
      map<pair<int, int>, double>::iterator iter2;
      double value;

      for(int i=0;i<fam_list[fam]->N;i++){
        trio_cat_1 = fam_list[fam]->trio_cat[i];
        mis_trio_cat_1 = fam_list[fam]->mis_trio_cat[i];
        dim_1 = missing_trio_list[mis_trio_cat_1]->dim;
        if(dim_1 > 0 && geno_trio_list[trio_cat_1]->info==1){
          nb2 = nb1;
          ind_1 = fam_list[fam]->fam_member[i][0];
          ind_2 = ind_1;
          iterk = (fam_list[fam]->kinship_coeff).find(pair<int, int>(ind_1,ind_2));
          if(iterk == (fam_list[fam]->kinship_coeff).end()){
            printf("No IBD coefficient for individual %d from family %d. Please check...\n\n", ind_1,fam_list[fam]->fam_id);
            exit(1);
          }

          for(int j=0;j<dim_1;j++)
            for(int k=0;k<dim_1;k++)
              covMatrix[nb1+j][nb2+k] = missing_trio_list[mis_trio_cat_1]->centmatrix[j][k];

          nb2+=dim_1;

          par_1 = 0;
          par_1+=missing_trio_list[mis_trio_cat_1]->typed_par1;
          par_1+=missing_trio_list[mis_trio_cat_1]->typed_par2;
          ind_1f = fam_list[fam]->fam_member[i][1];
          ind_1m = fam_list[fam]->fam_member[i][2];

          for(int j=i+1;j<fam_list[fam]->N;j++){
            trio_cat_2 = fam_list[fam]->trio_cat[j];
            mis_trio_cat_2 = fam_list[fam]->mis_trio_cat[j];
            dim_2 = missing_trio_list[mis_trio_cat_2]->dim;
            if(dim_2 > 0 && geno_trio_list[trio_cat_2]->info==1){
              ind_2 = fam_list[fam]->fam_member[j][0];
              iterk = (fam_list[fam]->kinship_coeff).find(pair<int, int>(ind_1,ind_2));
              if(iterk == (fam_list[fam]->kinship_coeff).end()){
                printf("No IBD coefficient between individual %d and individual %d from family %d. Please check...\n\n", ind_1,ind_2,fam_list[fam]->fam_id);
                exit(1);
              }

              if(iterk->second>1e-8){ // i is related to j
                par_2 = 0;
                par_2+=missing_trio_list[mis_trio_cat_2]->typed_par1;
                par_2+=missing_trio_list[mis_trio_cat_2]->typed_par2;
                ind_2f = fam_list[fam]->fam_member[j][1];
                ind_2m = fam_list[fam]->fam_member[j][2];

                kin = iterk->second;
                iter2 = (fam_list[fam]->IBD2_coeff).find(pair<int, int>(ind_1,ind_2));
                iter1 = (fam_list[fam]->IBD1_coeff).find(pair<int, int>(ind_1,ind_2));
                IBD_2 = iter2->second;
                IBD_1 = iter1->second;

                if(par_1==0 && par_2==0){
                  if(IBD_1>1e-8){
                    for(int k=0;k<h_size-1;k++)
                      memset(work1[k],0,(h_size-1)*sizeof(double));
                      
                    for(int k=0;k<dim_1;k++)
                      for(int l=0;l<dim_2;l++)
                        for(int m=0;m<h_size;m++)
                          work1[k][l]+=result_off_one[off_one[mis_trio_cat_1]][k][m]*result_off_one[off_one[mis_trio_cat_2]][l][m]*freq[m];
                          
                    for(int k=0;k<dim_1;k++)
                      for(int l=0;l<dim_2;l++)
                        covMatrix[nb1+k][nb2+l]+=IBD_1*work1[k][l];        
                  }
                  if(IBD_2>1e-8){
                    for(int k=0;k<h_size-1;k++)
                      memset(work2[k],0,(h_size-1)*sizeof(double));
                      
                    for(int k=0;k<dim_1;k++)
                      for(int l=0;l<dim_2;l++)
                        for(int m=0;m<col;m++)
                          work2[k][l]+=result_off_two[off_two[mis_trio_cat_1]][k][m]*result_off_two[off_two[mis_trio_cat_2]][l][m]*genofreq[m];
                          
                    for(int k=0;k<dim_1;k++)
                      for(int l=0;l<dim_2;l++)
                        covMatrix[nb1+k][nb2+l]+=IBD_2*work2[k][l];      
                  }
                }
                else if(par_1==0 && (par_2==1 || par_2==2)){
                  int order = fam_list[fam]->mis_par_order[j];
                  int ind_2par1 = ind_2f;
                  int ind_2par2 = ind_2m;
                  if(order==1){
                    ind_2par1 = ind_2m;
                    ind_2par2 = ind_2f;
                  }  
                  int diff = 0;
                  int index = relation(fam,ind_1,ind_2,ind_2par1,ind_2par2,kin,IBD_2,IBD_1,diff);
                  if(index==4){
                    for(int k=0;k<h_size-1;k++)
                      memset(work2[k],0,(h_size-1)*sizeof(double));

                    for(int k=0;k<dim_1;k++)
                      for(int l=0;l<dim_2;l++)
                        for(int m=0;m<col;m++)
                          work2[k][l]+=result_off_two[off_two[mis_trio_cat_1]][k][m]*result_par1_two[par1_two[mis_trio_cat_2]][l][m]*genofreq[m];
                          
                    for(int k=0;k<dim_1;k++)
                      for(int l=0;l<dim_2;l++)
                        covMatrix[nb1+k][nb2+l]+=work2[k][l];
                  }
                  else if(index==5){
                    for(int k=0;k<h_size-1;k++)
                      memset(work2[k],0,(h_size-1)*sizeof(double));
                      
                    for(int k=0;k<dim_1;k++)
                      for(int l=0;l<dim_2;l++)
                        for(int m=0;m<col;m++)
                          work2[k][l]+=result_off_two[off_two[mis_trio_cat_1]][k][m]*result_par2_two[par2_two[mis_trio_cat_2]][l][m]*genofreq[m];
                          
                    for(int k=0;k<dim_1;k++)
                      for(int l=0;l<dim_2;l++)
                        covMatrix[nb1+k][nb2+l]+=work2[k][l];  
                  }
                  else if(index==1){
                    iter2 = (fam_list[fam]->IBD2_coeff).find(pair<int, int>(ind_1,ind_2par1));
                    iter1 = (fam_list[fam]->IBD1_coeff).find(pair<int, int>(ind_1,ind_2par1));
                    IBD_2 = iter2->second;
                    IBD_1 = iter1->second;
                    
                    if(IBD_1>1e-8){
                      for(int k=0;k<h_size-1;k++)
                        memset(work1[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<h_size;m++)
                            work1[k][l]+=result_off_one[off_one[mis_trio_cat_1]][k][m]*result_par1_one[par1_one[mis_trio_cat_2]][l][m]*freq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=IBD_1*work1[k][l];        
                    }
                    if(IBD_2>1e-8){
                      for(int k=0;k<h_size-1;k++)
                        memset(work2[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<col;m++)
                            work2[k][l]+=result_off_two[off_two[mis_trio_cat_1]][k][m]*result_par1_two[par1_two[mis_trio_cat_2]][l][m]*genofreq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=IBD_2*work2[k][l];      
                    }                  
                  }
                  else if(index==2){
                    iter2 = (fam_list[fam]->IBD2_coeff).find(pair<int, int>(ind_1,ind_2par2));
                    iter1 = (fam_list[fam]->IBD1_coeff).find(pair<int, int>(ind_1,ind_2par2));
                    IBD_2 = iter2->second;
                    IBD_1 = iter1->second;
                    
                    if(IBD_1>1e-8){
                      for(int k=0;k<h_size-1;k++)
                        memset(work1[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<h_size;m++)
                            work1[k][l]+=result_off_one[off_one[mis_trio_cat_1]][k][m]*result_par2_one[par2_one[mis_trio_cat_2]][l][m]*freq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=IBD_1*work1[k][l];        
                    }
                    if(IBD_2>1e-8){
                      for(int k=0;k<h_size-1;k++)
                        memset(work2[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<col;m++)
                            work2[k][l]+=result_off_two[off_two[mis_trio_cat_1]][k][m]*result_par2_two[par2_two[mis_trio_cat_2]][l][m]*genofreq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=IBD_2*work2[k][l];      
                    }   
                  }
                  else if(index==3){
                    if(IBD_1>1e-8){
                      for(int k=0;k<h_size-1;k++)
                        memset(work1[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<h_size;m++)
                            work1[k][l]+=result_off_one[off_one[mis_trio_cat_1]][k][m]*result_off_one[off_one[mis_trio_cat_2]][l][m]*freq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=IBD_1*work1[k][l];        
                    }
                    if(IBD_2>1e-8){
                      for(int k=0;k<h_size-1;k++)
                        memset(work2[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<col;m++)
                            work2[k][l]+=result_off_two[off_two[mis_trio_cat_1]][k][m]*result_off_two[off_two[mis_trio_cat_2]][l][m]*genofreq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=IBD_2*work2[k][l];      
                    }   
                  }
                  else if(index==6){  // can not be sibs
                    if(diff==1)
                      printf("please check relationships\n");
                    
                    iter1 = (fam_list[fam]->IBD1_coeff).find(pair<int, int>(ind_1,ind_2par1));
                    IBD_1 = iter1->second;
                    
                    if(IBD_1>1e-8){
                      for(int k=0;k<h_size-1;k++)
                        memset(work1[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<h_size;m++)
                            work1[k][l]+=result_off_one[off_one[mis_trio_cat_1]][k][m]*result_par1_one[par1_one[mis_trio_cat_2]][l][m]*freq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=IBD_1*work1[k][l];        
                    }  
                    
                    iter1 = (fam_list[fam]->IBD1_coeff).find(pair<int, int>(ind_1,ind_2par2));
                    IBD_1 = iter1->second;
                    
                    if(IBD_1>1e-8){
                      for(int k=0;k<h_size-1;k++)
                        memset(work1[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<h_size;m++)
                            work1[k][l]+=result_off_one[off_one[mis_trio_cat_1]][k][m]*result_par2_one[par2_one[mis_trio_cat_2]][l][m]*freq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=IBD_1*work1[k][l];        
                    }               
                  } 
                  // index 7 & 8                        
                }
                else if((par_1==1 || par_1==2) && par_2==0){
                  int order = fam_list[fam]->mis_par_order[i];
                  int ind_1par1 = ind_1f;
                  int ind_1par2 = ind_1m;
                  if(order==1){
                    ind_1par1 = ind_1m;
                    ind_1par2 = ind_1f;
                  }  
                  int diff = 0;
                  int index = relation(fam,ind_2,ind_1,ind_1par1,ind_1par2,kin,IBD_2,IBD_1,diff);
                  if(index==4){
                    for(int k=0;k<h_size-1;k++)
                      memset(work2[k],0,(h_size-1)*sizeof(double));
                      
                    for(int k=0;k<dim_1;k++)
                      for(int l=0;l<dim_2;l++)
                        for(int m=0;m<col;m++)
                          work2[k][l]+=result_par1_two[par1_two[mis_trio_cat_1]][k][m]*result_off_two[off_two[mis_trio_cat_2]][l][m]*genofreq[m];
                          
                    for(int k=0;k<dim_1;k++)
                      for(int l=0;l<dim_2;l++)
                        covMatrix[nb1+k][nb2+l]+=work2[k][l];           
                  }
                  else if(index==5){
                    for(int k=0;k<h_size-1;k++)
                      memset(work2[k],0,(h_size-1)*sizeof(double));
                      
                    for(int k=0;k<dim_1;k++)
                      for(int l=0;l<dim_2;l++)
                        for(int m=0;m<col;m++)
                          work2[k][l]+=result_par2_two[par2_two[mis_trio_cat_1]][k][m]*result_off_two[off_two[mis_trio_cat_2]][l][m]*genofreq[m];
                          
                    for(int k=0;k<dim_1;k++)
                      for(int l=0;l<dim_2;l++)
                        covMatrix[nb1+k][nb2+l]+=work2[k][l];  
                  }
                  else if(index==1){
                    iter2 = (fam_list[fam]->IBD2_coeff).find(pair<int, int>(ind_1par1,ind_2));
                    iter1 = (fam_list[fam]->IBD1_coeff).find(pair<int, int>(ind_1par1,ind_2));
                    IBD_2 = iter2->second;
                    IBD_1 = iter1->second;
                    
                    if(IBD_1>1e-8){
                      for(int k=0;k<h_size-1;k++)
                        memset(work1[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<h_size;m++)
                            work1[k][l]+=result_par1_one[par1_one[mis_trio_cat_1]][k][m]*result_off_one[off_one[mis_trio_cat_2]][l][m]*freq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=IBD_1*work1[k][l];        
                    }
                    if(IBD_2>1e-8){
                      for(int k=0;k<h_size-1;k++)
                        memset(work2[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<col;m++)
                            work2[k][l]+=result_par1_two[par1_two[mis_trio_cat_1]][k][m]*result_off_two[off_two[mis_trio_cat_2]][l][m]*genofreq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=IBD_2*work2[k][l];      
                    }                  
                  }
                  else if(index==2){
                    iter2 = (fam_list[fam]->IBD2_coeff).find(pair<int, int>(ind_1par2,ind_2));
                    iter1 = (fam_list[fam]->IBD1_coeff).find(pair<int, int>(ind_1par2,ind_2));
                    IBD_2 = iter2->second;
                    IBD_1 = iter1->second;
                    
                    if(IBD_1>1e-8){
                      for(int k=0;k<h_size-1;k++)
                        memset(work1[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<h_size;m++)
                            work1[k][l]+=result_par2_one[par2_one[mis_trio_cat_1]][k][m]*result_off_one[off_one[mis_trio_cat_2]][l][m]*freq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=IBD_1*work1[k][l];        
                    }
                    if(IBD_2>1e-8){
                      for(int k=0;k<h_size-1;k++)
                        memset(work2[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<col;m++)
                            work2[k][l]+=result_par2_two[par2_two[mis_trio_cat_1]][k][m]*result_off_two[off_two[mis_trio_cat_2]][l][m]*genofreq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=IBD_2*work2[k][l];      
                    }   
                  }
                  else if(index==3){
                    if(IBD_1>1e-8){
                      for(int k=0;k<h_size-1;k++)
                        memset(work1[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<h_size;m++)
                            work1[k][l]+=result_off_one[off_one[mis_trio_cat_1]][k][m]*result_off_one[off_one[mis_trio_cat_2]][l][m]*freq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=IBD_1*work1[k][l];        
                    }
                    if(IBD_2>1e-8){
                      for(int k=0;k<h_size-1;k++)
                        memset(work2[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<col;m++)
                            work2[k][l]+=result_off_two[off_two[mis_trio_cat_1]][k][m]*result_off_two[off_two[mis_trio_cat_2]][l][m]*genofreq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=IBD_2*work2[k][l];      
                    }   
                  }
                  else if(index==6){  // can not be sibs
                    if(diff==1)
                      printf("please check relationships\n");
                    
                    iter1 = (fam_list[fam]->IBD1_coeff).find(pair<int, int>(ind_1par1,ind_2));
                    IBD_1 = iter1->second;
                    
                    if(IBD_1>1e-8){
                      for(int k=0;k<h_size-1;k++)
                        memset(work1[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<h_size;m++)
                            work1[k][l]+=result_par1_one[par1_one[mis_trio_cat_1]][k][m]*result_off_one[off_one[mis_trio_cat_2]][l][m]*freq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=IBD_1*work1[k][l];        
                    }  
                    
                    iter1 = (fam_list[fam]->IBD1_coeff).find(pair<int, int>(ind_1par2,ind_2));
                    IBD_1 = iter1->second;
                    
                    if(IBD_1>1e-8){
                      for(int k=0;k<h_size-1;k++)
                        memset(work1[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<h_size;m++)
                            work1[k][l]+=result_par2_one[par2_one[mis_trio_cat_1]][k][m]*result_off_one[off_one[mis_trio_cat_2]][l][m]*freq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=IBD_1*work1[k][l];        
                    }               
                  } 
                  // index 7 & 8                        
                }
                else if((par_1==1 || par_1==2) && (par_2==1 || par_2==2)){
                  int order_1 = fam_list[fam]->mis_par_order[i];
                  int ind_1par1 = ind_1f;
                  int ind_1par2 = ind_1m;
                  if(order_1==1){
                    ind_1par1 = ind_1m;
                    ind_1par2 = ind_1f;
                  }
                  int order_2 = fam_list[fam]->mis_par_order[j];
                  int ind_2par1 = ind_2f;
                  int ind_2par2 = ind_2m;
                  if(order_2==1){
                    ind_2par1 = ind_2m;
                    ind_2par2 = ind_2f;
                  }                    
                  int diff = 0;
                  int index = relation(fam,ind_1,ind_1par1,ind_1par2,ind_2,ind_2par1,ind_2par2,kin,IBD_2,IBD_1,diff);
                  if(index==13){
                    for(int k=0;k<h_size-1;k++)
                      memset(work2[k],0,(h_size-1)*sizeof(double));
                      
                    for(int k=0;k<dim_1;k++)
                      for(int l=0;l<dim_2;l++)
                        for(int m=0;m<col;m++)
                          work2[k][l]+=result_off_two[off_two[mis_trio_cat_1]][k][m]*result_par1_two[par1_two[mis_trio_cat_2]][l][m]*genofreq[m];
                          
                    for(int k=0;k<dim_1;k++)
                      for(int l=0;l<dim_2;l++)
                        covMatrix[nb1+k][nb2+l]+=work2[k][l];              
                  }
                  else if(index==14){
                    for(int k=0;k<h_size-1;k++)
                      memset(work2[k],0,(h_size-1)*sizeof(double));
                      
                    for(int k=0;k<dim_1;k++)
                      for(int l=0;l<dim_2;l++)
                        for(int m=0;m<col;m++)
                          work2[k][l]+=result_off_two[off_two[mis_trio_cat_1]][k][m]*result_par2_two[par2_two[mis_trio_cat_2]][l][m]*genofreq[m];
                          
                    for(int k=0;k<dim_1;k++)
                      for(int l=0;l<dim_2;l++)
                        covMatrix[nb1+k][nb2+l]+=work2[k][l];   
                  }
                  else if(index==15){
                    for(int k=0;k<h_size-1;k++)
                      memset(work2[k],0,(h_size-1)*sizeof(double));
                      
                    for(int k=0;k<dim_1;k++)
                      for(int l=0;l<dim_2;l++)
                        for(int m=0;m<col;m++)
                          work2[k][l]+=result_par1_two[par1_two[mis_trio_cat_1]][k][m]*result_off_two[off_two[mis_trio_cat_2]][l][m]*genofreq[m];
                          
                    for(int k=0;k<dim_1;k++)
                      for(int l=0;l<dim_2;l++)
                        covMatrix[nb1+k][nb2+l]+=work2[k][l];              
                  }
                  else if(index==16){
                    for(int k=0;k<h_size-1;k++)
                      memset(work2[k],0,(h_size-1)*sizeof(double));
                      
                    for(int k=0;k<dim_1;k++)
                      for(int l=0;l<dim_2;l++)
                        for(int m=0;m<col;m++)
                          work2[k][l]+=result_par2_two[par2_two[mis_trio_cat_1]][k][m]*result_off_two[off_two[mis_trio_cat_2]][l][m]*genofreq[m];
                          
                    for(int k=0;k<dim_1;k++)
                      for(int l=0;l<dim_2;l++)
                        covMatrix[nb1+k][nb2+l]+=work2[k][l];   
                  }
                  else if(index==17){
                    for(int k=0;k<h_size-1;k++)
                      memset(work2[k],0,(h_size-1)*sizeof(double));
                      
                    for(int k=0;k<dim_1;k++)
                      for(int l=0;l<dim_2;l++)
                        for(int m=0;m<col;m++)
                          for(int n=0;n<col;n++)
                            work2[k][l]+=result_parents_two[parents_two[mis_trio_cat_1]][k][m][n]*result_parents_two[parents_two[mis_trio_cat_2]][l][m][n]*genofreq[m]*genofreq[n];
                          
                    for(int k=0;k<dim_1;k++)
                      for(int l=0;l<dim_2;l++)
                        covMatrix[nb1+k][nb2+l]+=work2[k][l];     
                  }
                  else if(index==18){
                    if(ind_1par1==ind_2par1){
                      for(int k=0;k<h_size-1;k++)
                        memset(work2[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<col;m++)
                            work2[k][l]+=result_par1_two[par1_two[mis_trio_cat_1]][k][m]*result_par1_two[par1_two[mis_trio_cat_2]][l][m]*genofreq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=work2[k][l];                                          
                    }
                    else if(ind_1par1==ind_2par2){
                      for(int k=0;k<h_size-1;k++)
                        memset(work2[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<col;m++)
                            work2[k][l]+=result_par1_two[par1_two[mis_trio_cat_1]][k][m]*result_par2_two[par2_two[mis_trio_cat_2]][l][m]*genofreq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=work2[k][l];     
                    }
                    else if(ind_1par2==ind_2par1){
                      for(int k=0;k<h_size-1;k++)
                        memset(work2[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<col;m++)
                            work2[k][l]+=result_par2_two[par2_two[mis_trio_cat_1]][k][m]*result_par1_two[par1_two[mis_trio_cat_2]][l][m]*genofreq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=work2[k][l];      
                    }
                    else if(ind_1par2==ind_2par2){
                      for(int k=0;k<h_size-1;k++)
                        memset(work2[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<col;m++)
                            work2[k][l]+=result_par2_two[par2_two[mis_trio_cat_1]][k][m]*result_par2_two[par2_two[mis_trio_cat_2]][l][m]*genofreq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=work2[k][l];   
                    }     
                  }
                  else if(index==1){
                    iter2 = (fam_list[fam]->IBD2_coeff).find(pair<int, int>(ind_1par1,ind_2par1));
                    iter1 = (fam_list[fam]->IBD1_coeff).find(pair<int, int>(ind_1par1,ind_2par1));
                    IBD_2 = iter2->second;
                    IBD_1 = iter1->second;
                    
                    if(IBD_1>1e-8){
                      for(int k=0;k<h_size-1;k++)
                        memset(work1[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<h_size;m++)
                            work1[k][l]+=result_par1_one[par1_one[mis_trio_cat_1]][k][m]*result_par1_one[par1_one[mis_trio_cat_2]][l][m]*freq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=IBD_1*work1[k][l];        
                    }
                    if(IBD_2>1e-8){
                      for(int k=0;k<h_size-1;k++)
                        memset(work2[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<col;m++)
                            work2[k][l]+=result_par1_two[par1_two[mis_trio_cat_1]][k][m]*result_par1_two[par1_two[mis_trio_cat_2]][l][m]*genofreq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=IBD_2*work2[k][l];      
                    }   
                  }
                  else if(index==2){
                    iter2 = (fam_list[fam]->IBD2_coeff).find(pair<int, int>(ind_1par1,ind_2par2));
                    iter1 = (fam_list[fam]->IBD1_coeff).find(pair<int, int>(ind_1par1,ind_2par2));
                    IBD_2 = iter2->second;
                    IBD_1 = iter1->second;
                    
                    if(IBD_1>1e-8){
                      for(int k=0;k<h_size-1;k++)
                        memset(work1[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<h_size;m++)
                            work1[k][l]+=result_par1_one[par1_one[mis_trio_cat_1]][k][m]*result_par2_one[par2_one[mis_trio_cat_2]][l][m]*freq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=IBD_1*work1[k][l];        
                    }
                    if(IBD_2>1e-8){
                      for(int k=0;k<h_size-1;k++)
                        memset(work2[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<col;m++)
                            work2[k][l]+=result_par1_two[par1_two[mis_trio_cat_1]][k][m]*result_par2_two[par2_two[mis_trio_cat_2]][l][m]*genofreq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=IBD_2*work2[k][l];      
                    }   
                  }
                  else if(index==3){
                    iter2 = (fam_list[fam]->IBD2_coeff).find(pair<int, int>(ind_1par2,ind_2par1));
                    iter1 = (fam_list[fam]->IBD1_coeff).find(pair<int, int>(ind_1par2,ind_2par1));
                    IBD_2 = iter2->second;
                    IBD_1 = iter1->second;
                    
                    if(IBD_1>1e-8){
                      for(int k=0;k<h_size-1;k++)
                        memset(work1[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<h_size;m++)
                            work1[k][l]+=result_par2_one[par2_one[mis_trio_cat_1]][k][m]*result_par1_one[par1_one[mis_trio_cat_2]][l][m]*freq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=IBD_1*work1[k][l];        
                    }
                    if(IBD_2>1e-8){
                      for(int k=0;k<h_size-1;k++)
                        memset(work2[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<col;m++)
                            work2[k][l]+=result_par2_two[par2_two[mis_trio_cat_1]][k][m]*result_par1_two[par1_two[mis_trio_cat_2]][l][m]*genofreq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=IBD_2*work2[k][l];      
                    }   
                  }
                  else if(index==4){
                    iter2 = (fam_list[fam]->IBD2_coeff).find(pair<int, int>(ind_1par2,ind_2par2));
                    iter1 = (fam_list[fam]->IBD1_coeff).find(pair<int, int>(ind_1par2,ind_2par2));
                    IBD_2 = iter2->second;
                    IBD_1 = iter1->second;
                    
                    if(IBD_1>1e-8){
                      for(int k=0;k<h_size-1;k++)
                        memset(work1[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<h_size;m++)
                            work1[k][l]+=result_par2_one[par2_one[mis_trio_cat_1]][k][m]*result_par2_one[par2_one[mis_trio_cat_2]][l][m]*freq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=IBD_1*work1[k][l];        
                    }
                    if(IBD_2>1e-8){
                      for(int k=0;k<h_size-1;k++)
                        memset(work2[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<col;m++)
                            work2[k][l]+=result_par2_two[par2_two[mis_trio_cat_1]][k][m]*result_par2_two[par2_two[mis_trio_cat_2]][l][m]*genofreq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=IBD_2*work2[k][l];      
                    }   
                  }
                  else if(index==5){
                    iter2 = (fam_list[fam]->IBD2_coeff).find(pair<int, int>(ind_1par1,ind_2));
                    iter1 = (fam_list[fam]->IBD1_coeff).find(pair<int, int>(ind_1par1,ind_2));
                    IBD_2 = iter2->second;
                    IBD_1 = iter1->second;
                       
                    if(IBD_1>1e-8){
                      for(int k=0;k<h_size-1;k++)
                        memset(work1[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<h_size;m++)
                            work1[k][l]+=result_par1_one[par1_one[mis_trio_cat_1]][k][m]*result_off_one[off_one[mis_trio_cat_2]][l][m]*freq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=IBD_1*work1[k][l];        
                    }
                    if(IBD_2>1e-8){
                      for(int k=0;k<h_size-1;k++)
                        memset(work2[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<col;m++)
                            work2[k][l]+=result_par1_two[par1_two[mis_trio_cat_1]][k][m]*result_off_two[off_two[mis_trio_cat_2]][l][m]*genofreq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=IBD_2*work2[k][l];      
                    }     
                  }
                  else if(index==6){
                    iter2 = (fam_list[fam]->IBD2_coeff).find(pair<int, int>(ind_1par2,ind_2));
                    iter1 = (fam_list[fam]->IBD1_coeff).find(pair<int, int>(ind_1par2,ind_2));
                    IBD_2 = iter2->second;
                    IBD_1 = iter1->second;
                       
                    if(IBD_1>1e-8){
                      for(int k=0;k<h_size-1;k++)
                        memset(work1[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<h_size;m++)
                            work1[k][l]+=result_par2_one[par2_one[mis_trio_cat_1]][k][m]*result_off_one[off_one[mis_trio_cat_2]][l][m]*freq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=IBD_1*work1[k][l];        
                    }
                    if(IBD_2>1e-8){
                      for(int k=0;k<h_size-1;k++)
                        memset(work2[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<col;m++)
                            work2[k][l]+=result_par2_two[par2_two[mis_trio_cat_1]][k][m]*result_off_two[off_two[mis_trio_cat_2]][l][m]*genofreq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=IBD_2*work2[k][l];      
                    }     
                  }
                  else if(index==7){
                    iter2 = (fam_list[fam]->IBD2_coeff).find(pair<int, int>(ind_1,ind_2par1));
                    iter1 = (fam_list[fam]->IBD1_coeff).find(pair<int, int>(ind_1,ind_2par1));
                    IBD_2 = iter2->second;
                    IBD_1 = iter1->second;
                       
                    if(IBD_1>1e-8){
                      for(int k=0;k<h_size-1;k++)
                        memset(work1[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<h_size;m++)
                            work1[k][l]+=result_off_one[off_one[mis_trio_cat_1]][k][m]*result_par1_one[par1_one[mis_trio_cat_2]][l][m]*freq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=IBD_1*work1[k][l];        
                    }
                    if(IBD_2>1e-8){
                      for(int k=0;k<h_size-1;k++)
                        memset(work2[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<col;m++)
                            work2[k][l]+=result_off_two[off_two[mis_trio_cat_1]][k][m]*result_par1_two[par1_two[mis_trio_cat_2]][l][m]*genofreq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=IBD_2*work2[k][l];      
                    }     
                  }
                  else if(index==8){
                    iter2 = (fam_list[fam]->IBD2_coeff).find(pair<int, int>(ind_1,ind_2par2));
                    iter1 = (fam_list[fam]->IBD1_coeff).find(pair<int, int>(ind_1,ind_2par2));
                    IBD_2 = iter2->second;
                    IBD_1 = iter1->second;
                       
                    if(IBD_1>1e-8){
                      for(int k=0;k<h_size-1;k++)
                        memset(work1[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<h_size;m++)
                            work1[k][l]+=result_off_one[off_one[mis_trio_cat_1]][k][m]*result_par2_one[par2_one[mis_trio_cat_2]][l][m]*freq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=IBD_1*work1[k][l];        
                    }
                    if(IBD_2>1e-8){
                      for(int k=0;k<h_size-1;k++)
                        memset(work2[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<col;m++)
                            work2[k][l]+=result_off_two[off_two[mis_trio_cat_1]][k][m]*result_par2_two[par2_two[mis_trio_cat_2]][l][m]*genofreq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=IBD_2*work2[k][l];      
                    }     
                  }     
                  else if(index==9){
                    if(diff==1){  // j and if are sibs
                      for(int k=0;k<h_size-1;k++)
                        memset(work2[k],0,(h_size-1)*sizeof(double));
                      
                      loc=0;
                      giter = geno_freq_table.begin();
                      while(giter!=geno_freq_table.end()){
                        int h1 = (giter->first).first;
                        int h2 = (giter->first).second;
                        if(h1==h2){
                          for(int k=0;k<dim_1;k++)
                            for(int l=0;l<dim_2;l++)
                              work2[k][l]+=result_par1_two[par1_two[mis_trio_cat_1]][k][loc]*result_parents_one[parents_one[mis_trio_cat_2]][l][hap_map[h1]][hap_map[h2]]*genofreq[loc];
                        }
                        else{
                          for(int k=0;k<dim_1;k++)
                            for(int l=0;l<dim_2;l++)
                              work2[k][l]+=result_par1_two[par1_two[mis_trio_cat_1]][k][loc]*(result_parents_one[parents_one[mis_trio_cat_2]][l][hap_map[h1]][hap_map[h2]]+result_parents_one[parents_one[mis_trio_cat_2]][l][hap_map[h2]][hap_map[h1]])*genofreq[loc]/2.0;
                        }
                        loc++;
                        giter++;
                      }

                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=work2[k][l];
                    }
                    else{
                      iter1 = (fam_list[fam]->IBD1_coeff).find(pair<int, int>(ind_1par1,ind_2par1));
                      IBD_1 = iter1->second;
                    
                      if(IBD_1>1e-8){
                        for(int k=0;k<h_size-1;k++)
                          memset(work1[k],0,(h_size-1)*sizeof(double));
                      
                        for(int k=0;k<dim_1;k++)
                          for(int l=0;l<dim_2;l++)
                            for(int m=0;m<h_size;m++)
                              work1[k][l]+=result_par1_one[par1_one[mis_trio_cat_1]][k][m]*result_par1_one[par1_one[mis_trio_cat_2]][l][m]*freq[m];
                          
                        for(int k=0;k<dim_1;k++)
                          for(int l=0;l<dim_2;l++)
                            covMatrix[nb1+k][nb2+l]+=IBD_1*work1[k][l];
                      }
                    
                      iter1 = (fam_list[fam]->IBD1_coeff).find(pair<int, int>(ind_1par1,ind_2par2));
                      IBD_1 = iter1->second;
                    
                      if(IBD_1>1e-8){
                        for(int k=0;k<h_size-1;k++)
                          memset(work1[k],0,(h_size-1)*sizeof(double));
                      
                        for(int k=0;k<dim_1;k++)
                          for(int l=0;l<dim_2;l++)
                            for(int m=0;m<h_size;m++)
                              work1[k][l]+=result_par1_one[par1_one[mis_trio_cat_1]][k][m]*result_par2_one[par2_one[mis_trio_cat_2]][l][m]*freq[m];
                          
                        for(int k=0;k<dim_1;k++)
                          for(int l=0;l<dim_2;l++)
                            covMatrix[nb1+k][nb2+l]+=IBD_1*work1[k][l];
                      }
                    }                          
                  }
                  else if(index==10){
                    if(diff==1){  // j and im are sibs
                      for(int k=0;k<h_size-1;k++)
                        memset(work2[k],0,(h_size-1)*sizeof(double));
                      
                      loc=0;
                      giter = geno_freq_table.begin();
                      while(giter!=geno_freq_table.end()){
                        int h1 = (giter->first).first;
                        int h2 = (giter->first).second;
                        if(h1==h2){
                          for(int k=0;k<dim_1;k++)
                            for(int l=0;l<dim_2;l++)
                              work2[k][l]+=result_par2_two[par2_two[mis_trio_cat_1]][k][loc]*result_parents_one[parents_one[mis_trio_cat_2]][l][hap_map[h1]][hap_map[h2]]*genofreq[loc];
                        }
                        else{
                          for(int k=0;k<dim_1;k++)
                            for(int l=0;l<dim_2;l++)
                              work2[k][l]+=result_par2_two[par2_two[mis_trio_cat_1]][k][loc]*(result_parents_one[parents_one[mis_trio_cat_2]][l][hap_map[h1]][hap_map[h2]]+result_parents_one[parents_one[mis_trio_cat_2]][l][hap_map[h2]][hap_map[h1]])*genofreq[loc]/2.0;
                        }
                        loc++;
                        giter++;
                      }

                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=work2[k][l];
                    }
                    else{
                      iter1 = (fam_list[fam]->IBD1_coeff).find(pair<int, int>(ind_1par2,ind_2par1));
                      IBD_1 = iter1->second;
                    
                      if(IBD_1>1e-8){
                        for(int k=0;k<h_size-1;k++)
                          memset(work1[k],0,(h_size-1)*sizeof(double));
                      
                        for(int k=0;k<dim_1;k++)
                          for(int l=0;l<dim_2;l++)
                            for(int m=0;m<h_size;m++)
                              work1[k][l]+=result_par2_one[par2_one[mis_trio_cat_1]][k][m]*result_par1_one[par1_one[mis_trio_cat_2]][l][m]*freq[m];
                          
                        for(int k=0;k<dim_1;k++)
                          for(int l=0;l<dim_2;l++)
                            covMatrix[nb1+k][nb2+l]+=IBD_1*work1[k][l];
                      }
                    
                      iter1 = (fam_list[fam]->IBD1_coeff).find(pair<int, int>(ind_1par2,ind_2par2));
                      IBD_1 = iter1->second;
                    
                      if(IBD_1>1e-8){
                        for(int k=0;k<h_size-1;k++)
                          memset(work1[k],0,(h_size-1)*sizeof(double));
                      
                        for(int k=0;k<dim_1;k++)
                          for(int l=0;l<dim_2;l++)
                            for(int m=0;m<h_size;m++)
                              work1[k][l]+=result_par2_one[par2_one[mis_trio_cat_1]][k][m]*result_par2_one[par2_one[mis_trio_cat_2]][l][m]*freq[m];
                          
                        for(int k=0;k<dim_1;k++)
                          for(int l=0;l<dim_2;l++)
                            covMatrix[nb1+k][nb2+l]+=IBD_1*work1[k][l];
                      }
                    }                          
                  }
                  else if(index==11){
                    if(diff==1){  // i and jf are sibs
                      for(int k=0;k<h_size-1;k++)
                        memset(work2[k],0,(h_size-1)*sizeof(double));
                      
                      loc=0;
                      giter = geno_freq_table.begin();
                      while(giter!=geno_freq_table.end()){
                        int h1 = (giter->first).first;
                        int h2 = (giter->first).second;
                        if(h1==h2){
                          for(int k=0;k<dim_1;k++)
                            for(int l=0;l<dim_2;l++)
                              work2[k][l]+=result_parents_one[parents_one[mis_trio_cat_1]][k][hap_map[h1]][hap_map[h2]]*result_par1_two[par1_two[mis_trio_cat_2]][l][loc]*genofreq[loc];
                        }
                        else{
                          for(int k=0;k<dim_1;k++)
                            for(int l=0;l<dim_2;l++)
                              work2[k][l]+=(result_parents_one[parents_one[mis_trio_cat_1]][k][hap_map[h1]][hap_map[h2]]+result_parents_one[parents_one[mis_trio_cat_1]][k][hap_map[h2]][hap_map[h1]])*result_par1_two[par1_two[mis_trio_cat_2]][l][loc]*genofreq[loc]/2.0;
                        }
                        loc++;
                        giter++;
                      }

                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=work2[k][l];
                    }
                    else{
                      iter1 = (fam_list[fam]->IBD1_coeff).find(pair<int, int>(ind_1par1,ind_2par1));
                      IBD_1 = iter1->second;
                    
                      if(IBD_1>1e-8){
                        for(int k=0;k<h_size-1;k++)
                          memset(work1[k],0,(h_size-1)*sizeof(double));
                      
                        for(int k=0;k<dim_1;k++)
                          for(int l=0;l<dim_2;l++)
                            for(int m=0;m<h_size;m++)
                              work1[k][l]+=result_par1_one[par1_one[mis_trio_cat_1]][k][m]*result_par1_one[par1_one[mis_trio_cat_2]][l][m]*freq[m];
                          
                        for(int k=0;k<dim_1;k++)
                          for(int l=0;l<dim_2;l++)
                            covMatrix[nb1+k][nb2+l]+=IBD_1*work1[k][l];
                      }
                    
                      iter1 = (fam_list[fam]->IBD1_coeff).find(pair<int, int>(ind_1par2,ind_2par1));
                      IBD_1 = iter1->second;
                    
                      if(IBD_1>1e-8){
                        for(int k=0;k<h_size-1;k++)
                          memset(work1[k],0,(h_size-1)*sizeof(double));
                      
                        for(int k=0;k<dim_1;k++)
                          for(int l=0;l<dim_2;l++)
                            for(int m=0;m<h_size;m++)
                              work1[k][l]+=result_par2_one[par2_one[mis_trio_cat_1]][k][m]*result_par1_one[par1_one[mis_trio_cat_2]][l][m]*freq[m];
                          
                        for(int k=0;k<dim_1;k++)
                          for(int l=0;l<dim_2;l++)
                            covMatrix[nb1+k][nb2+l]+=IBD_1*work1[k][l];
                      }
                    }                          
                  }
                  else if(index==12){
                    if(diff==1){  // i and jm are sibs
                      for(int k=0;k<h_size-1;k++)
                        memset(work2[k],0,(h_size-1)*sizeof(double));
                      
                      loc=0;
                      giter = geno_freq_table.begin();
                      while(giter!=geno_freq_table.end()){
                        int h1 = (giter->first).first;
                        int h2 = (giter->first).second;
                        if(h1==h2){
                          for(int k=0;k<dim_1;k++)
                            for(int l=0;l<dim_2;l++)
                              work2[k][l]+=result_parents_one[parents_one[mis_trio_cat_1]][k][hap_map[h1]][hap_map[h2]]*result_par2_two[par2_two[mis_trio_cat_2]][l][loc]*genofreq[loc];
                        }
                        else{
                          for(int k=0;k<dim_1;k++)
                            for(int l=0;l<dim_2;l++)
                              work2[k][l]+=(result_parents_one[parents_one[mis_trio_cat_1]][k][hap_map[h1]][hap_map[h2]]+result_parents_one[parents_one[mis_trio_cat_1]][k][hap_map[h2]][hap_map[h1]])*result_par2_two[par2_two[mis_trio_cat_2]][l][loc]*genofreq[loc]/2.0;
                        }
                        loc++;
                        giter++;
                      }

                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=work2[k][l];
                    }
                    else{
                      iter1 = (fam_list[fam]->IBD1_coeff).find(pair<int, int>(ind_1par1,ind_2par2));
                      IBD_1 = iter1->second;
                    
                      if(IBD_1>1e-8){
                        for(int k=0;k<h_size-1;k++)
                          memset(work1[k],0,(h_size-1)*sizeof(double));
                      
                        for(int k=0;k<dim_1;k++)
                          for(int l=0;l<dim_2;l++)
                            for(int m=0;m<h_size;m++)
                              work1[k][l]+=result_par1_one[par1_one[mis_trio_cat_1]][k][m]*result_par2_one[par2_one[mis_trio_cat_2]][l][m]*freq[m];
                          
                        for(int k=0;k<dim_1;k++)
                          for(int l=0;l<dim_2;l++)
                            covMatrix[nb1+k][nb2+l]+=IBD_1*work1[k][l];
                      }
                    
                      iter1 = (fam_list[fam]->IBD1_coeff).find(pair<int, int>(ind_1par2,ind_2par2));
                      IBD_1 = iter1->second;
                    
                      if(IBD_1>1e-8){
                        for(int k=0;k<h_size-1;k++)
                          memset(work1[k],0,(h_size-1)*sizeof(double));
                      
                        for(int k=0;k<dim_1;k++)
                          for(int l=0;l<dim_2;l++)
                            for(int m=0;m<h_size;m++)
                              work1[k][l]+=result_par2_one[par2_one[mis_trio_cat_1]][k][m]*result_par2_one[par2_one[mis_trio_cat_2]][l][m]*freq[m];
                          
                        for(int k=0;k<dim_1;k++)
                          for(int l=0;l<dim_2;l++)
                            covMatrix[nb1+k][nb2+l]+=IBD_1*work1[k][l];
                      }
                    }                          
                  }   
                  // index 19 & 20
                }     
                for(int k=0;k<dim_1;k++)
                  for(int l=0;l<dim_2;l++)
                    covMatrix[nb2+l][nb1+k] = covMatrix[nb1+k][nb2+l];
              }  // i, j related
              nb2+=dim_2;
            }  // j is typed
          }  // complete j

          const vector<double> hap_num = geno_trio_list[trio_cat_1]->get_hap_num();
          for(int j=0;j<dim_1;j++){
            for(int k=0;k<h_size-1;k++)
              cholent[nb1+j][k] = missing_trio_list[mis_trio_cat_1]->derivative[j][k];
            cholent[nb1+j][h_size-1] = hap_num[j];
          }
          nb1+=dim_1;
        }  // i is typed
      }  // complete i

      if(cholesky(covMatrix,dim,cholent,h_size,chol,cholaug,0) != 1){
        //printf("cholesky decomposition of the cov matrix failed for family %d. \n", fam_list[fam]->fam_id);
        flag_pd = 0;
        //break;
      }

      if(flag_pd==1){
      double **s = self_colmultiply(cholaug,dim,1,h_size-1,1,h_size-1);
      double **t = self_colmultiply(cholaug,dim,1,h_size-1,h_size,h_size);

      for(int i=0;i<h_size-1;i++)
        for(int j=0;j<h_size-1;j++)
          denominator[i][j]+=s[i][j];

      for(int i=0;i<h_size-1;i++)
        numerator[i]+=t[i][0];

      clear_matrix(s,h_size-1,h_size-1);
      clear_matrix(t,h_size-1,1);
      }
      else{
      double var;
      double &determ = var;
      double **v = inverse(covMatrix,dim,determ);
      double **s = quadratic_product(cholent,dim,h_size,v,dim,dim,cholent,dim,h_size);

      for(int i=0;i<h_size-1;i++)
        for(int j=0;j<h_size-1;j++)
          denominator[i][j]+=s[i][j];

      for(int i=0;i<h_size-1;i++)
        numerator[i]+=s[i][h_size-1];

      clear_matrix(s,h_size,h_size);
      clear_matrix(v,dim,dim);
      }
    }
  }

  double **cholentmain = new double*[h_size-1];
  for(int i=0;i<h_size-1;i++){
    cholentmain[i] = new double[h_size];
    memset(cholentmain[i],0,h_size*sizeof(double));
  }

  double **cholmain = new double*[h_size-1];
  for(int i=0;i<h_size-1;i++){
    cholmain[i] = new double[h_size-1];
    memset(cholmain[i],0,(h_size-1)*sizeof(double));
  }

  double **cholaugmain = new double*[h_size-1];
  for(int i=0;i<h_size-1;i++){
    cholaugmain[i] = new double[h_size];
    memset(cholaugmain[i],0,h_size*sizeof(double));
  }

  for(int i=0;i<h_size-1;i++){
    cholentmain[i][i] = 1;
    cholentmain[i][h_size-1] = numerator[i];
  }

  if(cholesky(denominator,h_size-1,cholentmain,h_size,cholmain,cholaugmain,0) != 1){
    //printf("cholesky decomposition of the denominator matrix failed.\n");
    flag_pd = 0;
    //exit(1);
  }

  if(flag_pd==1){
  double **s = self_colmultiply(cholaugmain,h_size-1,1,h_size-1,h_size,h_size);
  iter=hap_freq_table.begin();
  double n=0;
  for(int i=0;i<h_size-1;i++){
    iter->second+=s[i][0];
    n+=iter->second;
    iter++;
  }
  iter->second = 1-n;

  clear_matrix(s,h_size-1,1);
  }
  else{
  double var;
  double &determ = var;
  double **v = inverse(denominator,h_size-1,determ);
  double **w = new double*[h_size-1];
  for(int i=0;i<h_size-1;i++){
    w[i] = new double[1];
    w[i][0] = numerator[i];
  }
  double **s = multiply(v,h_size-1,h_size-1,w,h_size-1,1);
  iter=hap_freq_table.begin();
  double n=0;
  for(int i=0;i<h_size-1;i++){
    iter->second+=s[i][0];
    n+=iter->second;
    iter++;
  }
  iter->second = 1-n;

  clear_matrix(s,h_size-1,1);
  clear_matrix(v,h_size-1,h_size-1);
  clear_matrix(w,h_size-1,1);
  }
  clear_matrix(cholentmain,h_size-1,h_size);
  clear_matrix(cholmain,h_size-1,h_size-1);
  clear_matrix(cholaugmain,h_size-1,h_size);

  update_geno_freq();
}


void IQLBEstimator::Score_test_new(map<int, double> & hap_weight){
     
  int h_size = hap_freq_table.size();
  int col = geno_freq_table.size();

  for(int i=0;i<h_size-1;i++)
    memset(rst[i],0,col*sizeof(double));  // record derivative for each pair of haplotypes

  int loc = 0;
  map<int, double>::iterator iter=hap_freq_table.begin();
  while(iter!=hap_freq_table.end()){
    freq[loc] = iter->second;
    loc++;
    iter++;
  }
  
  loc = 0;
  map<pair<int, int>, double>::iterator giter=geno_freq_table.begin();
  while(giter!=geno_freq_table.end()){
    genofreq[loc] = giter->second;
    loc++;
    giter++;
  }

  vector<double> h_weight;
  iter=hap_weight.begin();
  while(iter!=hap_weight.end()){
    h_weight.push_back(iter->second);
    iter++;
  }

  int h_last=0;
  iter = hap_freq_table.end();
  iter--;
  h_last=iter->first;

  giter=geno_freq_table.begin();
  while(giter!=geno_freq_table.end()){
    pair<int, int> conf = giter->first;
    loc = 0;
    iter = hap_freq_table.begin();
    while(loc<h_size-1){
      int h=iter->first;

      if(conf.first == conf.second && conf.first == h)
        rst[loc][hap_pair_map[conf]] = 2*hap_freq_table[h];
      else if(conf.first == conf.second && conf.first == h_last)
        rst[loc][hap_pair_map[conf]] = -2*hap_freq_table[h_last];
      else if(conf.first < conf.second && conf.first == h && conf.second != h_last)
        rst[loc][hap_pair_map[conf]] = 2*hap_freq_table[conf.second];
      else if(conf.first < conf.second && conf.second == h)
        rst[loc][hap_pair_map[conf]] = 2*hap_freq_table[conf.first];
      else if(conf.first == h && conf.second == h_last)
        rst[loc][hap_pair_map[conf]] = 2*(hap_freq_table[h_last]-hap_freq_table[h]);
      else if(conf.first != h && conf.first != h_last && conf.second == h_last)
        rst[loc][hap_pair_map[conf]] = -2*hap_freq_table[conf.first];

      loc++;
      iter++;
    }
    giter++;
  }

  for(int i=0;i<geno_trio_list.size();i++){
    if(geno_trio_list[i]->get_num()>0 && geno_trio_list[i]->info==1 && missing_trio_list[geno_missing_trio_map[i]]->dim>0)
      geno_trio_list[i]->expected_hap_num(hap_freq_table,freq,geno_freq_table,hap_pair_map,count,missing_trio_list[geno_missing_trio_map[i]]->sites);
  }

  for(int i=0;i<missing_trio_list.size();i++){
    if(missing_trio_list[i]->get_num()>0 && missing_trio_list[i]->dim>0){
      missing_trio_list[i]->set_weight(hap_freq_table,geno_freq_table,hap_pair_map,count);
      missing_trio_list[i]->compute_centmatrix(hap_freq_table,freq,geno_freq_table);
      missing_trio_list[i]->compute_derivative(hap_freq_table,geno_freq_table,hap_pair_map,count,rst,hst);
      missing_trio_list[i]->expected_hap_num(hap_freq_table,freq,geno_freq_table,hap_pair_map,count);
      if(missing_trio_list[i]->typed_par1==1 || missing_trio_list[i]->typed_par2==1)
        missing_trio_list[i]->compute_deriv1to3(hap_freq_table,geno_freq_table,hap_pair_map,hap_weight,count);
      else if(missing_trio_list[i]->typed_par1==0 && missing_trio_list[i]->typed_par2==0)
        missing_trio_list[i]->compute_deriv_r(freq,geno_freq_table,hap_pair_map,hap_weight,h_weight,count);
    }
  }

  int par1_one_size = par1_one.size();
  if(par1_one_size>0){
    map<int, int>::iterator pi=par1_one.begin();
    for(int i=0;i<par1_one_size;i++){
      pi->second = i;
      int mis = pi->first;
      int dim_1 = missing_trio_list[mis]->dim;
      for(int j=0;j<dim_1;j++)
        memset(result_par1_one[i][j],0,h_size*sizeof(double));
      missing_trio_list[mis]->cond_par1_one(result_par1_one[i],hap_freq_table,geno_freq_table,hap_map);  
      pi++;
    }
  }
  
  int par1_two_size = par1_two.size();
  if(par1_two_size>0){
    map<int, int>::iterator pi=par1_two.begin();
    for(int i=0;i<par1_two_size;i++){
      pi->second = i;
      int mis = pi->first;
      int dim_1 = missing_trio_list[mis]->dim;
      for(int j=0;j<dim_1;j++)
        memset(result_par1_two[i][j],0,col*sizeof(double));
      missing_trio_list[mis]->cond_par1_two(result_par1_two[i],hap_freq_table,geno_freq_table,hap_pair_map);  
      pi++;
    }
  }
  
  int par2_one_size = par2_one.size();
  if(par2_one_size>0){
    map<int, int>::iterator pi=par2_one.begin();
    for(int i=0;i<par2_one_size;i++){
      pi->second = i;
      int mis = pi->first;
      int dim_1 = missing_trio_list[mis]->dim;
      for(int j=0;j<dim_1;j++)
        memset(result_par2_one[i][j],0,h_size*sizeof(double));
      missing_trio_list[mis]->cond_par2_one(result_par2_one[i],hap_freq_table,geno_freq_table,hap_map);  
      pi++;
    }
  }
  
  int par2_two_size = par2_two.size();
  if(par2_two_size>0){
    map<int, int>::iterator pi=par2_two.begin();
    for(int i=0;i<par2_two_size;i++){
      pi->second = i;
      int mis = pi->first;
      int dim_1 = missing_trio_list[mis]->dim;
      for(int j=0;j<dim_1;j++)
        memset(result_par2_two[i][j],0,col*sizeof(double));
      missing_trio_list[mis]->cond_par2_two(result_par2_two[i],hap_freq_table,geno_freq_table,hap_pair_map);  
      pi++;
    }
  }  

  int off_one_size = off_one.size();
  if(off_one_size>0){
    map<int, int>::iterator pi=off_one.begin();
    for(int i=0;i<off_one_size;i++){
      pi->second = i;
      int mis = pi->first;
      int dim_1 = missing_trio_list[mis]->dim;
      for(int j=0;j<dim_1;j++)
        memset(result_off_one[i][j],0,h_size*sizeof(double));
      missing_trio_list[mis]->cond_off_one(result_off_one[i],hap_freq_table,geno_freq_table,hap_map);  
      pi++;
    }
  }
  
  int off_two_size = off_two.size();
  if(off_two_size>0){
    map<int, int>::iterator pi=off_two.begin();
    for(int i=0;i<off_two_size;i++){
      pi->second = i;
      int mis = pi->first;
      int dim_1 = missing_trio_list[mis]->dim;
      for(int j=0;j<dim_1;j++)
        memset(result_off_two[i][j],0,col*sizeof(double));
      missing_trio_list[mis]->cond_off_two(result_off_two[i],hap_freq_table,geno_freq_table,hap_pair_map);
      pi++;
    }
  }  

  int parents_one_size = parents_one.size();
  if(parents_one_size>0){
    map<int, int>::iterator pi=parents_one.begin();
    for(int i=0;i<parents_one_size;i++){
      pi->second = i;
      int mis = pi->first;
      int dim_1 = missing_trio_list[mis]->dim;
      for(int j=0;j<dim_1;j++){
        for(int k=0;k<h_size;k++)
          memset(result_parents_one[i][j][k],0,h_size*sizeof(double));  
      }
      missing_trio_list[mis]->cond_parents_one(result_parents_one[i],hap_freq_table,geno_freq_table,hap_map);
      pi++;
    }
  }
  
  int parents_two_size = parents_two.size();
  if(parents_two_size>0){
    map<int, int>::iterator pi=parents_two.begin();
    for(int i=0;i<parents_two_size;i++){
      pi->second = i;
      int mis = pi->first;
      int dim_1 = missing_trio_list[mis]->dim;
      for(int j=0;j<dim_1;j++){
        for(int k=0;k<col;k++)
          memset(result_parents_two[i][j][k],0,col*sizeof(double));
      }
      missing_trio_list[mis]->cond_parents_two(result_parents_two[i],hap_freq_table,geno_freq_table,hap_pair_map);
      pi++;
    }
  }
  
  for(int i=0;i<h_size-1;i++)
    memset(denominator[i],0,(h_size-1)*sizeof(double));  // F_p^T Omega^(-1) F_p

  double A1 = 0;  // F_r^T Omega^(-1) F_r

  double **A2 = new double*[h_size-1];
  for(int i=0;i<h_size-1;i++){
    A2[i] = new double[1];
    memset(A2[i],0,sizeof(double));
  }  // F_p^T Omega^(-1) F_r

  double numerator_1 = 0;  // F_r^T Omega^(-1) (Z-\mu)

  for(int i=0;i<tot_dim;i++){
    delete[] cholent[i];
    cholent[i] = new double[h_size+1];
    memset(cholent[i],0,(h_size+1)*sizeof(double));
  }

  for(int i=0;i<tot_dim;i++){
    delete[] cholaug[i];
    cholaug[i] = new double[h_size+1];
    memset(cholaug[i],0,(h_size+1)*sizeof(double));
  }

  for(int fam=0;fam<fam_list.size();fam++){
    if(fam_list[fam]->N_typed>0){
      int trio_cat;
      int mis_trio_cat;
      int dim=0;
      for(int i=0;i<fam_list[fam]->N;i++){
        trio_cat = fam_list[fam]->trio_cat[i];
        if(geno_trio_list[trio_cat]->info==1){
          mis_trio_cat = fam_list[fam]->mis_trio_cat[i];
          dim+=missing_trio_list[mis_trio_cat]->dim;
        }
      }

      if(dim<=0)
        continue;

      for(int i=0;i<tot_dim;i++){
        memset(cholent[i],0,(h_size+1)*sizeof(double));
        memset(chol[i],0,tot_dim*sizeof(double));
        memset(cholaug[i],0,(h_size+1)*sizeof(double));
        memset(covMatrix[i],0,tot_dim*sizeof(double));
      }

      // populate covariance matrix
      int nb1 = 0;
      int nb2 = 0;
      int ind_1, ind_2, ind_1f, ind_1m, ind_2f, ind_2m;
      int par_1, par_2;  // how many parents is typed
      int dim_1, dim_2;
      int trio_cat_1, trio_cat_2;
      int mis_trio_cat_1, mis_trio_cat_2;
      double IBD_2, IBD_1, kin;
      map<pair<int, int>, double>::iterator iterk;
      map<pair<int, int>, double>::iterator iter1;
      map<pair<int, int>, double>::iterator iter2;
      double value;
      double test_weight;
      double test_trio_weight1, test_trio_weight2, test_trio_weight3;
      
      for(int i=0;i<fam_list[fam]->N;i++){
        trio_cat_1 = fam_list[fam]->trio_cat[i];
        mis_trio_cat_1 = fam_list[fam]->mis_trio_cat[i];
        dim_1 = missing_trio_list[mis_trio_cat_1]->dim;
        if(dim_1 > 0 && geno_trio_list[trio_cat_1]->info==1){
          nb2 = nb1;
          ind_1 = fam_list[fam]->fam_member[i][0];
          ind_2 = ind_1;
          iterk = (fam_list[fam]->kinship_coeff).find(pair<int, int>(ind_1,ind_2));
          if(iterk == (fam_list[fam]->kinship_coeff).end()){
            printf("No IBD coefficient for individual %d from family %d. Please check...\n\n", ind_1,fam_list[fam]->fam_id);
            exit(1);
          }

          for(int j=0;j<dim_1;j++)
            for(int k=0;k<dim_1;k++)
              covMatrix[nb1+j][nb2+k] = missing_trio_list[mis_trio_cat_1]->centmatrix[j][k];

          nb2+=dim_1;

          par_1 = 0;
          par_1+=missing_trio_list[mis_trio_cat_1]->typed_par1;
          par_1+=missing_trio_list[mis_trio_cat_1]->typed_par2;
          ind_1f = fam_list[fam]->fam_member[i][1];
          ind_1m = fam_list[fam]->fam_member[i][2];

          for(int j=i+1;j<fam_list[fam]->N;j++){
            trio_cat_2 = fam_list[fam]->trio_cat[j];
            mis_trio_cat_2 = fam_list[fam]->mis_trio_cat[j];
            dim_2 = missing_trio_list[mis_trio_cat_2]->dim;
            if(dim_2 > 0 && geno_trio_list[trio_cat_2]->info==1){
              ind_2 = fam_list[fam]->fam_member[j][0];
              iterk = (fam_list[fam]->kinship_coeff).find(pair<int, int>(ind_1,ind_2));
              if(iterk == (fam_list[fam]->kinship_coeff).end()){
                printf("No IBD coefficient between individual %d and individual %d from family %d. Please check...\n\n", ind_1,ind_2,fam_list[fam]->fam_id);
                exit(1);
              }

              if(iterk->second>1e-8){ // i is related to j
                par_2 = 0;
                par_2+=missing_trio_list[mis_trio_cat_2]->typed_par1;
                par_2+=missing_trio_list[mis_trio_cat_2]->typed_par2;
                ind_2f = fam_list[fam]->fam_member[j][1];
                ind_2m = fam_list[fam]->fam_member[j][2];

                kin = iterk->second;
                iter2 = (fam_list[fam]->IBD2_coeff).find(pair<int, int>(ind_1,ind_2));
                iter1 = (fam_list[fam]->IBD1_coeff).find(pair<int, int>(ind_1,ind_2));
                IBD_2 = iter2->second;
                IBD_1 = iter1->second;

                if(par_1==0 && par_2==0){
                  if(IBD_1>1e-8){
                    for(int k=0;k<h_size-1;k++)
                      memset(work1[k],0,(h_size-1)*sizeof(double));
                      
                    for(int k=0;k<dim_1;k++)
                      for(int l=0;l<dim_2;l++)
                        for(int m=0;m<h_size;m++)
                          work1[k][l]+=result_off_one[off_one[mis_trio_cat_1]][k][m]*result_off_one[off_one[mis_trio_cat_2]][l][m]*freq[m];
                          
                    for(int k=0;k<dim_1;k++)
                      for(int l=0;l<dim_2;l++)
                        covMatrix[nb1+k][nb2+l]+=IBD_1*work1[k][l];        
                  }
                  if(IBD_2>1e-8){
                    for(int k=0;k<h_size-1;k++)
                      memset(work2[k],0,(h_size-1)*sizeof(double));
                      
                    for(int k=0;k<dim_1;k++)
                      for(int l=0;l<dim_2;l++)
                        for(int m=0;m<col;m++)
                          work2[k][l]+=result_off_two[off_two[mis_trio_cat_1]][k][m]*result_off_two[off_two[mis_trio_cat_2]][l][m]*genofreq[m];
                          
                    for(int k=0;k<dim_1;k++)
                      for(int l=0;l<dim_2;l++)
                        covMatrix[nb1+k][nb2+l]+=IBD_2*work2[k][l];      
                  }
                }
                else if(par_1==0 && (par_2==1 || par_2==2)){
                  int order = fam_list[fam]->mis_par_order[j];
                  int ind_2par1 = ind_2f;
                  int ind_2par2 = ind_2m;
                  if(order==1){
                    ind_2par1 = ind_2m;
                    ind_2par2 = ind_2f;
                  }  
                  int diff = 0;
                  int index = relation(fam,ind_1,ind_2,ind_2par1,ind_2par2,kin,IBD_2,IBD_1,diff);
                  if(index==4){
                    for(int k=0;k<h_size-1;k++)
                      memset(work2[k],0,(h_size-1)*sizeof(double));
                      
                    for(int k=0;k<dim_1;k++)
                      for(int l=0;l<dim_2;l++)
                        for(int m=0;m<col;m++)
                          work2[k][l]+=result_off_two[off_two[mis_trio_cat_1]][k][m]*result_par1_two[par1_two[mis_trio_cat_2]][l][m]*genofreq[m];
                          
                    for(int k=0;k<dim_1;k++)
                      for(int l=0;l<dim_2;l++)
                        covMatrix[nb1+k][nb2+l]+=work2[k][l];           
                  }
                  else if(index==5){
                    for(int k=0;k<h_size-1;k++)
                      memset(work2[k],0,(h_size-1)*sizeof(double));
                      
                    for(int k=0;k<dim_1;k++)
                      for(int l=0;l<dim_2;l++)
                        for(int m=0;m<col;m++)
                          work2[k][l]+=result_off_two[off_two[mis_trio_cat_1]][k][m]*result_par2_two[par2_two[mis_trio_cat_2]][l][m]*genofreq[m];
                          
                    for(int k=0;k<dim_1;k++)
                      for(int l=0;l<dim_2;l++)
                        covMatrix[nb1+k][nb2+l]+=work2[k][l];  
                  }
                  else if(index==1){
                    iter2 = (fam_list[fam]->IBD2_coeff).find(pair<int, int>(ind_1,ind_2par1));
                    iter1 = (fam_list[fam]->IBD1_coeff).find(pair<int, int>(ind_1,ind_2par1));
                    IBD_2 = iter2->second;
                    IBD_1 = iter1->second;
                    
                    if(IBD_1>1e-8){
                      for(int k=0;k<h_size-1;k++)
                        memset(work1[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<h_size;m++)
                            work1[k][l]+=result_off_one[off_one[mis_trio_cat_1]][k][m]*result_par1_one[par1_one[mis_trio_cat_2]][l][m]*freq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=IBD_1*work1[k][l];        
                    }
                    if(IBD_2>1e-8){
                      for(int k=0;k<h_size-1;k++)
                        memset(work2[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<col;m++)
                            work2[k][l]+=result_off_two[off_two[mis_trio_cat_1]][k][m]*result_par1_two[par1_two[mis_trio_cat_2]][l][m]*genofreq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=IBD_2*work2[k][l];      
                    }                  
                  }
                  else if(index==2){
                    iter2 = (fam_list[fam]->IBD2_coeff).find(pair<int, int>(ind_1,ind_2par2));
                    iter1 = (fam_list[fam]->IBD1_coeff).find(pair<int, int>(ind_1,ind_2par2));
                    IBD_2 = iter2->second;
                    IBD_1 = iter1->second;
                    
                    if(IBD_1>1e-8){
                      for(int k=0;k<h_size-1;k++)
                        memset(work1[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<h_size;m++)
                            work1[k][l]+=result_off_one[off_one[mis_trio_cat_1]][k][m]*result_par2_one[par2_one[mis_trio_cat_2]][l][m]*freq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=IBD_1*work1[k][l];        
                    }
                    if(IBD_2>1e-8){
                      for(int k=0;k<h_size-1;k++)
                        memset(work2[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<col;m++)
                            work2[k][l]+=result_off_two[off_two[mis_trio_cat_1]][k][m]*result_par2_two[par2_two[mis_trio_cat_2]][l][m]*genofreq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=IBD_2*work2[k][l];      
                    }   
                  }
                  else if(index==3){
                    if(IBD_1>1e-8){
                      for(int k=0;k<h_size-1;k++)
                        memset(work1[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<h_size;m++)
                            work1[k][l]+=result_off_one[off_one[mis_trio_cat_1]][k][m]*result_off_one[off_one[mis_trio_cat_2]][l][m]*freq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=IBD_1*work1[k][l];        
                    }
                    if(IBD_2>1e-8){
                      for(int k=0;k<h_size-1;k++)
                        memset(work2[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<col;m++)
                            work2[k][l]+=result_off_two[off_two[mis_trio_cat_1]][k][m]*result_off_two[off_two[mis_trio_cat_2]][l][m]*genofreq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=IBD_2*work2[k][l];      
                    }   
                  }
                  else if(index==6){  // can not be sibs
                    if(diff==1)
                      printf("please check relationships\n");
                    
                    iter1 = (fam_list[fam]->IBD1_coeff).find(pair<int, int>(ind_1,ind_2par1));
                    IBD_1 = iter1->second;
                    
                    if(IBD_1>1e-8){
                      for(int k=0;k<h_size-1;k++)
                        memset(work1[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<h_size;m++)
                            work1[k][l]+=result_off_one[off_one[mis_trio_cat_1]][k][m]*result_par1_one[par1_one[mis_trio_cat_2]][l][m]*freq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=IBD_1*work1[k][l];        
                    }  
                    
                    iter1 = (fam_list[fam]->IBD1_coeff).find(pair<int, int>(ind_1,ind_2par2));
                    IBD_1 = iter1->second;
                    
                    if(IBD_1>1e-8){
                      for(int k=0;k<h_size-1;k++)
                        memset(work1[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<h_size;m++)
                            work1[k][l]+=result_off_one[off_one[mis_trio_cat_1]][k][m]*result_par2_one[par2_one[mis_trio_cat_2]][l][m]*freq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=IBD_1*work1[k][l];        
                    }               
                  } 
                  // index 7 & 8                        
                }
                else if((par_1==1 || par_1==2) && par_2==0){
                  int order = fam_list[fam]->mis_par_order[i];
                  int ind_1par1 = ind_1f;
                  int ind_1par2 = ind_1m;
                  if(order==1){
                    ind_1par1 = ind_1m;
                    ind_1par2 = ind_1f;
                  }  
                  int diff = 0;
                  int index = relation(fam,ind_2,ind_1,ind_1par1,ind_1par2,kin,IBD_2,IBD_1,diff);
                  if(index==4){
                    for(int k=0;k<h_size-1;k++)
                      memset(work2[k],0,(h_size-1)*sizeof(double));
                      
                    for(int k=0;k<dim_1;k++)
                      for(int l=0;l<dim_2;l++)
                        for(int m=0;m<col;m++)
                          work2[k][l]+=result_par1_two[par1_two[mis_trio_cat_1]][k][m]*result_off_two[off_two[mis_trio_cat_2]][l][m]*genofreq[m];
                          
                    for(int k=0;k<dim_1;k++)
                      for(int l=0;l<dim_2;l++)
                        covMatrix[nb1+k][nb2+l]+=work2[k][l];           
                  }
                  else if(index==5){
                    for(int k=0;k<h_size-1;k++)
                      memset(work2[k],0,(h_size-1)*sizeof(double));
                      
                    for(int k=0;k<dim_1;k++)
                      for(int l=0;l<dim_2;l++)
                        for(int m=0;m<col;m++)
                          work2[k][l]+=result_par2_two[par2_two[mis_trio_cat_1]][k][m]*result_off_two[off_two[mis_trio_cat_2]][l][m]*genofreq[m];
                          
                    for(int k=0;k<dim_1;k++)
                      for(int l=0;l<dim_2;l++)
                        covMatrix[nb1+k][nb2+l]+=work2[k][l];  
                  }
                  else if(index==1){
                    iter2 = (fam_list[fam]->IBD2_coeff).find(pair<int, int>(ind_1par1,ind_2));
                    iter1 = (fam_list[fam]->IBD1_coeff).find(pair<int, int>(ind_1par1,ind_2));
                    IBD_2 = iter2->second;
                    IBD_1 = iter1->second;
                    
                    if(IBD_1>1e-8){
                      for(int k=0;k<h_size-1;k++)
                        memset(work1[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<h_size;m++)
                            work1[k][l]+=result_par1_one[par1_one[mis_trio_cat_1]][k][m]*result_off_one[off_one[mis_trio_cat_2]][l][m]*freq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=IBD_1*work1[k][l];        
                    }
                    if(IBD_2>1e-8){
                      for(int k=0;k<h_size-1;k++)
                        memset(work2[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<col;m++)
                            work2[k][l]+=result_par1_two[par1_two[mis_trio_cat_1]][k][m]*result_off_two[off_two[mis_trio_cat_2]][l][m]*genofreq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=IBD_2*work2[k][l];      
                    }                  
                  }
                  else if(index==2){
                    iter2 = (fam_list[fam]->IBD2_coeff).find(pair<int, int>(ind_1par2,ind_2));
                    iter1 = (fam_list[fam]->IBD1_coeff).find(pair<int, int>(ind_1par2,ind_2));
                    IBD_2 = iter2->second;
                    IBD_1 = iter1->second;
                    
                    if(IBD_1>1e-8){
                      for(int k=0;k<h_size-1;k++)
                        memset(work1[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<h_size;m++)
                            work1[k][l]+=result_par2_one[par2_one[mis_trio_cat_1]][k][m]*result_off_one[off_one[mis_trio_cat_2]][l][m]*freq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=IBD_1*work1[k][l];        
                    }
                    if(IBD_2>1e-8){
                      for(int k=0;k<h_size-1;k++)
                        memset(work2[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<col;m++)
                            work2[k][l]+=result_par2_two[par2_two[mis_trio_cat_1]][k][m]*result_off_two[off_two[mis_trio_cat_2]][l][m]*genofreq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=IBD_2*work2[k][l];      
                    }   
                  }
                  else if(index==3){
                    if(IBD_1>1e-8){
                      for(int k=0;k<h_size-1;k++)
                        memset(work1[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<h_size;m++)
                            work1[k][l]+=result_off_one[off_one[mis_trio_cat_1]][k][m]*result_off_one[off_one[mis_trio_cat_2]][l][m]*freq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=IBD_1*work1[k][l];        
                    }
                    if(IBD_2>1e-8){
                      for(int k=0;k<h_size-1;k++)
                        memset(work2[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<col;m++)
                            work2[k][l]+=result_off_two[off_two[mis_trio_cat_1]][k][m]*result_off_two[off_two[mis_trio_cat_2]][l][m]*genofreq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=IBD_2*work2[k][l];      
                    }   
                  }
                  else if(index==6){  // can not be sibs
                    if(diff==1)
                      printf("please check relationships\n");
                    
                    iter1 = (fam_list[fam]->IBD1_coeff).find(pair<int, int>(ind_1par1,ind_2));
                    IBD_1 = iter1->second;
                    
                    if(IBD_1>1e-8){
                      for(int k=0;k<h_size-1;k++)
                        memset(work1[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<h_size;m++)
                            work1[k][l]+=result_par1_one[par1_one[mis_trio_cat_1]][k][m]*result_off_one[off_one[mis_trio_cat_2]][l][m]*freq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=IBD_1*work1[k][l];        
                    }  
                    
                    iter1 = (fam_list[fam]->IBD1_coeff).find(pair<int, int>(ind_1par2,ind_2));
                    IBD_1 = iter1->second;
                    
                    if(IBD_1>1e-8){
                      for(int k=0;k<h_size-1;k++)
                        memset(work1[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<h_size;m++)
                            work1[k][l]+=result_par2_one[par2_one[mis_trio_cat_1]][k][m]*result_off_one[off_one[mis_trio_cat_2]][l][m]*freq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=IBD_1*work1[k][l];        
                    }               
                  } 
                  // index 7 & 8                        
                }
                else if((par_1==1 || par_1==2) && (par_2==1 || par_2==2)){
                  int order_1 = fam_list[fam]->mis_par_order[i];
                  int ind_1par1 = ind_1f;
                  int ind_1par2 = ind_1m;
                  if(order_1==1){
                    ind_1par1 = ind_1m;
                    ind_1par2 = ind_1f;
                  }
                  int order_2 = fam_list[fam]->mis_par_order[j];
                  int ind_2par1 = ind_2f;
                  int ind_2par2 = ind_2m;
                  if(order_2==1){
                    ind_2par1 = ind_2m;
                    ind_2par2 = ind_2f;
                  }                    
                  int diff = 0;
                  int index = relation(fam,ind_1,ind_1par1,ind_1par2,ind_2,ind_2par1,ind_2par2,kin,IBD_2,IBD_1,diff);
                  if(index==13){
                    for(int k=0;k<h_size-1;k++)
                      memset(work2[k],0,(h_size-1)*sizeof(double));
                      
                    for(int k=0;k<dim_1;k++)
                      for(int l=0;l<dim_2;l++)
                        for(int m=0;m<col;m++)
                          work2[k][l]+=result_off_two[off_two[mis_trio_cat_1]][k][m]*result_par1_two[par1_two[mis_trio_cat_2]][l][m]*genofreq[m];
                          
                    for(int k=0;k<dim_1;k++)
                      for(int l=0;l<dim_2;l++)
                        covMatrix[nb1+k][nb2+l]+=work2[k][l];              
                  }
                  else if(index==14){
                    for(int k=0;k<h_size-1;k++)
                      memset(work2[k],0,(h_size-1)*sizeof(double));
                      
                    for(int k=0;k<dim_1;k++)
                      for(int l=0;l<dim_2;l++)
                        for(int m=0;m<col;m++)
                          work2[k][l]+=result_off_two[off_two[mis_trio_cat_1]][k][m]*result_par2_two[par2_two[mis_trio_cat_2]][l][m]*genofreq[m];
                          
                    for(int k=0;k<dim_1;k++)
                      for(int l=0;l<dim_2;l++)
                        covMatrix[nb1+k][nb2+l]+=work2[k][l];   
                  }
                  else if(index==15){
                    for(int k=0;k<h_size-1;k++)
                      memset(work2[k],0,(h_size-1)*sizeof(double));
                      
                    for(int k=0;k<dim_1;k++)
                      for(int l=0;l<dim_2;l++)
                        for(int m=0;m<col;m++)
                          work2[k][l]+=result_par1_two[par1_two[mis_trio_cat_1]][k][m]*result_off_two[off_two[mis_trio_cat_2]][l][m]*genofreq[m];
                          
                    for(int k=0;k<dim_1;k++)
                      for(int l=0;l<dim_2;l++)
                        covMatrix[nb1+k][nb2+l]+=work2[k][l];              
                  }
                  else if(index==16){
                    for(int k=0;k<h_size-1;k++)
                      memset(work2[k],0,(h_size-1)*sizeof(double));
                      
                    for(int k=0;k<dim_1;k++)
                      for(int l=0;l<dim_2;l++)
                        for(int m=0;m<col;m++)
                          work2[k][l]+=result_par2_two[par2_two[mis_trio_cat_1]][k][m]*result_off_two[off_two[mis_trio_cat_2]][l][m]*genofreq[m];
                          
                    for(int k=0;k<dim_1;k++)
                      for(int l=0;l<dim_2;l++)
                        covMatrix[nb1+k][nb2+l]+=work2[k][l];   
                  }
                  else if(index==17){
                    for(int k=0;k<h_size-1;k++)
                      memset(work2[k],0,(h_size-1)*sizeof(double));
                      
                    for(int k=0;k<dim_1;k++)
                      for(int l=0;l<dim_2;l++)
                        for(int m=0;m<col;m++)
                          for(int n=0;n<col;n++)
                            work2[k][l]+=result_parents_two[parents_two[mis_trio_cat_1]][k][m][n]*result_parents_two[parents_two[mis_trio_cat_2]][l][m][n]*genofreq[m]*genofreq[n];
                          
                    for(int k=0;k<dim_1;k++)
                      for(int l=0;l<dim_2;l++)
                        covMatrix[nb1+k][nb2+l]+=work2[k][l];     
                  }
                  else if(index==18){
                    if(ind_1par1==ind_2par1){
                      for(int k=0;k<h_size-1;k++)
                        memset(work2[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<col;m++)
                            work2[k][l]+=result_par1_two[par1_two[mis_trio_cat_1]][k][m]*result_par1_two[par1_two[mis_trio_cat_2]][l][m]*genofreq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=work2[k][l];                                          
                    }
                    else if(ind_1par1==ind_2par2){
                      for(int k=0;k<h_size-1;k++)
                        memset(work2[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<col;m++)
                            work2[k][l]+=result_par1_two[par1_two[mis_trio_cat_1]][k][m]*result_par2_two[par2_two[mis_trio_cat_2]][l][m]*genofreq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=work2[k][l];     
                    }
                    else if(ind_1par2==ind_2par1){
                      for(int k=0;k<h_size-1;k++)
                        memset(work2[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<col;m++)
                            work2[k][l]+=result_par2_two[par2_two[mis_trio_cat_1]][k][m]*result_par1_two[par1_two[mis_trio_cat_2]][l][m]*genofreq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=work2[k][l];      
                    }
                    else if(ind_1par2==ind_2par2){
                      for(int k=0;k<h_size-1;k++)
                        memset(work2[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<col;m++)
                            work2[k][l]+=result_par2_two[par2_two[mis_trio_cat_1]][k][m]*result_par2_two[par2_two[mis_trio_cat_2]][l][m]*genofreq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=work2[k][l];   
                    }     
                  }
                  else if(index==1){
                    iter2 = (fam_list[fam]->IBD2_coeff).find(pair<int, int>(ind_1par1,ind_2par1));
                    iter1 = (fam_list[fam]->IBD1_coeff).find(pair<int, int>(ind_1par1,ind_2par1));
                    IBD_2 = iter2->second;
                    IBD_1 = iter1->second;
                    
                    if(IBD_1>1e-8){
                      for(int k=0;k<h_size-1;k++)
                        memset(work1[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<h_size;m++)
                            work1[k][l]+=result_par1_one[par1_one[mis_trio_cat_1]][k][m]*result_par1_one[par1_one[mis_trio_cat_2]][l][m]*freq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=IBD_1*work1[k][l];        
                    }
                    if(IBD_2>1e-8){
                      for(int k=0;k<h_size-1;k++)
                        memset(work2[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<col;m++)
                            work2[k][l]+=result_par1_two[par1_two[mis_trio_cat_1]][k][m]*result_par1_two[par1_two[mis_trio_cat_2]][l][m]*genofreq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=IBD_2*work2[k][l];      
                    }   
                  }
                  else if(index==2){
                    iter2 = (fam_list[fam]->IBD2_coeff).find(pair<int, int>(ind_1par1,ind_2par2));
                    iter1 = (fam_list[fam]->IBD1_coeff).find(pair<int, int>(ind_1par1,ind_2par2));
                    IBD_2 = iter2->second;
                    IBD_1 = iter1->second;
                    
                    if(IBD_1>1e-8){
                      for(int k=0;k<h_size-1;k++)
                        memset(work1[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<h_size;m++)
                            work1[k][l]+=result_par1_one[par1_one[mis_trio_cat_1]][k][m]*result_par2_one[par2_one[mis_trio_cat_2]][l][m]*freq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=IBD_1*work1[k][l];        
                    }
                    if(IBD_2>1e-8){
                      for(int k=0;k<h_size-1;k++)
                        memset(work2[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<col;m++)
                            work2[k][l]+=result_par1_two[par1_two[mis_trio_cat_1]][k][m]*result_par2_two[par2_two[mis_trio_cat_2]][l][m]*genofreq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=IBD_2*work2[k][l];      
                    }   
                  }
                  else if(index==3){
                    iter2 = (fam_list[fam]->IBD2_coeff).find(pair<int, int>(ind_1par2,ind_2par1));
                    iter1 = (fam_list[fam]->IBD1_coeff).find(pair<int, int>(ind_1par2,ind_2par1));
                    IBD_2 = iter2->second;
                    IBD_1 = iter1->second;
                    
                    if(IBD_1>1e-8){
                      for(int k=0;k<h_size-1;k++)
                        memset(work1[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<h_size;m++)
                            work1[k][l]+=result_par2_one[par2_one[mis_trio_cat_1]][k][m]*result_par1_one[par1_one[mis_trio_cat_2]][l][m]*freq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=IBD_1*work1[k][l];        
                    }
                    if(IBD_2>1e-8){
                      for(int k=0;k<h_size-1;k++)
                        memset(work2[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<col;m++)
                            work2[k][l]+=result_par2_two[par2_two[mis_trio_cat_1]][k][m]*result_par1_two[par1_two[mis_trio_cat_2]][l][m]*genofreq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=IBD_2*work2[k][l];      
                    }   
                  }
                  else if(index==4){
                    iter2 = (fam_list[fam]->IBD2_coeff).find(pair<int, int>(ind_1par2,ind_2par2));
                    iter1 = (fam_list[fam]->IBD1_coeff).find(pair<int, int>(ind_1par2,ind_2par2));
                    IBD_2 = iter2->second;
                    IBD_1 = iter1->second;
                    
                    if(IBD_1>1e-8){
                      for(int k=0;k<h_size-1;k++)
                        memset(work1[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<h_size;m++)
                            work1[k][l]+=result_par2_one[par2_one[mis_trio_cat_1]][k][m]*result_par2_one[par2_one[mis_trio_cat_2]][l][m]*freq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=IBD_1*work1[k][l];        
                    }
                    if(IBD_2>1e-8){
                      for(int k=0;k<h_size-1;k++)
                        memset(work2[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<col;m++)
                            work2[k][l]+=result_par2_two[par2_two[mis_trio_cat_1]][k][m]*result_par2_two[par2_two[mis_trio_cat_2]][l][m]*genofreq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=IBD_2*work2[k][l];      
                    }   
                  }
                  else if(index==5){
                    iter2 = (fam_list[fam]->IBD2_coeff).find(pair<int, int>(ind_1par1,ind_2));
                    iter1 = (fam_list[fam]->IBD1_coeff).find(pair<int, int>(ind_1par1,ind_2));
                    IBD_2 = iter2->second;
                    IBD_1 = iter1->second;
                       
                    if(IBD_1>1e-8){
                      for(int k=0;k<h_size-1;k++)
                        memset(work1[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<h_size;m++)
                            work1[k][l]+=result_par1_one[par1_one[mis_trio_cat_1]][k][m]*result_off_one[off_one[mis_trio_cat_2]][l][m]*freq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=IBD_1*work1[k][l];        
                    }
                    if(IBD_2>1e-8){
                      for(int k=0;k<h_size-1;k++)
                        memset(work2[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<col;m++)
                            work2[k][l]+=result_par1_two[par1_two[mis_trio_cat_1]][k][m]*result_off_two[off_two[mis_trio_cat_2]][l][m]*genofreq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=IBD_2*work2[k][l];      
                    }     
                  }
                  else if(index==6){
                    iter2 = (fam_list[fam]->IBD2_coeff).find(pair<int, int>(ind_1par2,ind_2));
                    iter1 = (fam_list[fam]->IBD1_coeff).find(pair<int, int>(ind_1par2,ind_2));
                    IBD_2 = iter2->second;
                    IBD_1 = iter1->second;
                       
                    if(IBD_1>1e-8){
                      for(int k=0;k<h_size-1;k++)
                        memset(work1[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<h_size;m++)
                            work1[k][l]+=result_par2_one[par2_one[mis_trio_cat_1]][k][m]*result_off_one[off_one[mis_trio_cat_2]][l][m]*freq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=IBD_1*work1[k][l];        
                    }
                    if(IBD_2>1e-8){
                      for(int k=0;k<h_size-1;k++)
                        memset(work2[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<col;m++)
                            work2[k][l]+=result_par2_two[par2_two[mis_trio_cat_1]][k][m]*result_off_two[off_two[mis_trio_cat_2]][l][m]*genofreq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=IBD_2*work2[k][l];      
                    }     
                  }
                  else if(index==7){
                    iter2 = (fam_list[fam]->IBD2_coeff).find(pair<int, int>(ind_1,ind_2par1));
                    iter1 = (fam_list[fam]->IBD1_coeff).find(pair<int, int>(ind_1,ind_2par1));
                    IBD_2 = iter2->second;
                    IBD_1 = iter1->second;
                       
                    if(IBD_1>1e-8){
                      for(int k=0;k<h_size-1;k++)
                        memset(work1[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<h_size;m++)
                            work1[k][l]+=result_off_one[off_one[mis_trio_cat_1]][k][m]*result_par1_one[par1_one[mis_trio_cat_2]][l][m]*freq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=IBD_1*work1[k][l];        
                    }
                    if(IBD_2>1e-8){
                      for(int k=0;k<h_size-1;k++)
                        memset(work2[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<col;m++)
                            work2[k][l]+=result_off_two[off_two[mis_trio_cat_1]][k][m]*result_par1_two[par1_two[mis_trio_cat_2]][l][m]*genofreq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=IBD_2*work2[k][l];      
                    }     
                  }
                  else if(index==8){
                    iter2 = (fam_list[fam]->IBD2_coeff).find(pair<int, int>(ind_1,ind_2par2));
                    iter1 = (fam_list[fam]->IBD1_coeff).find(pair<int, int>(ind_1,ind_2par2));
                    IBD_2 = iter2->second;
                    IBD_1 = iter1->second;
                       
                    if(IBD_1>1e-8){
                      for(int k=0;k<h_size-1;k++)
                        memset(work1[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<h_size;m++)
                            work1[k][l]+=result_off_one[off_one[mis_trio_cat_1]][k][m]*result_par2_one[par2_one[mis_trio_cat_2]][l][m]*freq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=IBD_1*work1[k][l];        
                    }
                    if(IBD_2>1e-8){
                      for(int k=0;k<h_size-1;k++)
                        memset(work2[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<col;m++)
                            work2[k][l]+=result_off_two[off_two[mis_trio_cat_1]][k][m]*result_par2_two[par2_two[mis_trio_cat_2]][l][m]*genofreq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=IBD_2*work2[k][l];      
                    }     
                  }     
                  else if(index==9){
                    if(diff==1){  // j and if are sibs
                      for(int k=0;k<h_size-1;k++)
                        memset(work2[k],0,(h_size-1)*sizeof(double));
                      
                      loc=0;
                      giter = geno_freq_table.begin();
                      while(giter!=geno_freq_table.end()){
                        int h1 = (giter->first).first;
                        int h2 = (giter->first).second;
                        if(h1==h2){
                          for(int k=0;k<dim_1;k++)
                            for(int l=0;l<dim_2;l++)
                              work2[k][l]+=result_par1_two[par1_two[mis_trio_cat_1]][k][loc]*result_parents_one[parents_one[mis_trio_cat_2]][l][hap_map[h1]][hap_map[h2]]*genofreq[loc];
                        }
                        else{
                          for(int k=0;k<dim_1;k++)
                            for(int l=0;l<dim_2;l++)
                              work2[k][l]+=result_par1_two[par1_two[mis_trio_cat_1]][k][loc]*(result_parents_one[parents_one[mis_trio_cat_2]][l][hap_map[h1]][hap_map[h2]]+result_parents_one[parents_one[mis_trio_cat_2]][l][hap_map[h2]][hap_map[h1]])*genofreq[loc]/2.0;
                        }
                        loc++;
                        giter++;
                      }

                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=work2[k][l];
                    }
                    else{
                      iter1 = (fam_list[fam]->IBD1_coeff).find(pair<int, int>(ind_1par1,ind_2par1));
                      IBD_1 = iter1->second;
                    
                      if(IBD_1>1e-8){
                        for(int k=0;k<h_size-1;k++)
                          memset(work1[k],0,(h_size-1)*sizeof(double));
                      
                        for(int k=0;k<dim_1;k++)
                          for(int l=0;l<dim_2;l++)
                            for(int m=0;m<h_size;m++)
                              work1[k][l]+=result_par1_one[par1_one[mis_trio_cat_1]][k][m]*result_par1_one[par1_one[mis_trio_cat_2]][l][m]*freq[m];
                          
                        for(int k=0;k<dim_1;k++)
                          for(int l=0;l<dim_2;l++)
                            covMatrix[nb1+k][nb2+l]+=IBD_1*work1[k][l];
                      }
                    
                      iter1 = (fam_list[fam]->IBD1_coeff).find(pair<int, int>(ind_1par1,ind_2par2));
                      IBD_1 = iter1->second;
                    
                      if(IBD_1>1e-8){
                        for(int k=0;k<h_size-1;k++)
                          memset(work1[k],0,(h_size-1)*sizeof(double));
                      
                        for(int k=0;k<dim_1;k++)
                          for(int l=0;l<dim_2;l++)
                            for(int m=0;m<h_size;m++)
                              work1[k][l]+=result_par1_one[par1_one[mis_trio_cat_1]][k][m]*result_par2_one[par2_one[mis_trio_cat_2]][l][m]*freq[m];
                          
                        for(int k=0;k<dim_1;k++)
                          for(int l=0;l<dim_2;l++)
                            covMatrix[nb1+k][nb2+l]+=IBD_1*work1[k][l];
                      }
                    }                          
                  }
                  else if(index==10){
                    if(diff==1){  // j and im are sibs
                      for(int k=0;k<h_size-1;k++)
                        memset(work2[k],0,(h_size-1)*sizeof(double));
                      
                      loc=0;
                      giter = geno_freq_table.begin();
                      while(giter!=geno_freq_table.end()){
                        int h1 = (giter->first).first;
                        int h2 = (giter->first).second;
                        if(h1==h2){
                          for(int k=0;k<dim_1;k++)
                            for(int l=0;l<dim_2;l++)
                              work2[k][l]+=result_par2_two[par2_two[mis_trio_cat_1]][k][loc]*result_parents_one[parents_one[mis_trio_cat_2]][l][hap_map[h1]][hap_map[h2]]*genofreq[loc];
                        }
                        else{
                          for(int k=0;k<dim_1;k++)
                            for(int l=0;l<dim_2;l++)
                              work2[k][l]+=result_par2_two[par2_two[mis_trio_cat_1]][k][loc]*(result_parents_one[parents_one[mis_trio_cat_2]][l][hap_map[h1]][hap_map[h2]]+result_parents_one[parents_one[mis_trio_cat_2]][l][hap_map[h2]][hap_map[h1]])*genofreq[loc]/2.0;
                        }
                        loc++;
                        giter++;
                      }

                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=work2[k][l];
                    }
                    else{
                      iter1 = (fam_list[fam]->IBD1_coeff).find(pair<int, int>(ind_1par2,ind_2par1));
                      IBD_1 = iter1->second;
                    
                      if(IBD_1>1e-8){
                        for(int k=0;k<h_size-1;k++)
                          memset(work1[k],0,(h_size-1)*sizeof(double));
                      
                        for(int k=0;k<dim_1;k++)
                          for(int l=0;l<dim_2;l++)
                            for(int m=0;m<h_size;m++)
                              work1[k][l]+=result_par2_one[par2_one[mis_trio_cat_1]][k][m]*result_par1_one[par1_one[mis_trio_cat_2]][l][m]*freq[m];
                          
                        for(int k=0;k<dim_1;k++)
                          for(int l=0;l<dim_2;l++)
                            covMatrix[nb1+k][nb2+l]+=IBD_1*work1[k][l];
                      }
                    
                      iter1 = (fam_list[fam]->IBD1_coeff).find(pair<int, int>(ind_1par2,ind_2par2));
                      IBD_1 = iter1->second;
                    
                      if(IBD_1>1e-8){
                        for(int k=0;k<h_size-1;k++)
                          memset(work1[k],0,(h_size-1)*sizeof(double));
                      
                        for(int k=0;k<dim_1;k++)
                          for(int l=0;l<dim_2;l++)
                            for(int m=0;m<h_size;m++)
                              work1[k][l]+=result_par2_one[par2_one[mis_trio_cat_1]][k][m]*result_par2_one[par2_one[mis_trio_cat_2]][l][m]*freq[m];
                          
                        for(int k=0;k<dim_1;k++)
                          for(int l=0;l<dim_2;l++)
                            covMatrix[nb1+k][nb2+l]+=IBD_1*work1[k][l];
                      }
                    }                          
                  }
                  else if(index==11){
                    if(diff==1){  // i and jf are sibs
                      for(int k=0;k<h_size-1;k++)
                        memset(work2[k],0,(h_size-1)*sizeof(double));
                      
                      loc=0;
                      giter = geno_freq_table.begin();
                      while(giter!=geno_freq_table.end()){
                        int h1 = (giter->first).first;
                        int h2 = (giter->first).second;
                        if(h1==h2){
                          for(int k=0;k<dim_1;k++)
                            for(int l=0;l<dim_2;l++)
                              work2[k][l]+=result_parents_one[parents_one[mis_trio_cat_1]][k][hap_map[h1]][hap_map[h2]]*result_par1_two[par1_two[mis_trio_cat_2]][l][loc]*genofreq[loc];
                        }
                        else{
                          for(int k=0;k<dim_1;k++)
                            for(int l=0;l<dim_2;l++)
                              work2[k][l]+=(result_parents_one[parents_one[mis_trio_cat_1]][k][hap_map[h1]][hap_map[h2]]+result_parents_one[parents_one[mis_trio_cat_1]][k][hap_map[h2]][hap_map[h1]])*result_par1_two[par1_two[mis_trio_cat_2]][l][loc]*genofreq[loc]/2.0;
                        }
                        loc++;
                        giter++;
                      }

                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=work2[k][l];
                    }
                    else{
                      iter1 = (fam_list[fam]->IBD1_coeff).find(pair<int, int>(ind_1par1,ind_2par1));
                      IBD_1 = iter1->second;
                    
                      if(IBD_1>1e-8){
                        for(int k=0;k<h_size-1;k++)
                          memset(work1[k],0,(h_size-1)*sizeof(double));
                      
                        for(int k=0;k<dim_1;k++)
                          for(int l=0;l<dim_2;l++)
                            for(int m=0;m<h_size;m++)
                              work1[k][l]+=result_par1_one[par1_one[mis_trio_cat_1]][k][m]*result_par1_one[par1_one[mis_trio_cat_2]][l][m]*freq[m];
                          
                        for(int k=0;k<dim_1;k++)
                          for(int l=0;l<dim_2;l++)
                            covMatrix[nb1+k][nb2+l]+=IBD_1*work1[k][l];
                      }
                    
                      iter1 = (fam_list[fam]->IBD1_coeff).find(pair<int, int>(ind_1par2,ind_2par1));
                      IBD_1 = iter1->second;
                    
                      if(IBD_1>1e-8){
                        for(int k=0;k<h_size-1;k++)
                          memset(work1[k],0,(h_size-1)*sizeof(double));
                      
                        for(int k=0;k<dim_1;k++)
                          for(int l=0;l<dim_2;l++)
                            for(int m=0;m<h_size;m++)
                              work1[k][l]+=result_par2_one[par2_one[mis_trio_cat_1]][k][m]*result_par1_one[par1_one[mis_trio_cat_2]][l][m]*freq[m];
                          
                        for(int k=0;k<dim_1;k++)
                          for(int l=0;l<dim_2;l++)
                            covMatrix[nb1+k][nb2+l]+=IBD_1*work1[k][l];
                      }
                    }                          
                  }
                  else if(index==12){
                    if(diff==1){  // i and jm are sibs
                      for(int k=0;k<h_size-1;k++)
                        memset(work2[k],0,(h_size-1)*sizeof(double));
                      
                      loc=0;
                      giter = geno_freq_table.begin();
                      while(giter!=geno_freq_table.end()){
                        int h1 = (giter->first).first;
                        int h2 = (giter->first).second;
                        if(h1==h2){
                          for(int k=0;k<dim_1;k++)
                            for(int l=0;l<dim_2;l++)
                              work2[k][l]+=result_parents_one[parents_one[mis_trio_cat_1]][k][hap_map[h1]][hap_map[h2]]*result_par2_two[par2_two[mis_trio_cat_2]][l][loc]*genofreq[loc];
                        }
                        else{
                          for(int k=0;k<dim_1;k++)
                            for(int l=0;l<dim_2;l++)
                              work2[k][l]+=(result_parents_one[parents_one[mis_trio_cat_1]][k][hap_map[h1]][hap_map[h2]]+result_parents_one[parents_one[mis_trio_cat_1]][k][hap_map[h2]][hap_map[h1]])*result_par2_two[par2_two[mis_trio_cat_2]][l][loc]*genofreq[loc]/2.0;
                        }
                        loc++;
                        giter++;
                      }

                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=work2[k][l];
                    }
                    else{
                      iter1 = (fam_list[fam]->IBD1_coeff).find(pair<int, int>(ind_1par1,ind_2par2));
                      IBD_1 = iter1->second;
                    
                      if(IBD_1>1e-8){
                        for(int k=0;k<h_size-1;k++)
                          memset(work1[k],0,(h_size-1)*sizeof(double));
                      
                        for(int k=0;k<dim_1;k++)
                          for(int l=0;l<dim_2;l++)
                            for(int m=0;m<h_size;m++)
                              work1[k][l]+=result_par1_one[par1_one[mis_trio_cat_1]][k][m]*result_par2_one[par2_one[mis_trio_cat_2]][l][m]*freq[m];
                          
                        for(int k=0;k<dim_1;k++)
                          for(int l=0;l<dim_2;l++)
                            covMatrix[nb1+k][nb2+l]+=IBD_1*work1[k][l];
                      }
                    
                      iter1 = (fam_list[fam]->IBD1_coeff).find(pair<int, int>(ind_1par2,ind_2par2));
                      IBD_1 = iter1->second;
                    
                      if(IBD_1>1e-8){
                        for(int k=0;k<h_size-1;k++)
                          memset(work1[k],0,(h_size-1)*sizeof(double));
                      
                        for(int k=0;k<dim_1;k++)
                          for(int l=0;l<dim_2;l++)
                            for(int m=0;m<h_size;m++)
                              work1[k][l]+=result_par2_one[par2_one[mis_trio_cat_1]][k][m]*result_par2_one[par2_one[mis_trio_cat_2]][l][m]*freq[m];
                          
                        for(int k=0;k<dim_1;k++)
                          for(int l=0;l<dim_2;l++)
                            covMatrix[nb1+k][nb2+l]+=IBD_1*work1[k][l];
                      }
                    }                          
                  }   
                  // index 19 & 20
                }     
                for(int k=0;k<dim_1;k++)
                  for(int l=0;l<dim_2;l++)
                    covMatrix[nb2+l][nb1+k] = covMatrix[nb1+k][nb2+l];
              }  // i, j related
              nb2+=dim_2;
            }  // j is typed
          }  // complete j

          const vector<double> hap_num = geno_trio_list[trio_cat_1]->get_hap_num();
          test_weight = fam_list[fam]->test_weight[i];

          if(missing_trio_list[mis_trio_cat_1]->typed_par1==1 || missing_trio_list[mis_trio_cat_1]->typed_par2==1){
            test_trio_weight1 = fam_list[fam]->test_trio_weight1[i];
            test_trio_weight2 = fam_list[fam]->test_trio_weight2[i];
            test_trio_weight3 = fam_list[fam]->test_trio_weight3[i];
            for(int j=0;j<dim_1;j++){
              for(int k=0;k<h_size-1;k++)
                cholent[nb1+j][k] = missing_trio_list[mis_trio_cat_1]->derivative[j][k];
              cholent[nb1+j][h_size-1] = hap_num[j];
              cholent[nb1+j][h_size] = test_trio_weight1*(missing_trio_list[mis_trio_cat_1]->deriv1[j])+test_trio_weight2*(missing_trio_list[mis_trio_cat_1]->deriv2[j])+test_trio_weight3*(missing_trio_list[mis_trio_cat_1]->deriv3[j]);
              cholent[nb1+j][h_size]+=test_weight*freq[missing_trio_list[mis_trio_cat_1]->sites[j]]*h_weight[missing_trio_list[mis_trio_cat_1]->sites[j]];
            }
          }
          else{
            for(int j=0;j<dim_1;j++){
              for(int k=0;k<h_size-1;k++)
                cholent[nb1+j][k] = missing_trio_list[mis_trio_cat_1]->derivative[j][k];
              cholent[nb1+j][h_size-1] = hap_num[j];
              cholent[nb1+j][h_size] = test_weight*missing_trio_list[mis_trio_cat_1]->deriv_r[j];
            }
          }
          nb1+=dim_1;
        }  // i is typed
      }  // complete i

      if(cholesky(covMatrix,dim,cholent,h_size+1,chol,cholaug,0) != 1){
        //printf("cholesky decomposition of the cov matrix failed for family %d. \n", fam_list[fam]->fam_id);
        flag_pd = 0;
        //break;
      }

      if(flag_pd==1){
      double **s = self_colmultiply(cholaug,dim,1,h_size-1,1,h_size-1);
      double **t = self_colmultiply(cholaug,dim,h_size+1,h_size+1,h_size+1,h_size+1);
      double **u = self_colmultiply(cholaug,dim,1,h_size-1,h_size+1,h_size+1);
      double **v = self_colmultiply(cholaug,dim,h_size+1,h_size+1,h_size,h_size);

      for(int i=0;i<h_size-1;i++)
        for(int j=0;j<h_size-1;j++)
          denominator[i][j]+=s[i][j];

      A1+=t[0][0];

      for(int i=0;i<h_size-1;i++)
        A2[i][0]+=u[i][0];

      numerator_1+=v[0][0];

      clear_matrix(s,h_size-1,h_size-1);
      clear_matrix(t,1,1);
      clear_matrix(u,h_size-1,1);
      clear_matrix(v,1,1);
      }
      else{
      double var;
      double &determ = var;
      double **v = inverse(covMatrix,dim,determ);
      double **s = quadratic_product(cholent,dim,h_size+1,v,dim,dim,cholent,dim,h_size+1);

      for(int i=0;i<h_size-1;i++)
        for(int j=0;j<h_size-1;j++)
          denominator[i][j]+=s[i][j];

      A1+=s[h_size][h_size];

      for(int i=0;i<h_size-1;i++)
        A2[i][0]+=s[i][h_size];

      numerator_1+=s[h_size][h_size-1];

      clear_matrix(s,h_size+1,h_size+1);
      clear_matrix(v,dim,dim);
      }
    }
  }

  h_weight.clear();

  double **cholmain = new double*[h_size-1];
  for(int i=0;i<h_size-1;i++){
    cholmain[i] = new double[h_size-1];
    memset(cholmain[i],0,(h_size-1)*sizeof(double));
  }

  double **cholaugmain = new double*[h_size-1];
  for(int i=0;i<h_size-1;i++){
    cholaugmain[i] = new double[1];
    memset(cholaugmain[i],0,sizeof(double));
  }

  if(cholesky(denominator,h_size-1,A2,1,cholmain,cholaugmain,0) != 1){
    //printf("cholesky decomposition of the denominator matrix failed.\n");
    flag_pd==0;
    //exit(1);
  }

  if(flag_pd==1){
  double **s = self_colmultiply(cholaugmain,h_size-1,1,1,1,1);

  double orig = A1-s[0][0];
  clear_matrix(s,1,1);
  
  test = numerator_1*numerator_1/orig;
  }
  else{
  double var;
  double &determ = var;
  double **v = inverse(denominator,h_size-1,determ);
  double **s = quadratic_product(A2,h_size-1,1,v,h_size-1,h_size-1,A2,h_size-1,1);

  double orig = A1-s[0][0];
  test = numerator_1*numerator_1/orig;

  clear_matrix(s,1,1);
  clear_matrix(v,h_size-1,h_size-1);
  }
  
  clear_matrix(cholmain,h_size-1,h_size-1);
  clear_matrix(cholaugmain,h_size-1,1);
  
  pvalue = pochisq(test,1);

  clear_matrix(A2,h_size-1,1);
}


int IQLBEstimator::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)
{
  // assume ind_1 != ind_2
  if(ind_1==ind_2f)
    return 4;
  else if(ind_1==ind_2m)
    return 5;
  else if(fabs(IBD_2-0.0625)<1e-8 && fabs(IBD_1-0.375)<1e-8) // double first cousin
    return 7;
  else if(fabs(IBD_2-0.125)<1e-8 && fabs(IBD_1-0.5)<1e-8) // half-sib plus first cousin
    return 8;
  else if(fabs(IBD_2-0.25)<1e-8 && fabs(IBD_1-0.5)<1e-8) // sib
  {
    diff = 1;
    return 6;
  }
  else  // four persons
  {
    double value_f, value_m;
    map<pair<int, int>, double>::iterator iter;
    iter = (fam_list[fam]->kinship_coeff).find(pair<int, int>(ind_1,ind_2f));
    if(iter == (fam_list[fam]->kinship_coeff).end())
    {
      printf("No kinship coefficient between individual %d and individual %d from family %d. Please check...\n\n", ind_1,ind_2f,fam_list[fam]->fam_id);
      exit(1);
    }
    value_f = iter->second;
    iter = (fam_list[fam]->kinship_coeff).find(pair<int, int>(ind_1,ind_2m));
    if(iter == (fam_list[fam]->kinship_coeff).end())
    {
      printf("No kinship coefficient between individual %d and individual %d from family %d. Please check...\n\n", ind_1,ind_2m,fam_list[fam]->fam_id);
      exit(1);
    }
    value_m = iter->second;

    if(value_m<1e-8) // i is unrelated to jm
      return 1;
    else if(value_f<1e-8) // i is unrelated to jf
      return 2;
    else if(value_o>value_f && value_o>value_m)  // j is an ancestor of i
      return 3;
    else if(IBD_2<1e-8) // jf and jm are great...great grandparents of i
    {
      double gen = log(2.0/IBD_1)/log(2.0);
      diff = (gen-floor(gen)>=0.5) ? int(ceil(gen)) : int(floor(gen));
      return 6;
    }
    return 0;
  }
}


int IQLBEstimator::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)
{
  // assume ind_1 != ind_2
  if(ind_1==ind_2f)
    return 13;
  else if(ind_1==ind_2m)
    return 14;
  else if(ind_2==ind_1f)
    return 15;
  else if(ind_2==ind_1m)
    return 16;
  else if(ind_1f==ind_2f && ind_1m==ind_2m) // sib
    return 17;
  else if((ind_1f==ind_2f || ind_1m==ind_2m || ind_1f==ind_2m || ind_1m==ind_2f) && IBD_2<1e-8 && fabs(IBD_1-0.5)<1e-8) // half-sib
    return 18;
  else if(fabs(IBD_2-0.0625)<1e-8 && fabs(IBD_1-0.375)<1e-8) // double first cousin
    return 19;
  else if(fabs(IBD_2-0.125)<1e-8 && fabs(IBD_1-0.5)<1e-8) // half-sib plus first cousin
    return 20;
  else
  {
    double value_of, value_om, value_fo, value_mo, value_ff, value_fm, value_mf, value_mm;
    map<pair<int, int>, double>::iterator iter;
    iter = (fam_list[fam]->kinship_coeff).find(pair<int, int>(ind_1,ind_2f));
    if(iter == (fam_list[fam]->kinship_coeff).end())
    {
      printf("No kinship coefficient between individual %d and individual %d from family %d. Please check...\n\n", ind_1,ind_2f,fam_list[fam]->fam_id);
      exit(1);
    }
    value_of = iter->second;
    iter = (fam_list[fam]->kinship_coeff).find(pair<int, int>(ind_1,ind_2m));
    if(iter == (fam_list[fam]->kinship_coeff).end())
    {
      printf("No kinship coefficient between individual %d and individual %d from family %d. Please check...\n\n", ind_1,ind_2m,fam_list[fam]->fam_id);
      exit(1);
    }
    value_om = iter->second;
    iter = (fam_list[fam]->kinship_coeff).find(pair<int, int>(ind_2,ind_1f));
    if(iter == (fam_list[fam]->kinship_coeff).end())
    {
      printf("No kinship coefficient between individual %d and individual %d from family %d. Please check...\n\n", ind_2,ind_1f,fam_list[fam]->fam_id);
      exit(1);
    }
    value_fo = iter->second;
    iter = (fam_list[fam]->kinship_coeff).find(pair<int, int>(ind_2,ind_1m));
    if(iter == (fam_list[fam]->kinship_coeff).end())
    {
      printf("No kinship coefficient between individual %d and individual %d from family %d. Please check...\n\n", ind_2,ind_1m,fam_list[fam]->fam_id);
      exit(1);
    }
    value_mo = iter->second;
    iter = (fam_list[fam]->kinship_coeff).find(pair<int, int>(ind_1f,ind_2f));
    if(iter == (fam_list[fam]->kinship_coeff).end())
    {
      printf("No kinship coefficient between individual %d and individual %d from family %d. Please check...\n\n", ind_1f,ind_2f,fam_list[fam]->fam_id);
      exit(1);
    }
    value_ff = iter->second;
    iter = (fam_list[fam]->kinship_coeff).find(pair<int, int>(ind_1f,ind_2m));
    if(iter == (fam_list[fam]->kinship_coeff).end())
    {
      printf("No kinship coefficient between individual %d and individual %d from family %d. Please check...\n\n", ind_1f,ind_2m,fam_list[fam]->fam_id);
      exit(1);
    }
    value_fm = iter->second;
    iter = (fam_list[fam]->kinship_coeff).find(pair<int, int>(ind_1m,ind_2f));
    if(iter == (fam_list[fam]->kinship_coeff).end())
    {
      printf("No kinship coefficient between individual %d and individual %d from family %d. Please check...\n\n", ind_1m,ind_2f,fam_list[fam]->fam_id);
      exit(1);
    }
    value_mf = iter->second;
    iter = (fam_list[fam]->kinship_coeff).find(pair<int, int>(ind_1m,ind_2m));
    if(iter == (fam_list[fam]->kinship_coeff).end())
    {
      printf("No kinship coefficient between individual %d and individual %d from family %d. Please check...\n\n", ind_1m,ind_2m,fam_list[fam]->fam_id);
      exit(1);
    }
    value_mm = iter->second;

    if(value_om<1e-8 && value_mo<1e-8) // if is related to jf
      return 1;
    else if(value_of<1e-8 && value_mo<1e-8) // if is related to jm
      return 2;
    else if(value_om<1e-8 && value_fo<1e-8) // im is related to jf
      return 3;
    else if(value_of<1e-8 && value_fo<1e-8) // im is related to jm
      return 4;

    else if(value_mo<1e-8 && value_fo>value_ff && value_fo>value_fm)  // j is an ancestor of if
      return 5;
    else if(value_fo<1e-8 && value_mo>value_mf && value_mo>value_mm)  // j is an ancestor of im
      return 6;
    else if(value_om<1e-8 && value_of>value_ff && value_of>value_mf)  // i is an ancestor of jf
      return 7;
    else if(value_of<1e-8 && value_om>value_fm && value_om>value_mm)  // i is an ancestor of jm
      return 8;

    else if(value_mo<1e-8 && IBD_2<1e-8) // jf and jm are great...great grandparents of i, if
    {
      double gen = log(2.0/IBD_1)/log(2.0)-1;
      diff = (gen-floor(gen)>=0.5) ? int(ceil(gen)) : int(floor(gen));
      return 9;
    }
    else if(value_fo<1e-8 && IBD_2<1e-8) // jf and jm are great...great grandparents of i, im
    {
      double gen = log(2.0/IBD_1)/log(2.0)-1;
      diff = (gen-floor(gen)>=0.5) ? int(ceil(gen)) : int(floor(gen));
      return 10;
    }
    else if(value_om<1e-8 && IBD_2<1e-8) // if and im are great...great grandparents of j, jf
    {
      double gen = log(2.0/IBD_1)/log(2.0)-1;
      diff = (gen-floor(gen)>=0.5) ? int(ceil(gen)) : int(floor(gen));
      return 11;
    }
    else if(value_of<1e-8 && IBD_2<1e-8) // if and im are great...great grandparents of j, jm
    {
      double gen = log(2.0/IBD_1)/log(2.0)-1;
      diff = (gen-floor(gen)>=0.5) ? int(ceil(gen)) : int(floor(gen));
      return 12;
    }
    return 0;
  }
}


void IQLBEstimator::estimate_subgroup(){
     
  full_hap_freq_table = hap_freq_table;
  
  // affected  
  // unaffected  
  // unknown
  aff=0;
  unaff=0;
  unknown=0;

  int trio_cat_1;
  int mis_trio_cat_1;
  int dim_1;
  for(int fam=0;fam<fam_list.size();fam++){
    if(fam_list[fam]->N_typed>0){
      for(int i=0;i<fam_list[fam]->N;i++){
        trio_cat_1 = fam_list[fam]->trio_cat[i];
        mis_trio_cat_1 = fam_list[fam]->mis_trio_cat[i];
        dim_1 = missing_trio_list[mis_trio_cat_1]->dim;
        if(dim_1 > 0 && geno_trio_list[trio_cat_1]->info==1){
          if(fam_list[fam]->fam_member[i][4]==0)
            unknown++;
          else if(fam_list[fam]->fam_member[i][4]==1)
            unaff++;
          else if(fam_list[fam]->fam_member[i][4]==2)
            aff++;
        }              
      }                                   
    }
  }
  
  if(aff>=15){
    reset_freq_table(full_hap_freq_table);
    int count=0;
    int flag;
    while(1){
      map<int, double> of = hap_freq_table;
      NR_Iteration_subgroup(2);
      map<int, double> nf = hap_freq_table;

      double diff = 0;
      flag = 1;
      for(int i=0;i<of.size();i++){
        diff += (nf[i]-of[i])*(nf[i]-of[i]);
        if(nf[i]<0){
          flag = 0;
          break;
        }
      }

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

      count++;
    }
    aff_hap_freq_table = hap_freq_table;
  }
  if(unaff>=15){
    reset_freq_table(full_hap_freq_table);
    int count=0;
    int flag;
    while(1){
      map<int, double> of = hap_freq_table;
      NR_Iteration_subgroup(1);
      map<int, double> nf = hap_freq_table;

      double diff = 0;
      flag = 1;
      for(int i=0;i<of.size();i++){
        diff += (nf[i]-of[i])*(nf[i]-of[i]);
        if(nf[i]<0){
          flag = 0;
          break;
        }
      }

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

      count++;
    }
    unaff_hap_freq_table = hap_freq_table;
  }
  if(unknown>=15){
    reset_freq_table(full_hap_freq_table);
    int count=0;
    int flag;
    while(1){
      map<int, double> of = hap_freq_table;
      NR_Iteration_subgroup(0);
      map<int, double> nf = hap_freq_table;

      double diff = 0;
      flag = 1;
      for(int i=0;i<of.size();i++){
        diff += (nf[i]-of[i])*(nf[i]-of[i]);
        if(nf[i]<0){
          flag = 0;
          break;
        }
      }

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

      count++;
    }
    unknown_hap_freq_table = hap_freq_table;
  }
  
  hap_freq_table = full_hap_freq_table;                                    
}     


void IQLBEstimator::NR_Iteration_subgroup(int group){

  int h_size = hap_freq_table.size();
  int col = geno_freq_table.size();

  for(int i=0;i<h_size-1;i++)
    memset(rst[i],0,col*sizeof(double));  // record derivative for each pair of haplotypes

  int loc = 0;
  map<int, double>::iterator iter=hap_freq_table.begin();
  while(iter!=hap_freq_table.end()){
    freq[loc] = iter->second;
    loc++;
    iter++;
  }
  
  loc = 0;
  map<pair<int, int>, double>::iterator giter=geno_freq_table.begin();
  while(giter!=geno_freq_table.end()){
    genofreq[loc] = giter->second;
    loc++;
    giter++;
  }

  int h_last=0;
  iter = hap_freq_table.end();
  iter--;
  h_last=iter->first;

  giter=geno_freq_table.begin();
  while(giter!=geno_freq_table.end()){
    pair<int, int> conf = giter->first;
    loc = 0;
    iter = hap_freq_table.begin();
    while(loc<h_size-1){
      int h=iter->first;

      if(conf.first == conf.second && conf.first == h)
        rst[loc][hap_pair_map[conf]] = 2*hap_freq_table[h];
      else if(conf.first == conf.second && conf.first == h_last)
        rst[loc][hap_pair_map[conf]] = -2*hap_freq_table[h_last];
      else if(conf.first < conf.second && conf.first == h && conf.second != h_last)
        rst[loc][hap_pair_map[conf]] = 2*hap_freq_table[conf.second];
      else if(conf.first < conf.second && conf.second == h)
        rst[loc][hap_pair_map[conf]] = 2*hap_freq_table[conf.first];
      else if(conf.first == h && conf.second == h_last)
        rst[loc][hap_pair_map[conf]] = 2*(hap_freq_table[h_last]-hap_freq_table[h]);
      else if(conf.first != h && conf.first != h_last && conf.second == h_last)
        rst[loc][hap_pair_map[conf]] = -2*hap_freq_table[conf.first];

      loc++;
      iter++;
    }
    giter++;
  }

  for(int i=0;i<geno_trio_list.size();i++){
    if(geno_trio_list[i]->get_num()>0 && geno_trio_list[i]->info==1 && missing_trio_list[geno_missing_trio_map[i]]->dim>0)
      geno_trio_list[i]->expected_hap_num(hap_freq_table,freq,geno_freq_table,hap_pair_map,count,missing_trio_list[geno_missing_trio_map[i]]->sites);
  }

  for(int i=0;i<missing_trio_list.size();i++){
    if(missing_trio_list[i]->get_num()>0 && missing_trio_list[i]->dim>0){
      missing_trio_list[i]->set_weight(hap_freq_table,geno_freq_table,hap_pair_map,count);
      missing_trio_list[i]->compute_centmatrix(hap_freq_table,freq,geno_freq_table);
      missing_trio_list[i]->compute_derivative(hap_freq_table,geno_freq_table,hap_pair_map,count,rst,hst);
      missing_trio_list[i]->expected_hap_num(hap_freq_table,freq,geno_freq_table,hap_pair_map,count);
    }
  }

  int par1_one_size = par1_one.size();
  if(par1_one_size>0){
    map<int, int>::iterator pi=par1_one.begin();
    for(int i=0;i<par1_one_size;i++){
      pi->second = i;
      int mis = pi->first;
      int dim_1 = missing_trio_list[mis]->dim;
      for(int j=0;j<dim_1;j++)
        memset(result_par1_one[i][j],0,h_size*sizeof(double));
      missing_trio_list[mis]->cond_par1_one(result_par1_one[i],hap_freq_table,geno_freq_table,hap_map);
      pi++;
    }
  }

  int par1_two_size = par1_two.size();
  if(par1_two_size>0){
    map<int, int>::iterator pi=par1_two.begin();
    for(int i=0;i<par1_two_size;i++){
      pi->second = i;
      int mis = pi->first;
      int dim_1 = missing_trio_list[mis]->dim;
      for(int j=0;j<dim_1;j++)
        memset(result_par1_two[i][j],0,col*sizeof(double));
      missing_trio_list[mis]->cond_par1_two(result_par1_two[i],hap_freq_table,geno_freq_table,hap_pair_map);  
      pi++;
    }
  }
  
  int par2_one_size = par2_one.size();
  if(par2_one_size>0){
    map<int, int>::iterator pi=par2_one.begin();
    for(int i=0;i<par2_one_size;i++){
      pi->second = i;
      int mis = pi->first;
      int dim_1 = missing_trio_list[mis]->dim;
      for(int j=0;j<dim_1;j++)
        memset(result_par2_one[i][j],0,h_size*sizeof(double));
      missing_trio_list[mis]->cond_par2_one(result_par2_one[i],hap_freq_table,geno_freq_table,hap_map);  
      pi++;
    }
  }
  
  int par2_two_size = par2_two.size();
  if(par2_two_size>0){
    map<int, int>::iterator pi=par2_two.begin();
    for(int i=0;i<par2_two_size;i++){
      pi->second = i;
      int mis = pi->first;
      int dim_1 = missing_trio_list[mis]->dim;
      for(int j=0;j<dim_1;j++)
        memset(result_par2_two[i][j],0,col*sizeof(double));
      missing_trio_list[mis]->cond_par2_two(result_par2_two[i],hap_freq_table,geno_freq_table,hap_pair_map);  
      pi++;
    }
  }  

  int off_one_size = off_one.size();
  if(off_one_size>0){
    map<int, int>::iterator pi=off_one.begin();
    for(int i=0;i<off_one_size;i++){
      pi->second = i;
      int mis = pi->first;
      int dim_1 = missing_trio_list[mis]->dim;
      for(int j=0;j<dim_1;j++)
        memset(result_off_one[i][j],0,h_size*sizeof(double));
      missing_trio_list[mis]->cond_off_one(result_off_one[i],hap_freq_table,geno_freq_table,hap_map);  
      pi++;
    }
  }
  
  int off_two_size = off_two.size();
  if(off_two_size>0){
    map<int, int>::iterator pi=off_two.begin();
    for(int i=0;i<off_two_size;i++){
      pi->second = i;
      int mis = pi->first;
      int dim_1 = missing_trio_list[mis]->dim;
      for(int j=0;j<dim_1;j++)
        memset(result_off_two[i][j],0,col*sizeof(double));
      missing_trio_list[mis]->cond_off_two(result_off_two[i],hap_freq_table,geno_freq_table,hap_pair_map);
      pi++;
    }
  }  

  int parents_one_size = parents_one.size();
  if(parents_one_size>0){
    map<int, int>::iterator pi=parents_one.begin();
    for(int i=0;i<parents_one_size;i++){
      pi->second = i;
      int mis = pi->first;
      int dim_1 = missing_trio_list[mis]->dim;
      for(int j=0;j<dim_1;j++){
        for(int k=0;k<h_size;k++)
          memset(result_parents_one[i][j][k],0,h_size*sizeof(double));  
      }
      missing_trio_list[mis]->cond_parents_one(result_parents_one[i],hap_freq_table,geno_freq_table,hap_map);
      pi++;
    }
  }
  
  int parents_two_size = parents_two.size();
  if(parents_two_size>0){
    map<int, int>::iterator pi=parents_two.begin();
    for(int i=0;i<parents_two_size;i++){
      pi->second = i;
      int mis = pi->first;
      int dim_1 = missing_trio_list[mis]->dim;
      for(int j=0;j<dim_1;j++){
        for(int k=0;k<col;k++)
          memset(result_parents_two[i][j][k],0,col*sizeof(double));
      }
      missing_trio_list[mis]->cond_parents_two(result_parents_two[i],hap_freq_table,geno_freq_table,hap_pair_map);
      pi++;
    }
  }

  for(int i=0;i<h_size-1;i++)
    memset(denominator[i],0,(h_size-1)*sizeof(double));  // F_p^T Omega^(-1) F_p

  memset(numerator,0,(h_size-1)*sizeof(double));  // F_p^T Omega^(-1) (Z-\mu)

  for(int fam=0;fam<fam_list.size();fam++){
    if(fam_list[fam]->N_typed>0){
      int trio_cat;
      int mis_trio_cat;
      int dim=0;
      for(int i=0;i<fam_list[fam]->N;i++){
        trio_cat = fam_list[fam]->trio_cat[i];
        if(geno_trio_list[trio_cat]->info==1 && fam_list[fam]->fam_member[i][4]==group){
          mis_trio_cat = fam_list[fam]->mis_trio_cat[i];
          dim+=missing_trio_list[mis_trio_cat]->dim;
        }
      }

      if(dim<=0)
        continue;

      for(int i=0;i<tot_dim;i++){
        memset(cholent[i],0,h_size*sizeof(double));
        memset(chol[i],0,tot_dim*sizeof(double));
        memset(cholaug[i],0,h_size*sizeof(double));
        memset(covMatrix[i],0,tot_dim*sizeof(double));
      }

      // populate covariance matrix
      int nb1 = 0;
      int nb2 = 0;
      int ind_1, ind_2, ind_1f, ind_1m, ind_2f, ind_2m;
      int par_1, par_2;  // how many parents is typed
      int dim_1, dim_2;
      int trio_cat_1, trio_cat_2;
      int mis_trio_cat_1, mis_trio_cat_2;
      double IBD_2, IBD_1, kin;
      map<pair<int, int>, double>::iterator iterk;
      map<pair<int, int>, double>::iterator iter1;
      map<pair<int, int>, double>::iterator iter2;
      double value;

      for(int i=0;i<fam_list[fam]->N;i++){
        trio_cat_1 = fam_list[fam]->trio_cat[i];
        mis_trio_cat_1 = fam_list[fam]->mis_trio_cat[i];
        dim_1 = missing_trio_list[mis_trio_cat_1]->dim;
        if(dim_1 > 0 && geno_trio_list[trio_cat_1]->info==1 && fam_list[fam]->fam_member[i][4]==group){
          nb2 = nb1;
          ind_1 = fam_list[fam]->fam_member[i][0];
          ind_2 = ind_1;
          iterk = (fam_list[fam]->kinship_coeff).find(pair<int, int>(ind_1,ind_2));
          if(iterk == (fam_list[fam]->kinship_coeff).end()){
            printf("No IBD coefficient for individual %d from family %d. Please check...\n\n", ind_1,fam_list[fam]->fam_id);
            exit(1);
          }

          for(int j=0;j<dim_1;j++)
            for(int k=0;k<dim_1;k++)
              covMatrix[nb1+j][nb2+k] = missing_trio_list[mis_trio_cat_1]->centmatrix[j][k];

          nb2+=dim_1;

          par_1 = 0;
          par_1+=missing_trio_list[mis_trio_cat_1]->typed_par1;
          par_1+=missing_trio_list[mis_trio_cat_1]->typed_par2;
          ind_1f = fam_list[fam]->fam_member[i][1];
          ind_1m = fam_list[fam]->fam_member[i][2];

          for(int j=i+1;j<fam_list[fam]->N;j++){
            trio_cat_2 = fam_list[fam]->trio_cat[j];
            mis_trio_cat_2 = fam_list[fam]->mis_trio_cat[j];
            dim_2 = missing_trio_list[mis_trio_cat_2]->dim;
            if(dim_2 > 0 && geno_trio_list[trio_cat_2]->info==1 && fam_list[fam]->fam_member[j][4]==group){
              ind_2 = fam_list[fam]->fam_member[j][0];
              iterk = (fam_list[fam]->kinship_coeff).find(pair<int, int>(ind_1,ind_2));
              if(iterk == (fam_list[fam]->kinship_coeff).end()){
                printf("No IBD coefficient between individual %d and individual %d from family %d. Please check...\n\n", ind_1,ind_2,fam_list[fam]->fam_id);
                exit(1);
              }

              if(iterk->second>1e-8){ // i is related to j
                par_2 = 0;
                par_2+=missing_trio_list[mis_trio_cat_2]->typed_par1;
                par_2+=missing_trio_list[mis_trio_cat_2]->typed_par2;
                ind_2f = fam_list[fam]->fam_member[j][1];
                ind_2m = fam_list[fam]->fam_member[j][2];

                kin = iterk->second;
                iter2 = (fam_list[fam]->IBD2_coeff).find(pair<int, int>(ind_1,ind_2));
                iter1 = (fam_list[fam]->IBD1_coeff).find(pair<int, int>(ind_1,ind_2));
                IBD_2 = iter2->second;
                IBD_1 = iter1->second;

                if(par_1==0 && par_2==0){
                  if(IBD_1>1e-8){
                    for(int k=0;k<h_size-1;k++)
                      memset(work1[k],0,(h_size-1)*sizeof(double));
                      
                    for(int k=0;k<dim_1;k++)
                      for(int l=0;l<dim_2;l++)
                        for(int m=0;m<h_size;m++)
                          work1[k][l]+=result_off_one[off_one[mis_trio_cat_1]][k][m]*result_off_one[off_one[mis_trio_cat_2]][l][m]*freq[m];
                          
                    for(int k=0;k<dim_1;k++)
                      for(int l=0;l<dim_2;l++)
                        covMatrix[nb1+k][nb2+l]+=IBD_1*work1[k][l];        
                  }
                  if(IBD_2>1e-8){
                    for(int k=0;k<h_size-1;k++)
                      memset(work2[k],0,(h_size-1)*sizeof(double));
                      
                    for(int k=0;k<dim_1;k++)
                      for(int l=0;l<dim_2;l++)
                        for(int m=0;m<col;m++)
                          work2[k][l]+=result_off_two[off_two[mis_trio_cat_1]][k][m]*result_off_two[off_two[mis_trio_cat_2]][l][m]*genofreq[m];
                          
                    for(int k=0;k<dim_1;k++)
                      for(int l=0;l<dim_2;l++)
                        covMatrix[nb1+k][nb2+l]+=IBD_2*work2[k][l];      
                  }
                }
                else if(par_1==0 && (par_2==1 || par_2==2)){
                  int order = fam_list[fam]->mis_par_order[j];
                  int ind_2par1 = ind_2f;
                  int ind_2par2 = ind_2m;
                  if(order==1){
                    ind_2par1 = ind_2m;
                    ind_2par2 = ind_2f;
                  }  
                  int diff = 0;
                  int index = relation(fam,ind_1,ind_2,ind_2par1,ind_2par2,kin,IBD_2,IBD_1,diff);
                  if(index==4){
                    for(int k=0;k<h_size-1;k++)
                      memset(work2[k],0,(h_size-1)*sizeof(double));

                    for(int k=0;k<dim_1;k++)
                      for(int l=0;l<dim_2;l++)
                        for(int m=0;m<col;m++)
                          work2[k][l]+=result_off_two[off_two[mis_trio_cat_1]][k][m]*result_par1_two[par1_two[mis_trio_cat_2]][l][m]*genofreq[m];
                          
                    for(int k=0;k<dim_1;k++)
                      for(int l=0;l<dim_2;l++)
                        covMatrix[nb1+k][nb2+l]+=work2[k][l];
                  }
                  else if(index==5){
                    for(int k=0;k<h_size-1;k++)
                      memset(work2[k],0,(h_size-1)*sizeof(double));
                      
                    for(int k=0;k<dim_1;k++)
                      for(int l=0;l<dim_2;l++)
                        for(int m=0;m<col;m++)
                          work2[k][l]+=result_off_two[off_two[mis_trio_cat_1]][k][m]*result_par2_two[par2_two[mis_trio_cat_2]][l][m]*genofreq[m];
                          
                    for(int k=0;k<dim_1;k++)
                      for(int l=0;l<dim_2;l++)
                        covMatrix[nb1+k][nb2+l]+=work2[k][l];  
                  }
                  else if(index==1){
                    iter2 = (fam_list[fam]->IBD2_coeff).find(pair<int, int>(ind_1,ind_2par1));
                    iter1 = (fam_list[fam]->IBD1_coeff).find(pair<int, int>(ind_1,ind_2par1));
                    IBD_2 = iter2->second;
                    IBD_1 = iter1->second;
                    
                    if(IBD_1>1e-8){
                      for(int k=0;k<h_size-1;k++)
                        memset(work1[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<h_size;m++)
                            work1[k][l]+=result_off_one[off_one[mis_trio_cat_1]][k][m]*result_par1_one[par1_one[mis_trio_cat_2]][l][m]*freq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=IBD_1*work1[k][l];        
                    }
                    if(IBD_2>1e-8){
                      for(int k=0;k<h_size-1;k++)
                        memset(work2[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<col;m++)
                            work2[k][l]+=result_off_two[off_two[mis_trio_cat_1]][k][m]*result_par1_two[par1_two[mis_trio_cat_2]][l][m]*genofreq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=IBD_2*work2[k][l];      
                    }                  
                  }
                  else if(index==2){
                    iter2 = (fam_list[fam]->IBD2_coeff).find(pair<int, int>(ind_1,ind_2par2));
                    iter1 = (fam_list[fam]->IBD1_coeff).find(pair<int, int>(ind_1,ind_2par2));
                    IBD_2 = iter2->second;
                    IBD_1 = iter1->second;
                    
                    if(IBD_1>1e-8){
                      for(int k=0;k<h_size-1;k++)
                        memset(work1[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<h_size;m++)
                            work1[k][l]+=result_off_one[off_one[mis_trio_cat_1]][k][m]*result_par2_one[par2_one[mis_trio_cat_2]][l][m]*freq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=IBD_1*work1[k][l];        
                    }
                    if(IBD_2>1e-8){
                      for(int k=0;k<h_size-1;k++)
                        memset(work2[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<col;m++)
                            work2[k][l]+=result_off_two[off_two[mis_trio_cat_1]][k][m]*result_par2_two[par2_two[mis_trio_cat_2]][l][m]*genofreq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=IBD_2*work2[k][l];      
                    }   
                  }
                  else if(index==3){
                    if(IBD_1>1e-8){
                      for(int k=0;k<h_size-1;k++)
                        memset(work1[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<h_size;m++)
                            work1[k][l]+=result_off_one[off_one[mis_trio_cat_1]][k][m]*result_off_one[off_one[mis_trio_cat_2]][l][m]*freq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=IBD_1*work1[k][l];        
                    }
                    if(IBD_2>1e-8){
                      for(int k=0;k<h_size-1;k++)
                        memset(work2[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<col;m++)
                            work2[k][l]+=result_off_two[off_two[mis_trio_cat_1]][k][m]*result_off_two[off_two[mis_trio_cat_2]][l][m]*genofreq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=IBD_2*work2[k][l];      
                    }   
                  }
                  else if(index==6){  // can not be sibs
                    if(diff==1)
                      printf("please check relationships\n");
                    
                    iter1 = (fam_list[fam]->IBD1_coeff).find(pair<int, int>(ind_1,ind_2par1));
                    IBD_1 = iter1->second;
                    
                    if(IBD_1>1e-8){
                      for(int k=0;k<h_size-1;k++)
                        memset(work1[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<h_size;m++)
                            work1[k][l]+=result_off_one[off_one[mis_trio_cat_1]][k][m]*result_par1_one[par1_one[mis_trio_cat_2]][l][m]*freq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=IBD_1*work1[k][l];        
                    }  
                    
                    iter1 = (fam_list[fam]->IBD1_coeff).find(pair<int, int>(ind_1,ind_2par2));
                    IBD_1 = iter1->second;
                    
                    if(IBD_1>1e-8){
                      for(int k=0;k<h_size-1;k++)
                        memset(work1[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<h_size;m++)
                            work1[k][l]+=result_off_one[off_one[mis_trio_cat_1]][k][m]*result_par2_one[par2_one[mis_trio_cat_2]][l][m]*freq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=IBD_1*work1[k][l];        
                    }               
                  } 
                  // index 7 & 8                        
                }
                else if((par_1==1 || par_1==2) && par_2==0){
                  int order = fam_list[fam]->mis_par_order[i];
                  int ind_1par1 = ind_1f;
                  int ind_1par2 = ind_1m;
                  if(order==1){
                    ind_1par1 = ind_1m;
                    ind_1par2 = ind_1f;
                  }  
                  int diff = 0;
                  int index = relation(fam,ind_2,ind_1,ind_1par1,ind_1par2,kin,IBD_2,IBD_1,diff);
                  if(index==4){
                    for(int k=0;k<h_size-1;k++)
                      memset(work2[k],0,(h_size-1)*sizeof(double));
                      
                    for(int k=0;k<dim_1;k++)
                      for(int l=0;l<dim_2;l++)
                        for(int m=0;m<col;m++)
                          work2[k][l]+=result_par1_two[par1_two[mis_trio_cat_1]][k][m]*result_off_two[off_two[mis_trio_cat_2]][l][m]*genofreq[m];
                          
                    for(int k=0;k<dim_1;k++)
                      for(int l=0;l<dim_2;l++)
                        covMatrix[nb1+k][nb2+l]+=work2[k][l];           
                  }
                  else if(index==5){
                    for(int k=0;k<h_size-1;k++)
                      memset(work2[k],0,(h_size-1)*sizeof(double));
                      
                    for(int k=0;k<dim_1;k++)
                      for(int l=0;l<dim_2;l++)
                        for(int m=0;m<col;m++)
                          work2[k][l]+=result_par2_two[par2_two[mis_trio_cat_1]][k][m]*result_off_two[off_two[mis_trio_cat_2]][l][m]*genofreq[m];
                          
                    for(int k=0;k<dim_1;k++)
                      for(int l=0;l<dim_2;l++)
                        covMatrix[nb1+k][nb2+l]+=work2[k][l];  
                  }
                  else if(index==1){
                    iter2 = (fam_list[fam]->IBD2_coeff).find(pair<int, int>(ind_1par1,ind_2));
                    iter1 = (fam_list[fam]->IBD1_coeff).find(pair<int, int>(ind_1par1,ind_2));
                    IBD_2 = iter2->second;
                    IBD_1 = iter1->second;
                    
                    if(IBD_1>1e-8){
                      for(int k=0;k<h_size-1;k++)
                        memset(work1[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<h_size;m++)
                            work1[k][l]+=result_par1_one[par1_one[mis_trio_cat_1]][k][m]*result_off_one[off_one[mis_trio_cat_2]][l][m]*freq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=IBD_1*work1[k][l];        
                    }
                    if(IBD_2>1e-8){
                      for(int k=0;k<h_size-1;k++)
                        memset(work2[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<col;m++)
                            work2[k][l]+=result_par1_two[par1_two[mis_trio_cat_1]][k][m]*result_off_two[off_two[mis_trio_cat_2]][l][m]*genofreq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=IBD_2*work2[k][l];      
                    }                  
                  }
                  else if(index==2){
                    iter2 = (fam_list[fam]->IBD2_coeff).find(pair<int, int>(ind_1par2,ind_2));
                    iter1 = (fam_list[fam]->IBD1_coeff).find(pair<int, int>(ind_1par2,ind_2));
                    IBD_2 = iter2->second;
                    IBD_1 = iter1->second;
                    
                    if(IBD_1>1e-8){
                      for(int k=0;k<h_size-1;k++)
                        memset(work1[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<h_size;m++)
                            work1[k][l]+=result_par2_one[par2_one[mis_trio_cat_1]][k][m]*result_off_one[off_one[mis_trio_cat_2]][l][m]*freq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=IBD_1*work1[k][l];        
                    }
                    if(IBD_2>1e-8){
                      for(int k=0;k<h_size-1;k++)
                        memset(work2[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<col;m++)
                            work2[k][l]+=result_par2_two[par2_two[mis_trio_cat_1]][k][m]*result_off_two[off_two[mis_trio_cat_2]][l][m]*genofreq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=IBD_2*work2[k][l];      
                    }   
                  }
                  else if(index==3){
                    if(IBD_1>1e-8){
                      for(int k=0;k<h_size-1;k++)
                        memset(work1[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<h_size;m++)
                            work1[k][l]+=result_off_one[off_one[mis_trio_cat_1]][k][m]*result_off_one[off_one[mis_trio_cat_2]][l][m]*freq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=IBD_1*work1[k][l];        
                    }
                    if(IBD_2>1e-8){
                      for(int k=0;k<h_size-1;k++)
                        memset(work2[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<col;m++)
                            work2[k][l]+=result_off_two[off_two[mis_trio_cat_1]][k][m]*result_off_two[off_two[mis_trio_cat_2]][l][m]*genofreq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=IBD_2*work2[k][l];      
                    }   
                  }
                  else if(index==6){  // can not be sibs
                    if(diff==1)
                      printf("please check relationships\n");
                    
                    iter1 = (fam_list[fam]->IBD1_coeff).find(pair<int, int>(ind_1par1,ind_2));
                    IBD_1 = iter1->second;
                    
                    if(IBD_1>1e-8){
                      for(int k=0;k<h_size-1;k++)
                        memset(work1[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<h_size;m++)
                            work1[k][l]+=result_par1_one[par1_one[mis_trio_cat_1]][k][m]*result_off_one[off_one[mis_trio_cat_2]][l][m]*freq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=IBD_1*work1[k][l];        
                    }  
                    
                    iter1 = (fam_list[fam]->IBD1_coeff).find(pair<int, int>(ind_1par2,ind_2));
                    IBD_1 = iter1->second;
                    
                    if(IBD_1>1e-8){
                      for(int k=0;k<h_size-1;k++)
                        memset(work1[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<h_size;m++)
                            work1[k][l]+=result_par2_one[par2_one[mis_trio_cat_1]][k][m]*result_off_one[off_one[mis_trio_cat_2]][l][m]*freq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=IBD_1*work1[k][l];        
                    }               
                  } 
                  // index 7 & 8                        
                }
                else if((par_1==1 || par_1==2) && (par_2==1 || par_2==2)){
                  int order_1 = fam_list[fam]->mis_par_order[i];
                  int ind_1par1 = ind_1f;
                  int ind_1par2 = ind_1m;
                  if(order_1==1){
                    ind_1par1 = ind_1m;
                    ind_1par2 = ind_1f;
                  }
                  int order_2 = fam_list[fam]->mis_par_order[j];
                  int ind_2par1 = ind_2f;
                  int ind_2par2 = ind_2m;
                  if(order_2==1){
                    ind_2par1 = ind_2m;
                    ind_2par2 = ind_2f;
                  }                    
                  int diff = 0;
                  int index = relation(fam,ind_1,ind_1par1,ind_1par2,ind_2,ind_2par1,ind_2par2,kin,IBD_2,IBD_1,diff);
                  if(index==13){
                    for(int k=0;k<h_size-1;k++)
                      memset(work2[k],0,(h_size-1)*sizeof(double));
                      
                    for(int k=0;k<dim_1;k++)
                      for(int l=0;l<dim_2;l++)
                        for(int m=0;m<col;m++)
                          work2[k][l]+=result_off_two[off_two[mis_trio_cat_1]][k][m]*result_par1_two[par1_two[mis_trio_cat_2]][l][m]*genofreq[m];
                          
                    for(int k=0;k<dim_1;k++)
                      for(int l=0;l<dim_2;l++)
                        covMatrix[nb1+k][nb2+l]+=work2[k][l];              
                  }
                  else if(index==14){
                    for(int k=0;k<h_size-1;k++)
                      memset(work2[k],0,(h_size-1)*sizeof(double));
                      
                    for(int k=0;k<dim_1;k++)
                      for(int l=0;l<dim_2;l++)
                        for(int m=0;m<col;m++)
                          work2[k][l]+=result_off_two[off_two[mis_trio_cat_1]][k][m]*result_par2_two[par2_two[mis_trio_cat_2]][l][m]*genofreq[m];
                          
                    for(int k=0;k<dim_1;k++)
                      for(int l=0;l<dim_2;l++)
                        covMatrix[nb1+k][nb2+l]+=work2[k][l];   
                  }
                  else if(index==15){
                    for(int k=0;k<h_size-1;k++)
                      memset(work2[k],0,(h_size-1)*sizeof(double));
                      
                    for(int k=0;k<dim_1;k++)
                      for(int l=0;l<dim_2;l++)
                        for(int m=0;m<col;m++)
                          work2[k][l]+=result_par1_two[par1_two[mis_trio_cat_1]][k][m]*result_off_two[off_two[mis_trio_cat_2]][l][m]*genofreq[m];
                          
                    for(int k=0;k<dim_1;k++)
                      for(int l=0;l<dim_2;l++)
                        covMatrix[nb1+k][nb2+l]+=work2[k][l];              
                  }
                  else if(index==16){
                    for(int k=0;k<h_size-1;k++)
                      memset(work2[k],0,(h_size-1)*sizeof(double));
                      
                    for(int k=0;k<dim_1;k++)
                      for(int l=0;l<dim_2;l++)
                        for(int m=0;m<col;m++)
                          work2[k][l]+=result_par2_two[par2_two[mis_trio_cat_1]][k][m]*result_off_two[off_two[mis_trio_cat_2]][l][m]*genofreq[m];
                          
                    for(int k=0;k<dim_1;k++)
                      for(int l=0;l<dim_2;l++)
                        covMatrix[nb1+k][nb2+l]+=work2[k][l];   
                  }
                  else if(index==17){
                    for(int k=0;k<h_size-1;k++)
                      memset(work2[k],0,(h_size-1)*sizeof(double));
                      
                    for(int k=0;k<dim_1;k++)
                      for(int l=0;l<dim_2;l++)
                        for(int m=0;m<col;m++)
                          for(int n=0;n<col;n++)
                            work2[k][l]+=result_parents_two[parents_two[mis_trio_cat_1]][k][m][n]*result_parents_two[parents_two[mis_trio_cat_2]][l][m][n]*genofreq[m]*genofreq[n];
                          
                    for(int k=0;k<dim_1;k++)
                      for(int l=0;l<dim_2;l++)
                        covMatrix[nb1+k][nb2+l]+=work2[k][l];     
                  }
                  else if(index==18){
                    if(ind_1par1==ind_2par1){
                      for(int k=0;k<h_size-1;k++)
                        memset(work2[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<col;m++)
                            work2[k][l]+=result_par1_two[par1_two[mis_trio_cat_1]][k][m]*result_par1_two[par1_two[mis_trio_cat_2]][l][m]*genofreq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=work2[k][l];                                          
                    }
                    else if(ind_1par1==ind_2par2){
                      for(int k=0;k<h_size-1;k++)
                        memset(work2[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<col;m++)
                            work2[k][l]+=result_par1_two[par1_two[mis_trio_cat_1]][k][m]*result_par2_two[par2_two[mis_trio_cat_2]][l][m]*genofreq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=work2[k][l];     
                    }
                    else if(ind_1par2==ind_2par1){
                      for(int k=0;k<h_size-1;k++)
                        memset(work2[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<col;m++)
                            work2[k][l]+=result_par2_two[par2_two[mis_trio_cat_1]][k][m]*result_par1_two[par1_two[mis_trio_cat_2]][l][m]*genofreq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=work2[k][l];      
                    }
                    else if(ind_1par2==ind_2par2){
                      for(int k=0;k<h_size-1;k++)
                        memset(work2[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<col;m++)
                            work2[k][l]+=result_par2_two[par2_two[mis_trio_cat_1]][k][m]*result_par2_two[par2_two[mis_trio_cat_2]][l][m]*genofreq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=work2[k][l];   
                    }     
                  }
                  else if(index==1){
                    iter2 = (fam_list[fam]->IBD2_coeff).find(pair<int, int>(ind_1par1,ind_2par1));
                    iter1 = (fam_list[fam]->IBD1_coeff).find(pair<int, int>(ind_1par1,ind_2par1));
                    IBD_2 = iter2->second;
                    IBD_1 = iter1->second;
                    
                    if(IBD_1>1e-8){
                      for(int k=0;k<h_size-1;k++)
                        memset(work1[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<h_size;m++)
                            work1[k][l]+=result_par1_one[par1_one[mis_trio_cat_1]][k][m]*result_par1_one[par1_one[mis_trio_cat_2]][l][m]*freq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=IBD_1*work1[k][l];        
                    }
                    if(IBD_2>1e-8){
                      for(int k=0;k<h_size-1;k++)
                        memset(work2[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<col;m++)
                            work2[k][l]+=result_par1_two[par1_two[mis_trio_cat_1]][k][m]*result_par1_two[par1_two[mis_trio_cat_2]][l][m]*genofreq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=IBD_2*work2[k][l];      
                    }   
                  }
                  else if(index==2){
                    iter2 = (fam_list[fam]->IBD2_coeff).find(pair<int, int>(ind_1par1,ind_2par2));
                    iter1 = (fam_list[fam]->IBD1_coeff).find(pair<int, int>(ind_1par1,ind_2par2));
                    IBD_2 = iter2->second;
                    IBD_1 = iter1->second;
                    
                    if(IBD_1>1e-8){
                      for(int k=0;k<h_size-1;k++)
                        memset(work1[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<h_size;m++)
                            work1[k][l]+=result_par1_one[par1_one[mis_trio_cat_1]][k][m]*result_par2_one[par2_one[mis_trio_cat_2]][l][m]*freq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=IBD_1*work1[k][l];        
                    }
                    if(IBD_2>1e-8){
                      for(int k=0;k<h_size-1;k++)
                        memset(work2[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<col;m++)
                            work2[k][l]+=result_par1_two[par1_two[mis_trio_cat_1]][k][m]*result_par2_two[par2_two[mis_trio_cat_2]][l][m]*genofreq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=IBD_2*work2[k][l];      
                    }   
                  }
                  else if(index==3){
                    iter2 = (fam_list[fam]->IBD2_coeff).find(pair<int, int>(ind_1par2,ind_2par1));
                    iter1 = (fam_list[fam]->IBD1_coeff).find(pair<int, int>(ind_1par2,ind_2par1));
                    IBD_2 = iter2->second;
                    IBD_1 = iter1->second;
                    
                    if(IBD_1>1e-8){
                      for(int k=0;k<h_size-1;k++)
                        memset(work1[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<h_size;m++)
                            work1[k][l]+=result_par2_one[par2_one[mis_trio_cat_1]][k][m]*result_par1_one[par1_one[mis_trio_cat_2]][l][m]*freq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=IBD_1*work1[k][l];        
                    }
                    if(IBD_2>1e-8){
                      for(int k=0;k<h_size-1;k++)
                        memset(work2[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<col;m++)
                            work2[k][l]+=result_par2_two[par2_two[mis_trio_cat_1]][k][m]*result_par1_two[par1_two[mis_trio_cat_2]][l][m]*genofreq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=IBD_2*work2[k][l];      
                    }   
                  }
                  else if(index==4){
                    iter2 = (fam_list[fam]->IBD2_coeff).find(pair<int, int>(ind_1par2,ind_2par2));
                    iter1 = (fam_list[fam]->IBD1_coeff).find(pair<int, int>(ind_1par2,ind_2par2));
                    IBD_2 = iter2->second;
                    IBD_1 = iter1->second;
                    
                    if(IBD_1>1e-8){
                      for(int k=0;k<h_size-1;k++)
                        memset(work1[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<h_size;m++)
                            work1[k][l]+=result_par2_one[par2_one[mis_trio_cat_1]][k][m]*result_par2_one[par2_one[mis_trio_cat_2]][l][m]*freq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=IBD_1*work1[k][l];        
                    }
                    if(IBD_2>1e-8){
                      for(int k=0;k<h_size-1;k++)
                        memset(work2[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<col;m++)
                            work2[k][l]+=result_par2_two[par2_two[mis_trio_cat_1]][k][m]*result_par2_two[par2_two[mis_trio_cat_2]][l][m]*genofreq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=IBD_2*work2[k][l];      
                    }   
                  }
                  else if(index==5){
                    iter2 = (fam_list[fam]->IBD2_coeff).find(pair<int, int>(ind_1par1,ind_2));
                    iter1 = (fam_list[fam]->IBD1_coeff).find(pair<int, int>(ind_1par1,ind_2));
                    IBD_2 = iter2->second;
                    IBD_1 = iter1->second;
                       
                    if(IBD_1>1e-8){
                      for(int k=0;k<h_size-1;k++)
                        memset(work1[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<h_size;m++)
                            work1[k][l]+=result_par1_one[par1_one[mis_trio_cat_1]][k][m]*result_off_one[off_one[mis_trio_cat_2]][l][m]*freq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=IBD_1*work1[k][l];        
                    }
                    if(IBD_2>1e-8){
                      for(int k=0;k<h_size-1;k++)
                        memset(work2[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<col;m++)
                            work2[k][l]+=result_par1_two[par1_two[mis_trio_cat_1]][k][m]*result_off_two[off_two[mis_trio_cat_2]][l][m]*genofreq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=IBD_2*work2[k][l];      
                    }     
                  }
                  else if(index==6){
                    iter2 = (fam_list[fam]->IBD2_coeff).find(pair<int, int>(ind_1par2,ind_2));
                    iter1 = (fam_list[fam]->IBD1_coeff).find(pair<int, int>(ind_1par2,ind_2));
                    IBD_2 = iter2->second;
                    IBD_1 = iter1->second;
                       
                    if(IBD_1>1e-8){
                      for(int k=0;k<h_size-1;k++)
                        memset(work1[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<h_size;m++)
                            work1[k][l]+=result_par2_one[par2_one[mis_trio_cat_1]][k][m]*result_off_one[off_one[mis_trio_cat_2]][l][m]*freq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=IBD_1*work1[k][l];        
                    }
                    if(IBD_2>1e-8){
                      for(int k=0;k<h_size-1;k++)
                        memset(work2[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<col;m++)
                            work2[k][l]+=result_par2_two[par2_two[mis_trio_cat_1]][k][m]*result_off_two[off_two[mis_trio_cat_2]][l][m]*genofreq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=IBD_2*work2[k][l];      
                    }     
                  }
                  else if(index==7){
                    iter2 = (fam_list[fam]->IBD2_coeff).find(pair<int, int>(ind_1,ind_2par1));
                    iter1 = (fam_list[fam]->IBD1_coeff).find(pair<int, int>(ind_1,ind_2par1));
                    IBD_2 = iter2->second;
                    IBD_1 = iter1->second;
                       
                    if(IBD_1>1e-8){
                      for(int k=0;k<h_size-1;k++)
                        memset(work1[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<h_size;m++)
                            work1[k][l]+=result_off_one[off_one[mis_trio_cat_1]][k][m]*result_par1_one[par1_one[mis_trio_cat_2]][l][m]*freq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=IBD_1*work1[k][l];        
                    }
                    if(IBD_2>1e-8){
                      for(int k=0;k<h_size-1;k++)
                        memset(work2[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<col;m++)
                            work2[k][l]+=result_off_two[off_two[mis_trio_cat_1]][k][m]*result_par1_two[par1_two[mis_trio_cat_2]][l][m]*genofreq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=IBD_2*work2[k][l];      
                    }     
                  }
                  else if(index==8){
                    iter2 = (fam_list[fam]->IBD2_coeff).find(pair<int, int>(ind_1,ind_2par2));
                    iter1 = (fam_list[fam]->IBD1_coeff).find(pair<int, int>(ind_1,ind_2par2));
                    IBD_2 = iter2->second;
                    IBD_1 = iter1->second;
                       
                    if(IBD_1>1e-8){
                      for(int k=0;k<h_size-1;k++)
                        memset(work1[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<h_size;m++)
                            work1[k][l]+=result_off_one[off_one[mis_trio_cat_1]][k][m]*result_par2_one[par2_one[mis_trio_cat_2]][l][m]*freq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=IBD_1*work1[k][l];        
                    }
                    if(IBD_2>1e-8){
                      for(int k=0;k<h_size-1;k++)
                        memset(work2[k],0,(h_size-1)*sizeof(double));
                      
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          for(int m=0;m<col;m++)
                            work2[k][l]+=result_off_two[off_two[mis_trio_cat_1]][k][m]*result_par2_two[par2_two[mis_trio_cat_2]][l][m]*genofreq[m];
                          
                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=IBD_2*work2[k][l];      
                    }     
                  }     
                  else if(index==9){
                    if(diff==1){  // j and if are sibs
                      for(int k=0;k<h_size-1;k++)
                        memset(work2[k],0,(h_size-1)*sizeof(double));
                      
                      loc=0;
                      giter = geno_freq_table.begin();
                      while(giter!=geno_freq_table.end()){
                        int h1 = (giter->first).first;
                        int h2 = (giter->first).second;
                        if(h1==h2){
                          for(int k=0;k<dim_1;k++)
                            for(int l=0;l<dim_2;l++)
                              work2[k][l]+=result_par1_two[par1_two[mis_trio_cat_1]][k][loc]*result_parents_one[parents_one[mis_trio_cat_2]][l][hap_map[h1]][hap_map[h2]]*genofreq[loc];
                        }
                        else{
                          for(int k=0;k<dim_1;k++)
                            for(int l=0;l<dim_2;l++)
                              work2[k][l]+=result_par1_two[par1_two[mis_trio_cat_1]][k][loc]*(result_parents_one[parents_one[mis_trio_cat_2]][l][hap_map[h1]][hap_map[h2]]+result_parents_one[parents_one[mis_trio_cat_2]][l][hap_map[h2]][hap_map[h1]])*genofreq[loc]/2.0;
                        }
                        loc++;
                        giter++;
                      }

                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=work2[k][l];
                    }
                    else{
                      iter1 = (fam_list[fam]->IBD1_coeff).find(pair<int, int>(ind_1par1,ind_2par1));
                      IBD_1 = iter1->second;
                    
                      if(IBD_1>1e-8){
                        for(int k=0;k<h_size-1;k++)
                          memset(work1[k],0,(h_size-1)*sizeof(double));
                      
                        for(int k=0;k<dim_1;k++)
                          for(int l=0;l<dim_2;l++)
                            for(int m=0;m<h_size;m++)
                              work1[k][l]+=result_par1_one[par1_one[mis_trio_cat_1]][k][m]*result_par1_one[par1_one[mis_trio_cat_2]][l][m]*freq[m];
                          
                        for(int k=0;k<dim_1;k++)
                          for(int l=0;l<dim_2;l++)
                            covMatrix[nb1+k][nb2+l]+=IBD_1*work1[k][l];
                      }
                    
                      iter1 = (fam_list[fam]->IBD1_coeff).find(pair<int, int>(ind_1par1,ind_2par2));
                      IBD_1 = iter1->second;
                    
                      if(IBD_1>1e-8){
                        for(int k=0;k<h_size-1;k++)
                          memset(work1[k],0,(h_size-1)*sizeof(double));
                      
                        for(int k=0;k<dim_1;k++)
                          for(int l=0;l<dim_2;l++)
                            for(int m=0;m<h_size;m++)
                              work1[k][l]+=result_par1_one[par1_one[mis_trio_cat_1]][k][m]*result_par2_one[par2_one[mis_trio_cat_2]][l][m]*freq[m];
                          
                        for(int k=0;k<dim_1;k++)
                          for(int l=0;l<dim_2;l++)
                            covMatrix[nb1+k][nb2+l]+=IBD_1*work1[k][l];
                      }
                    }                          
                  }
                  else if(index==10){
                    if(diff==1){  // j and im are sibs
                      for(int k=0;k<h_size-1;k++)
                        memset(work2[k],0,(h_size-1)*sizeof(double));
                      
                      loc=0;
                      giter = geno_freq_table.begin();
                      while(giter!=geno_freq_table.end()){
                        int h1 = (giter->first).first;
                        int h2 = (giter->first).second;
                        if(h1==h2){
                          for(int k=0;k<dim_1;k++)
                            for(int l=0;l<dim_2;l++)
                              work2[k][l]+=result_par2_two[par2_two[mis_trio_cat_1]][k][loc]*result_parents_one[parents_one[mis_trio_cat_2]][l][hap_map[h1]][hap_map[h2]]*genofreq[loc];
                        }
                        else{
                          for(int k=0;k<dim_1;k++)
                            for(int l=0;l<dim_2;l++)
                              work2[k][l]+=result_par2_two[par2_two[mis_trio_cat_1]][k][loc]*(result_parents_one[parents_one[mis_trio_cat_2]][l][hap_map[h1]][hap_map[h2]]+result_parents_one[parents_one[mis_trio_cat_2]][l][hap_map[h2]][hap_map[h1]])*genofreq[loc]/2.0;
                        }
                        loc++;
                        giter++;
                      }

                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=work2[k][l];
                    }
                    else{
                      iter1 = (fam_list[fam]->IBD1_coeff).find(pair<int, int>(ind_1par2,ind_2par1));
                      IBD_1 = iter1->second;
                    
                      if(IBD_1>1e-8){
                        for(int k=0;k<h_size-1;k++)
                          memset(work1[k],0,(h_size-1)*sizeof(double));
                      
                        for(int k=0;k<dim_1;k++)
                          for(int l=0;l<dim_2;l++)
                            for(int m=0;m<h_size;m++)
                              work1[k][l]+=result_par2_one[par2_one[mis_trio_cat_1]][k][m]*result_par1_one[par1_one[mis_trio_cat_2]][l][m]*freq[m];
                          
                        for(int k=0;k<dim_1;k++)
                          for(int l=0;l<dim_2;l++)
                            covMatrix[nb1+k][nb2+l]+=IBD_1*work1[k][l];
                      }
                    
                      iter1 = (fam_list[fam]->IBD1_coeff).find(pair<int, int>(ind_1par2,ind_2par2));
                      IBD_1 = iter1->second;
                    
                      if(IBD_1>1e-8){
                        for(int k=0;k<h_size-1;k++)
                          memset(work1[k],0,(h_size-1)*sizeof(double));
                      
                        for(int k=0;k<dim_1;k++)
                          for(int l=0;l<dim_2;l++)
                            for(int m=0;m<h_size;m++)
                              work1[k][l]+=result_par2_one[par2_one[mis_trio_cat_1]][k][m]*result_par2_one[par2_one[mis_trio_cat_2]][l][m]*freq[m];
                          
                        for(int k=0;k<dim_1;k++)
                          for(int l=0;l<dim_2;l++)
                            covMatrix[nb1+k][nb2+l]+=IBD_1*work1[k][l];
                      }
                    }                          
                  }
                  else if(index==11){
                    if(diff==1){  // i and jf are sibs
                      for(int k=0;k<h_size-1;k++)
                        memset(work2[k],0,(h_size-1)*sizeof(double));
                      
                      loc=0;
                      giter = geno_freq_table.begin();
                      while(giter!=geno_freq_table.end()){
                        int h1 = (giter->first).first;
                        int h2 = (giter->first).second;
                        if(h1==h2){
                          for(int k=0;k<dim_1;k++)
                            for(int l=0;l<dim_2;l++)
                              work2[k][l]+=result_parents_one[parents_one[mis_trio_cat_1]][k][hap_map[h1]][hap_map[h2]]*result_par1_two[par1_two[mis_trio_cat_2]][l][loc]*genofreq[loc];
                        }
                        else{
                          for(int k=0;k<dim_1;k++)
                            for(int l=0;l<dim_2;l++)
                              work2[k][l]+=(result_parents_one[parents_one[mis_trio_cat_1]][k][hap_map[h1]][hap_map[h2]]+result_parents_one[parents_one[mis_trio_cat_1]][k][hap_map[h2]][hap_map[h1]])*result_par1_two[par1_two[mis_trio_cat_2]][l][loc]*genofreq[loc]/2.0;
                        }
                        loc++;
                        giter++;
                      }

                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=work2[k][l];
                    }
                    else{
                      iter1 = (fam_list[fam]->IBD1_coeff).find(pair<int, int>(ind_1par1,ind_2par1));
                      IBD_1 = iter1->second;
                    
                      if(IBD_1>1e-8){
                        for(int k=0;k<h_size-1;k++)
                          memset(work1[k],0,(h_size-1)*sizeof(double));
                      
                        for(int k=0;k<dim_1;k++)
                          for(int l=0;l<dim_2;l++)
                            for(int m=0;m<h_size;m++)
                              work1[k][l]+=result_par1_one[par1_one[mis_trio_cat_1]][k][m]*result_par1_one[par1_one[mis_trio_cat_2]][l][m]*freq[m];
                          
                        for(int k=0;k<dim_1;k++)
                          for(int l=0;l<dim_2;l++)
                            covMatrix[nb1+k][nb2+l]+=IBD_1*work1[k][l];
                      }
                    
                      iter1 = (fam_list[fam]->IBD1_coeff).find(pair<int, int>(ind_1par2,ind_2par1));
                      IBD_1 = iter1->second;
                    
                      if(IBD_1>1e-8){
                        for(int k=0;k<h_size-1;k++)
                          memset(work1[k],0,(h_size-1)*sizeof(double));
                      
                        for(int k=0;k<dim_1;k++)
                          for(int l=0;l<dim_2;l++)
                            for(int m=0;m<h_size;m++)
                              work1[k][l]+=result_par2_one[par2_one[mis_trio_cat_1]][k][m]*result_par1_one[par1_one[mis_trio_cat_2]][l][m]*freq[m];
                          
                        for(int k=0;k<dim_1;k++)
                          for(int l=0;l<dim_2;l++)
                            covMatrix[nb1+k][nb2+l]+=IBD_1*work1[k][l];
                      }
                    }                          
                  }
                  else if(index==12){
                    if(diff==1){  // i and jm are sibs
                      for(int k=0;k<h_size-1;k++)
                        memset(work2[k],0,(h_size-1)*sizeof(double));
                      
                      loc=0;
                      giter = geno_freq_table.begin();
                      while(giter!=geno_freq_table.end()){
                        int h1 = (giter->first).first;
                        int h2 = (giter->first).second;
                        if(h1==h2){
                          for(int k=0;k<dim_1;k++)
                            for(int l=0;l<dim_2;l++)
                              work2[k][l]+=result_parents_one[parents_one[mis_trio_cat_1]][k][hap_map[h1]][hap_map[h2]]*result_par2_two[par2_two[mis_trio_cat_2]][l][loc]*genofreq[loc];
                        }
                        else{
                          for(int k=0;k<dim_1;k++)
                            for(int l=0;l<dim_2;l++)
                              work2[k][l]+=(result_parents_one[parents_one[mis_trio_cat_1]][k][hap_map[h1]][hap_map[h2]]+result_parents_one[parents_one[mis_trio_cat_1]][k][hap_map[h2]][hap_map[h1]])*result_par2_two[par2_two[mis_trio_cat_2]][l][loc]*genofreq[loc]/2.0;
                        }
                        loc++;
                        giter++;
                      }

                      for(int k=0;k<dim_1;k++)
                        for(int l=0;l<dim_2;l++)
                          covMatrix[nb1+k][nb2+l]+=work2[k][l];
                    }
                    else{
                      iter1 = (fam_list[fam]->IBD1_coeff).find(pair<int, int>(ind_1par1,ind_2par2));
                      IBD_1 = iter1->second;
                    
                      if(IBD_1>1e-8){
                        for(int k=0;k<h_size-1;k++)
                          memset(work1[k],0,(h_size-1)*sizeof(double));
                      
                        for(int k=0;k<dim_1;k++)
                          for(int l=0;l<dim_2;l++)
                            for(int m=0;m<h_size;m++)
                              work1[k][l]+=result_par1_one[par1_one[mis_trio_cat_1]][k][m]*result_par2_one[par2_one[mis_trio_cat_2]][l][m]*freq[m];
                          
                        for(int k=0;k<dim_1;k++)
                          for(int l=0;l<dim_2;l++)
                            covMatrix[nb1+k][nb2+l]+=IBD_1*work1[k][l];
                      }
                    
                      iter1 = (fam_list[fam]->IBD1_coeff).find(pair<int, int>(ind_1par2,ind_2par2));
                      IBD_1 = iter1->second;
                    
                      if(IBD_1>1e-8){
                        for(int k=0;k<h_size-1;k++)
                          memset(work1[k],0,(h_size-1)*sizeof(double));
                      
                        for(int k=0;k<dim_1;k++)
                          for(int l=0;l<dim_2;l++)
                            for(int m=0;m<h_size;m++)
                              work1[k][l]+=result_par2_one[par2_one[mis_trio_cat_1]][k][m]*result_par2_one[par2_one[mis_trio_cat_2]][l][m]*freq[m];
                          
                        for(int k=0;k<dim_1;k++)
                          for(int l=0;l<dim_2;l++)
                            covMatrix[nb1+k][nb2+l]+=IBD_1*work1[k][l];
                      }
                    }                          
                  }   
                  // index 19 & 20
                }     
                for(int k=0;k<dim_1;k++)
                  for(int l=0;l<dim_2;l++)
                    covMatrix[nb2+l][nb1+k] = covMatrix[nb1+k][nb2+l];
              }  // i, j related
              nb2+=dim_2;
            }  // j is typed
          }  // complete j

          const vector<double> hap_num = geno_trio_list[trio_cat_1]->get_hap_num();
          for(int j=0;j<dim_1;j++){
            for(int k=0;k<h_size-1;k++)
              cholent[nb1+j][k] = missing_trio_list[mis_trio_cat_1]->derivative[j][k];
            cholent[nb1+j][h_size-1] = hap_num[j];
          }
          nb1+=dim_1;
        }  // i is typed
      }  // complete i

      if(cholesky(covMatrix,dim,cholent,h_size,chol,cholaug,0) != 1){
        //printf("cholesky decomposition of the cov matrix failed for family %d. \n", fam_list[fam]->fam_id);
        flag_pd = 0;
        //break;
      }

      if(flag_pd==1){
      double **s = self_colmultiply(cholaug,dim,1,h_size-1,1,h_size-1);
      double **t = self_colmultiply(cholaug,dim,1,h_size-1,h_size,h_size);

      for(int i=0;i<h_size-1;i++)
        for(int j=0;j<h_size-1;j++)
          denominator[i][j]+=s[i][j];

      for(int i=0;i<h_size-1;i++)
        numerator[i]+=t[i][0];

      clear_matrix(s,h_size-1,h_size-1);
      clear_matrix(t,h_size-1,1);
      }
      else{
      double var;
      double &determ = var;
      double **v = inverse(covMatrix,dim,determ);
      double **s = quadratic_product(cholent,dim,h_size,v,dim,dim,cholent,dim,h_size);

      for(int i=0;i<h_size-1;i++)
        for(int j=0;j<h_size-1;j++)
          denominator[i][j]+=s[i][j];

      for(int i=0;i<h_size-1;i++)
        numerator[i]+=s[i][h_size-1];

      clear_matrix(s,h_size,h_size);
      clear_matrix(v,dim,dim);
      }
    }
  }

  double **cholentmain = new double*[h_size-1];
  for(int i=0;i<h_size-1;i++){
    cholentmain[i] = new double[h_size];
    memset(cholentmain[i],0,h_size*sizeof(double));
  }

  double **cholmain = new double*[h_size-1];
  for(int i=0;i<h_size-1;i++){
    cholmain[i] = new double[h_size-1];
    memset(cholmain[i],0,(h_size-1)*sizeof(double));
  }

  double **cholaugmain = new double*[h_size-1];
  for(int i=0;i<h_size-1;i++){
    cholaugmain[i] = new double[h_size];
    memset(cholaugmain[i],0,h_size*sizeof(double));
  }

  for(int i=0;i<h_size-1;i++){
    cholentmain[i][i] = 1;
    cholentmain[i][h_size-1] = numerator[i];
  }

  if(cholesky(denominator,h_size-1,cholentmain,h_size,cholmain,cholaugmain,0) != 1){
    //printf("cholesky decomposition of the denominator matrix failed.\n");
    flag_pd = 0;
    //exit(1);
  }

  if(flag_pd==1){
  double **s = self_colmultiply(cholaugmain,h_size-1,1,h_size-1,h_size,h_size);
  iter=hap_freq_table.begin();
  double n=0;
  for(int i=0;i<h_size-1;i++){
    iter->second+=s[i][0];
    n+=iter->second;
    iter++;
  }
  iter->second = 1-n;

  clear_matrix(s,h_size-1,1);
  }
  else{
  double var;
  double &determ = var;
  double **v = inverse(denominator,h_size-1,determ);
  double **w = new double*[h_size-1];
  for(int i=0;i<h_size-1;i++){
    w[i] = new double[1];
    w[i][0] = numerator[i];
  }
  double **s = multiply(v,h_size-1,h_size-1,w,h_size-1,1);
  iter=hap_freq_table.begin();
  double n=0;
  for(int i=0;i<h_size-1;i++){
    iter->second+=s[i][0];
    n+=iter->second;
    iter++;
  }
  iter->second = 1-n;

  clear_matrix(s,h_size-1,1);
  clear_matrix(v,h_size-1,h_size-1);
  clear_matrix(w,h_size-1,1);
  }
  clear_matrix(cholentmain,h_size-1,h_size);
  clear_matrix(cholmain,h_size-1,h_size-1);
  clear_matrix(cholaugmain,h_size-1,h_size);

  update_geno_freq();
}


int compare_equ(int *a, int *b, int loci_num)
{
  for(int k=0;k<loci_num;k++)
  {
    if((a[k]==-1 && b[k]!=-1) || (b[k]==-1 && a[k]!=-1))
      return 0;
  }
  return 1;
}