#include "missingtrio.h"


Missingtrio::Missingtrio(Missing * M_offspring, Missing * M_par1, Missing * M_par2, int typed_o, int typed_par1, int typed_par2, const map<int, double> & hap_list, const vector<char *> & content, int size)
{
  if(size==0)
    return;

  this->M_offspring = M_offspring;
  this->M_par1 = M_par1;
  this->M_par2 = M_par2;

  this->typed_o = typed_o;
  this->typed_par1 = typed_par1;
  this->typed_par2 = typed_par2;
  
  this->size = size;  // # of SNPs
  this->count=0;
  this->ind_clus=0;
  this->ind_1df=0;

  set_vec();

  keep(hap_list,content);
  dim = sites.size();

  int row=hap_list.size();
  if(row==0)
    row=1<<size;
  int col = missing_pat_trio_genotype.size();

  if(dim>0)
  {
    weight = new double*[dim];
    for(int i=0;i<dim;i++)
    {
      weight[i] = new double[col];
      memset(weight[i],0,col*sizeof(double));
    }

    gfreq = new double[col];
    memset(gfreq,0,col*sizeof(double));

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

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

    if(typed_par1==1 || typed_par2==1)
    {
      deriv1 = new double*[dim];
      for(int i=0;i<dim;i++)
      {
        deriv1[i] = new double[row-1];
        memset(deriv1[i],0,(row-1)*sizeof(double));
      }

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

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


Missingtrio::Missingtrio(Missing * M_offspring, int typed_o, int typed_par1, int typed_par2, const map<int, double> & hap_list, const vector<char *> & content, int size)
{
  if(size==0)
    return;
  if(typed_par1==1)
    printf("Constructor of Missingtrio for founder fails. Please check\n");
  if(typed_par2==1)
    printf("Constructor of Missingtrio for founder fails. Please check\n");

  this->M_offspring = M_offspring;

  this->typed_o = typed_o;
  this->typed_par1 = typed_par1;
  this->typed_par2 = typed_par2;
  
  this->size = size;  // # of SNPs
  this->count=0;
  this->ind_clus=0;
  this->ind_1df=0;

  set_vec();

  keep(hap_list,content);
  dim = sites.size();

  int row=hap_list.size();
  if(row==0)
    row=1<<size;
  int col = missing_pat_trio_genotype.size();

  if(dim>0)
  {
    weight = new double*[dim];
    for(int i=0;i<dim;i++)
    {
      weight[i] = new double[col];
      memset(weight[i],0,col*sizeof(double));
    }

    gfreq = new double[col];
    memset(gfreq,0,col*sizeof(double));

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

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

void Missingtrio::set_vec()
{
  if(typed_par1==1 || typed_par2==1)
  {
    const vector<Genotype *> offspring = M_offspring->get_missing_pat_genotype();
    const vector<Genotype *> father = M_par1->get_missing_pat_genotype();
    const vector<Genotype *> mother = M_par2->get_missing_pat_genotype();

    for(int i=0;i<offspring.size();i++)
      for(int j=0;j<father.size();j++)
        for(int k=0;k<mother.size();k++)
        {
          Genotypetrio *trio_geno = new Genotypetrio(offspring[i],father[j],mother[k],typed_o,typed_par1,typed_par2,size);
          const vector<inheritance> iv = trio_geno->get_inh_vec();
          if(iv.size()>0)
            missing_pat_trio_genotype.push_back(trio_geno);
          else
            delete trio_geno;
        }
  }
  else if(typed_par1==0 && typed_par2==0)
  {
    const vector<Genotype *> offspring = M_offspring->get_missing_pat_genotype();

    for(int i=0;i<offspring.size();i++)
    {
      Genotypetrio *trio_geno = new Genotypetrio(offspring[i],typed_o,typed_par1,typed_par2,size);
      const vector<inheritance> iv = trio_geno->get_inh_vec();
      if(iv.size()>0)
        missing_pat_trio_genotype.push_back(trio_geno);
      else
        delete trio_geno;
    }
  }  // both parents are untyped
}


void Missingtrio::set_weight(map<int, double> & hap_freq_table, map<pair<int, int>, double> & geno_freq_table, map<pair<int, int>, int> & hap_pair_map, int **count)
{     
  int col = missing_pat_trio_genotype.size();

  for(int i=0;i<dim;i++)
    memset(weight[i],0,col*sizeof(double));
  memset(gfreq,0,col*sizeof(double));

  double value;
  int loc;

  if(typed_par1==1 && typed_par2==1)
  {
    for(int i=0;i<missing_pat_trio_genotype.size();i++)
    {
      const vector<inheritance> iv = missing_pat_trio_genotype[i]->get_inh_vec();
      if(iv.size()>0)
      {
        if(missing_pat_trio_genotype[i]->phase==1)
        {
          for(int j=0;j<iv.size();j++)
          {
            pair<int, int> hp_offspring = iv[j].hp_offspring;
            pair<int, int> hp_father = iv[j].hp_par1;
            pair<int, int> hp_mother = iv[j].hp_par2;

            gfreq[i]+=geno_freq_table[hp_father]*geno_freq_table[hp_mother]*transProb(hp_offspring,hp_father,hp_mother,hap_freq_table,geno_freq_table,1);
          }
          
          pair<int, int> hp_offspring = iv[0].hp_offspring;
          loc = hap_pair_map[hp_offspring];
          for(int j=0;j<dim;j++)
          {
            if(count[sites[j]][loc]>0)
              weight[j][i] = double(count[sites[j]][loc])/2.0;
          }
        }
        else
        {
          for(int j=0;j<iv.size();j++)
          {
            pair<int, int> hp_offspring = iv[j].hp_offspring;
            pair<int, int> hp_father = iv[j].hp_par1;
            pair<int, int> hp_mother = iv[j].hp_par2;

            value = geno_freq_table[hp_father]*geno_freq_table[hp_mother]*transProb(hp_offspring,hp_father,hp_mother,hap_freq_table,geno_freq_table,1);
            gfreq[i]+=value;

            loc = hap_pair_map[hp_offspring];
            for(int k=0;k<dim;k++)
            {
              if(count[sites[k]][loc]>0)
                weight[k][i]+=count[sites[k]][loc]*value;
            }
          }

          for(int j=0;j<dim;j++)
            weight[j][i]/=(2.0*gfreq[i]);
        }  // offspring haplotype phase is undetermined
      }  // trio hap config is possible
    }  // complete i
  }
  else if(typed_par1==1 && typed_par2==0)
  {
    for(int i=0;i<missing_pat_trio_genotype.size();i++)
    {
      const vector<inheritance> iv = missing_pat_trio_genotype[i]->get_inh_vec();
      if(iv.size()>0)
      {
        if(missing_pat_trio_genotype[i]->phase==1)
        {
          for(int j=0;j<iv.size();j++)
          {
            pair<int, int> hp_offspring = iv[j].hp_offspring;
            pair<int, int> hp_father = iv[j].hp_par1;

            gfreq[i]+=geno_freq_table[hp_father]*transProbonep(hp_offspring,hp_father,hap_freq_table,geno_freq_table,1);
          }
          
          pair<int, int> hp_offspring = iv[0].hp_offspring;
          loc = hap_pair_map[hp_offspring];
          for(int j=0;j<dim;j++)
          {
            if(count[sites[j]][loc]>0)
              weight[j][i] = double(count[sites[j]][loc])/2.0;
          }
        }
        else
        {
          for(int j=0;j<iv.size();j++)
          {
            pair<int, int> hp_offspring = iv[j].hp_offspring;
            pair<int, int> hp_father = iv[j].hp_par1;

            value = geno_freq_table[hp_father]*transProbonep(hp_offspring,hp_father,hap_freq_table,geno_freq_table,1);
            gfreq[i]+=value;

            loc = hap_pair_map[hp_offspring];
            for(int k=0;k<dim;k++)
            {
              if(count[sites[k]][loc]>0)
                weight[k][i]+=count[sites[k]][loc]*value;
            }
          }

          for(int j=0;j<dim;j++)
            weight[j][i]/=(2.0*gfreq[i]);
        }  // offspring haplotype phase is undetermined
      }  // trio hap config is possible
    }  // complete i
  }
  else if(typed_par1==0 && typed_par2==0)
  {
    for(int i=0;i<missing_pat_trio_genotype.size();i++)
    {
      const vector<inheritance> iv = missing_pat_trio_genotype[i]->get_inh_vec();
      if(iv.size()>0)
      {
        if(missing_pat_trio_genotype[i]->phase==1)
        {
          for(int j=0;j<iv.size();j++)
          {
            pair<int, int> hp_offspring = iv[j].hp_offspring;

            gfreq[i]+=geno_freq_table[hp_offspring];
          }
          
          pair<int, int> hp_offspring = iv[0].hp_offspring;
          loc = hap_pair_map[hp_offspring];
          for(int j=0;j<dim;j++)
          {
            if(count[sites[j]][loc]>0)
              weight[j][i] = double(count[sites[j]][loc])/2.0;
          }
        }
        else
        {
          for(int j=0;j<iv.size();j++)
          {
            pair<int, int> hp_offspring = iv[j].hp_offspring;

            gfreq[i]+=geno_freq_table[hp_offspring];

            loc = hap_pair_map[hp_offspring];
            for(int k=0;k<dim;k++)
            {
              if(count[sites[k]][loc]>0)
                weight[k][i]+=count[sites[k]][loc]*geno_freq_table[hp_offspring];
            }
          }

          for(int j=0;j<dim;j++)
            weight[j][i]/=(2.0*gfreq[i]);
        }  // offspring haplotype phase is undetermined
      }  // trio hap config is possible
    }  // complete i
  }
}


void Missingtrio::compute_centmatrix(map<int, double> & hap_freq_table, double *freq, map<pair<int, int>, double> & geno_freq_table){
     
  for(int i=0;i<dim;i++)
    memset(centmatrix[i],0,dim*sizeof(double));

  nonzero.clear();
  for(int i=0;i<missing_pat_trio_genotype.size();i++){
    for(int j=0;j<dim;j++){
      if(weight[j][i]>0){
        nonzero.push_back(i);
        break;
      }
    }
  }

  for(int i=0;i<nonzero.size();i++){
    for(int j=0;j<dim;j++){
      centmatrix[j][j]+=weight[j][nonzero[i]]*weight[j][nonzero[i]]*gfreq[nonzero[i]];
      for(int k=j+1;k<dim;k++)
        centmatrix[j][k]+=weight[j][nonzero[i]]*weight[k][nonzero[i]]*gfreq[nonzero[i]];
    }
  }

  for(int i=0;i<dim;i++){
    centmatrix[i][i]-=freq[sites[i]]*freq[sites[i]];
    for(int j=i+1;j<dim;j++){
      centmatrix[i][j]-=freq[sites[i]]*freq[sites[j]];
      centmatrix[j][i] = centmatrix[i][j];
    }
  }
}


void Missingtrio::compute_derivative(map<int, double> & hap_freq_table, map<pair<int, int>, double> & geno_freq_table, map<pair<int, int>, int> & hap_pair_map, int **count, double **rst, double ***hst){
     
  int row = hap_freq_table.size();

  for(int i=0;i<dim;i++)
    memset(derivative[i],0,(row-1)*sizeof(double));

  double value;
  int loc;
  int loc_f, loc_m;

  double *matrix1 = new double[row-1];
  memset(matrix1,0,(row-1)*sizeof(double));
  double **matrix2 = new double*[dim];
  for(int j=0;j<dim;j++){
    matrix2[j] = new double[row-1];
    memset(matrix2[j],0,(row-1)*sizeof(double));
  }

  if(typed_par1==1 && typed_par2==1)
  {
    for(int i=0;i<missing_pat_trio_genotype.size();i++)
    {
      const vector<inheritance> iv = missing_pat_trio_genotype[i]->get_inh_vec();
      if(iv.size()>0)
      {
        if(missing_pat_trio_genotype[i]->phase==0)
        {
          memset(matrix1,0,(row-1)*sizeof(double));
          for(int j=0;j<dim;j++)
            memset(matrix2[j],0,(row-1)*sizeof(double));

          for(int j=0;j<iv.size();j++)
          {
            pair<int, int> hp_offspring = iv[j].hp_offspring;
            pair<int, int> hp_father = iv[j].hp_par1;
            pair<int, int> hp_mother = iv[j].hp_par2;

            loc = hap_pair_map[hp_offspring];
            loc_f = hap_pair_map[hp_father];
            loc_m = hap_pair_map[hp_mother];
            for(int k=0;k<row-1;k++)
            {
              value = (rst[k][loc_f]*geno_freq_table[hp_mother]+geno_freq_table[hp_father]*rst[k][loc_m])*transProb(hp_offspring,hp_father,hp_mother,hap_freq_table,geno_freq_table,1);
              matrix1[k]+=value;
              for(int l=0;l<dim;l++)
              {
                if(count[sites[l]][loc]>0)
                  matrix2[l][k]+=count[sites[l]][loc]*value;
              }
            }
          }

          for(int j=0;j<dim;j++)
            for(int k=0;k<row-1;k++)
              derivative[j][k]-=(matrix2[j][k]/2.0-weight[j][i]*matrix1[k]);
        }
      }
    }
  }
  else if(typed_par1==1 && typed_par2==0)
  {
    for(int i=0;i<missing_pat_trio_genotype.size();i++)
    {
      const vector<inheritance> iv = missing_pat_trio_genotype[i]->get_inh_vec();
      if(iv.size()>0)
      {
        if(missing_pat_trio_genotype[i]->phase==0)
        {
          memset(matrix1,0,(row-1)*sizeof(double));
          for(int j=0;j<dim;j++)
            memset(matrix2[j],0,(row-1)*sizeof(double));

          for(int j=0;j<iv.size();j++)
          {
            pair<int, int> hp_offspring = iv[j].hp_offspring;
            pair<int, int> hp_father = iv[j].hp_par1;

            loc = hap_pair_map[hp_offspring];
            loc_f = hap_pair_map[hp_father];
            for(int k=0;k<row-1;k++)
            {
              value = rst[k][loc_f]*transProbonep(hp_offspring,hp_father,hap_freq_table,geno_freq_table,1)+geno_freq_table[hp_father]*hst[k][loc][loc_f];
              matrix1[k]+=value;
              for(int l=0;l<dim;l++)
              {
                if(count[sites[l]][loc]>0)
                  matrix2[l][k]+=count[sites[l]][loc]*value;
              }
            }
          }

          for(int j=0;j<dim;j++)
            for(int k=0;k<row-1;k++)
              derivative[j][k]-=(matrix2[j][k]/2.0-weight[j][i]*matrix1[k]);
        }
      }
    }
  }
  else if(typed_par1==0 && typed_par2==0)
  {
    for(int i=0;i<missing_pat_trio_genotype.size();i++)
    {
      const vector<inheritance> iv = missing_pat_trio_genotype[i]->get_inh_vec();
      if(iv.size()>0)
      {
        if(missing_pat_trio_genotype[i]->phase==0)
        {
          memset(matrix1,0,(row-1)*sizeof(double));
          for(int j=0;j<dim;j++)
            memset(matrix2[j],0,(row-1)*sizeof(double));

          for(int j=0;j<iv.size();j++)
          {
            pair<int, int> hp_offspring = iv[j].hp_offspring;

            loc = hap_pair_map[hp_offspring];
            for(int k=0;k<row-1;k++)
              matrix1[k]+=rst[k][loc];

            for(int k=0;k<dim;k++)
            {
              if(count[sites[k]][loc]>0)
              {
                for(int l=0;l<row-1;l++)
                  matrix2[k][l]+=count[sites[k]][loc]*rst[l][loc];
              }
            }
          }

          for(int j=0;j<dim;j++)
            for(int k=0;k<row-1;k++)
              derivative[j][k]-=(matrix2[j][k]/2.0-weight[j][i]*matrix1[k]);
        }
      }
    }
  }

  for(int i=0;i<dim;i++)
    derivative[i][sites[i]]+=1;

  delete[] matrix1;
  clear_matrix(matrix2,dim,row-1);
}


void Missingtrio::compute_deriv1to3(map<int, double> & hap_freq_table, map<pair<int, int>, double> & geno_freq_table, map<pair<int, int>, int> & hap_pair_map, int **count, double **rst, double ***drst)
{
  int row = hap_freq_table.size();

  for(int i=0;i<dim;i++)
  {
    memset(deriv1[i],0,(row-1)*sizeof(double));
    memset(deriv2[i],0,(row-1)*sizeof(double));
    memset(deriv3[i],0,(row-1)*sizeof(double));
  }

  double value;
  double value1, value2, value3;
  int loc;
  int loc_f, loc_m;

  double *matrix1 = new double[row-1];
  memset(matrix1,0,(row-1)*sizeof(double));
  double **matrix2 = new double*[dim];
  for(int j=0;j<dim;j++)
  {
    matrix2[j] = new double[row-1];
    memset(matrix2[j],0,(row-1)*sizeof(double));
  }
  double *matrix3 = new double[row-1];
  memset(matrix3,0,(row-1)*sizeof(double));
  double **matrix4 = new double*[dim];
  for(int j=0;j<dim;j++)
  {
    matrix4[j] = new double[row-1];
    memset(matrix4[j],0,(row-1)*sizeof(double));
  }
  double *matrix5 = new double[row-1];
  memset(matrix5,0,(row-1)*sizeof(double));
  double **matrix6 = new double*[dim];
  for(int j=0;j<dim;j++)
  {
    matrix6[j] = new double[row-1];
    memset(matrix6[j],0,(row-1)*sizeof(double));
  }

  if(typed_par1==1 && typed_par2==1)
  {
    for(int i=0;i<missing_pat_trio_genotype.size();i++)
    {
      const vector<inheritance> iv = missing_pat_trio_genotype[i]->get_inh_vec();
      if(iv.size()>0)
      {
        if(missing_pat_trio_genotype[i]->phase==0)
        {
          memset(matrix1,0,(row-1)*sizeof(double));
          for(int j=0;j<dim;j++)
            memset(matrix2[j],0,(row-1)*sizeof(double));
          memset(matrix3,0,(row-1)*sizeof(double));
          for(int j=0;j<dim;j++)
            memset(matrix4[j],0,(row-1)*sizeof(double));
          memset(matrix5,0,(row-1)*sizeof(double));
          for(int j=0;j<dim;j++)
            memset(matrix6[j],0,(row-1)*sizeof(double));

          for(int j=0;j<iv.size();j++)
          {
            pair<int, int> hp_offspring = iv[j].hp_offspring;
            pair<int, int> hp_father = iv[j].hp_par1;
            pair<int, int> hp_mother = iv[j].hp_par2;

            loc = hap_pair_map[hp_offspring];
            loc_f = hap_pair_map[hp_father];
            loc_m = hap_pair_map[hp_mother];

            value = transProb(hp_offspring,hp_father,hp_mother,hap_freq_table,geno_freq_table,1);
  
            for(int k=0;k<row-1;k++)
            {
              value1 = geno_freq_table[hp_father]*geno_freq_table[hp_mother]*value*rst[k][loc]/geno_freq_table[hp_offspring];
              value2 = rst[k][loc_f]*geno_freq_table[hp_mother]*value;
              value3 = geno_freq_table[hp_father]*rst[k][loc_m]*value;
              matrix1[k]+=value1;
              matrix3[k]+=value2;
              matrix5[k]+=value3;
              for(int l=0;l<dim;l++)
              {
                if(count[sites[l]][loc]>0)
                {
                  matrix2[l][k]+=count[sites[l]][loc]*value1;
                  matrix4[l][k]+=count[sites[l]][loc]*value2;
                  matrix6[l][k]+=count[sites[l]][loc]*value3;
                }
              }
            }
          }

          for(int j=0;j<dim;j++)
            for(int k=0;k<row-1;k++)
            {
              deriv1[j][k]-=(matrix2[j][k]/2.0-weight[j][i]*matrix1[k]);
              deriv2[j][k]-=(matrix4[j][k]/2.0-weight[j][i]*matrix3[k]);
              deriv3[j][k]-=(matrix6[j][k]/2.0-weight[j][i]*matrix5[k]);
            }
        }
      }
    }
  }
  else if(typed_par1==1 && typed_par2==0)
  {
    for(int i=0;i<missing_pat_trio_genotype.size();i++)
    {
      const vector<inheritance> iv = missing_pat_trio_genotype[i]->get_inh_vec();
      if(iv.size()>0)
      {
        if(missing_pat_trio_genotype[i]->phase==0)
        {
          memset(matrix1,0,(row-1)*sizeof(double));
          for(int j=0;j<dim;j++)
            memset(matrix2[j],0,(row-1)*sizeof(double));
          memset(matrix3,0,(row-1)*sizeof(double));
          for(int j=0;j<dim;j++)
            memset(matrix4[j],0,(row-1)*sizeof(double));
          memset(matrix5,0,(row-1)*sizeof(double));
          for(int j=0;j<dim;j++)
            memset(matrix6[j],0,(row-1)*sizeof(double));

          for(int j=0;j<iv.size();j++)
          {
            pair<int, int> hp_offspring = iv[j].hp_offspring;
            pair<int, int> hp_father = iv[j].hp_par1;

            loc = hap_pair_map[hp_offspring];
            loc_f = hap_pair_map[hp_father];

            value = transProbonep(hp_offspring,hp_father,hap_freq_table,geno_freq_table,1);
  
            for(int k=0;k<row-1;k++)
            {
              value1 = geno_freq_table[hp_father]*value*rst[k][loc]/geno_freq_table[hp_offspring];
              value2 = rst[k][loc_f]*value;
              value3 = geno_freq_table[hp_father]*value*drst[k][loc][loc_f];
              matrix1[k]+=value1;
              matrix3[k]+=value2;
              matrix5[k]+=value3;
              for(int l=0;l<dim;l++)
              {
                if(count[sites[l]][loc]>0)
                {
                  matrix2[l][k]+=count[sites[l]][loc]*value1;
                  matrix4[l][k]+=count[sites[l]][loc]*value2;
                  matrix6[l][k]+=count[sites[l]][loc]*value3;
                }
              }
            }
          }

          for(int j=0;j<dim;j++)
            for(int k=0;k<row-1;k++)
            {
              deriv1[j][k]-=(matrix2[j][k]/2.0-weight[j][i]*matrix1[k]);
              deriv2[j][k]-=(matrix4[j][k]/2.0-weight[j][i]*matrix3[k]);
              deriv3[j][k]-=(matrix6[j][k]/2.0-weight[j][i]*matrix5[k]);
            }
        }
      }
    }
  }
  
  delete[] matrix1;
  clear_matrix(matrix2,dim,row-1);
  delete[] matrix3;
  clear_matrix(matrix4,dim,row-1);
  delete[] matrix5;
  clear_matrix(matrix6,dim,row-1);
}


void Missingtrio::compute_derivative_cluster(map<int, double> & hap_freq_table, double **mul_clus, int col_clus){
     
  if(dim>0){
    ind_clus = 1;
    int row = hap_freq_table.size();
    derivative_clus = multiply(derivative,dim,row-1,mul_clus,row-1,col_clus);
  }
}


void Missingtrio::compute_deriv1to3_cluster(map<int, double> & hap_freq_table, double **mul_clus, int col_clus){
     
  if(dim>0){
    ind_clus = 1;
    int row = hap_freq_table.size();
    deriv1_clus = multiply(deriv1,dim,row-1,mul_clus,row-1,col_clus);
    deriv2_clus = multiply(deriv2,dim,row-1,mul_clus,row-1,col_clus);
    deriv3_clus = multiply(deriv3,dim,row-1,mul_clus,row-1,col_clus);
  }
}


void Missingtrio::compute_derivative_1df(map<int, double> & hap_freq_table, double **mul_1df, int col_1df){
     
  if(dim>0){
    ind_1df = 1;
    int row = hap_freq_table.size();
    derivative_1df = multiply(derivative,dim,row-1,mul_1df,row-1,col_1df);
  }
}


void Missingtrio::compute_deriv1to3_1df(map<int, double> & hap_freq_table, double **mul_1df, int col_1df){
     
  if(dim>0){
    ind_1df = 1;
    int row = hap_freq_table.size();
    deriv1_1df = multiply(deriv1,dim,row-1,mul_1df,row-1,col_1df);
    deriv2_1df = multiply(deriv2,dim,row-1,mul_1df,row-1,col_1df);
    deriv3_1df = multiply(deriv3,dim,row-1,mul_1df,row-1,col_1df);
  }
}


void Missingtrio::keep(const map<int, double> & hap_list, const vector<char *> & content)
{
  sites.clear();
/*
  if((typed_par1==1 || typed_par2==1) && typed_o==1)
  {
    int tot = 1<<size;

    int full_dim = hap_list.size();
    if(full_dim==0)
      full_dim = tot;

    if(M_offspring->dim==full_dim-1)
    {
      for(int i=0;i<full_dim-1;i++)
        sites.push_back(i);
    }
    else
    {
      vector<int> missing_site;

      if(typed_par1==1 && typed_par2==1)
      {
        int *a_o=M_offspring->get_gcode();
        int *a_f=M_par1->get_gcode();
        int *a_m=M_par2->get_gcode();

        for(int i=0;i<size;i++)
        {
          if(a_o[i]==-1 && a_f[i]==-1 && a_m[i]==-1)
            missing_site.push_back(i);
        }
      }
      else if(typed_par1==1 && typed_par2==0)
      {
        int *a_o=M_offspring->get_gcode();
        int *a_f=M_par1->get_gcode();

        for(int i=0;i<size;i++)
        {
          if(a_o[i]==-1 && a_f[i]==-1)
            missing_site.push_back(i);
        }
      }

      if(missing_site.size()==0)
      {
        for(int i=0;i<full_dim-1;i++)
          sites.push_back(i);
      }
      else if(missing_site.size()>0 && missing_site.size()<size)
      {
        map<int, int> hcom_hinc;

        if(hap_list.size()==0 || hap_list.size()==tot)
        {
          int *hcode = new int[size];
          memset(hcode,0,size*sizeof(int));

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

          for(int i=0;i<tot;i++)
          {
            int h_full = 0;
            int h_part = 0;
            int *a = (*vec)[i];
            for(int j=0;j<size;j++)
            {
              if(a[j]==1)
              {
                h_full+=(1<<j);
                if(!find_element(missing_site,j))
                  h_part+=(1<<j);
              }
            }

            hcom_hinc[h_full] = h_part;
          }

          delete[] hcode;

          for(int k=0;k<tot;k++)
            delete[] (*vec)[k];

          vec->clear();
          delete vec;
        }
        else
        {
          int hap_num = content.size();

          for(int i=0;i<hap_num;i++)
          {
            int h_full = 0;
            int h_part = 0;
            for(int j=0;j<size;j++)
            {
              if(content[i][j]=='1')
              {
                h_full+=(1<<j);
                if(!find_element(missing_site,j))
                  h_part+=(1<<j);
              }
            }

            hcom_hinc[h_full] = h_part;
          }
        }

        map<int, int> e_map;
        int loc=0;
        map<int, int>::iterator iter = hcom_hinc.begin();
        while(iter!=hcom_hinc.end())
        {
          int h_part = iter->second;
          if(e_map.find(h_part)==e_map.end())
            e_map[h_part] = loc;

          loc++;
          iter++;
        }

        loc=0;
        iter = e_map.begin();
        while(loc!=e_map.size()-1)
        {
          sites.push_back(iter->second);
          iter++;
          loc++;
        }
      }
      missing_site.clear();
    }
  }
  else
  {
    for(int i=0;i<M_offspring->dim;i++)
      sites.push_back(M_offspring->sites[i]);
  }
*/
  for(int i=0;i<M_offspring->dim;i++)
      sites.push_back(M_offspring->sites[i]);
}


vector<int*> *Missingtrio::complete_missing_haplotype(int *hcode, int index)
{
  int code = hcode[index];

  if(index==size-1)
  {
    vector<int*> *vec = new vector<int*>();
    
    if(code==-1)
    {
      int *a1 = new int[size];
      memset(a1,0,size*sizeof(int));
       
      a1[size-1] = code;
      vec->push_back(a1);
    }
    else if(code==0)
    {
      int *a1 = new int[size];
      int *a2 = new int[size];
      memset(a1,0,size*sizeof(int));
      memset(a2,0,size*sizeof(int));
      a1[size-1]=0;
      a2[size-1]=1;
      vec->push_back(a1);
      vec->push_back(a2);
    }

    return vec;
  }
  
  vector<int*> *vec = complete_missing_haplotype(hcode,index+1);

  int pos = index;
  
  if(code==-1)
  {
    for(int i=0;i<vec->size();i++)
      (*vec)[i][pos] = code;
    
    return vec;  
  }
  
  // else if code == 0
  int vsize = vec->size();
  for(int i=0;i<vsize;i++)
  {
    int *a1 = new int[size];
    memset(a1,0,size*sizeof(int));

    for(int j=pos+1;j<size;j++)
      a1[j]=(*vec)[i][j];
    
    a1[pos]=1;
    (*vec)[i][pos]=0;
    vec->push_back(a1);
  }
  
  return vec;
}


void Missingtrio::show_gstat()  // print trio genotype patten
{
  printf("\nOffspring has genotype patten: ");
  int *a=M_offspring->get_gcode();
  for(int i=0;i<size;i++)
  {
    if(a[i]!=-1)
      printf("%d", a[i]);
    else
      printf("*");
    if(i!=size-1)
      printf("-");
  }

  printf(" Parent 1 has genotype patten: ");
  if(typed_par1==1)
  {
    a=M_par1->get_gcode();
    for(int i=0;i<size;i++)
    {
      if(a[i]!=-1)
        printf("%d", a[i]);
      else
        printf("*");
      if(i!=size-1)
        printf("-");
    }
  }
  else
  {
    for(int i=0;i<size;i++)
    {
      printf("*");
      if(i!=size-1)
        printf("-");
    }
  }
  
  printf(" Parent 2 has genotype patten: ");
  if(typed_par2==1)
  {
    a=M_par2->get_gcode();
    for(int i=0;i<size;i++)
    {
      if(a[i]!=-1)
        printf("%d", a[i]);
      else
        printf("*");
      if(i!=size-1)
        printf("-");
    }
  }
  else
  {
    for(int i=0;i<size;i++)
    {
      printf("*");
      if(i!=size-1)
        printf("-");
    }
  }

  printf("\n");
}


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

  if(dim>0)
  {
    for(int i=0;i<dim;i++)
    {
      delete[] weight[i];
      delete[] centmatrix[i];
      delete[] derivative[i];
    }

    delete[] weight;
    delete[] gfreq;
    delete[] centmatrix;
    delete[] derivative;

    sites.clear();
    nonzero.clear();

    if(typed_par1==1 || typed_par2==1)
    {
      for(int i=0;i<dim;i++)
      {
        delete[] deriv1[i];
        delete[] deriv2[i];
        delete[] deriv3[i];
      }

      delete[] deriv1;
      delete[] deriv2;
      delete[] deriv3;
    }

    if(ind_clus==1)
    {
      for(int i=0;i<dim;i++)
        delete[] derivative_clus[i];
      delete[] derivative_clus;

      if(typed_par1==1 || typed_par2==1)
      {
        for(int i=0;i<dim;i++)
        {
          delete[] deriv1_clus[i];
          delete[] deriv2_clus[i];
          delete[] deriv3_clus[i];
        }

        delete[] deriv1_clus;
        delete[] deriv2_clus;
        delete[] deriv3_clus;
      }
    }

    if(ind_1df==1)
    {
      for(int i=0;i<dim;i++)
        delete[] derivative_1df[i];
      delete[] derivative_1df;

      if(typed_par1==1 || typed_par2==1)
      {
        for(int i=0;i<dim;i++)
        {
          delete[] deriv1_1df[i];
          delete[] deriv2_1df[i];
          delete[] deriv3_1df[i];
        }

        delete[] deriv1_1df;
        delete[] deriv2_1df;
        delete[] deriv3_1df;
      }
    }
  }
}


// calculate Z-h
void Missingtrio::expected_hap_num(map<int, double> & hap_freq_table, double *freq, map<pair<int, int>, double> & geno_freq_table, map<pair<int, int>, int> & hap_pair_map, int **count){
     
  for(int i=0;i<missing_pat_trio_genotype.size();i++)   
    missing_pat_trio_genotype[i]->expected_hap_num(hap_freq_table,freq,geno_freq_table,hap_pair_map,count,sites);
}


// calculate E(Z-h|event)
void Missingtrio::cond_par1_one(double **result, map<int, double> & hap_freq_table, map<pair<int, int>, double> & geno_freq_table, map<int, int> & hap_map){
// result is of dim*h_size
  int row = hap_freq_table.size();
  
  double *matrix1 = new double[row];
  memset(matrix1,0,row*sizeof(double));
  
  if(typed_par1==1 && typed_par2==1)
  {
    for(int i=0;i<missing_pat_trio_genotype.size();i++)
    {
      memset(matrix1,0,row*sizeof(double));
      const vector<inheritance> iv = missing_pat_trio_genotype[i]->get_inh_vec();
      if(iv.size()>0)
      {
        for(int j=0;j<iv.size();j++)
        {
          pair<int, int> hp_offspring = iv[j].hp_offspring;
          pair<int, int> hp_father = iv[j].hp_par1;
          pair<int, int> hp_mother = iv[j].hp_par2;
          
          if(hp_father.first==hp_father.second)
            matrix1[hap_map[hp_father.first]] += hap_freq_table[hp_father.first]*geno_freq_table[hp_mother]*transProb(hp_offspring,hp_father,hp_mother,hap_freq_table,geno_freq_table,1);
          else{
            matrix1[hap_map[hp_father.first]] += hap_freq_table[hp_father.second]*geno_freq_table[hp_mother]*transProb(hp_offspring,hp_father,hp_mother,hap_freq_table,geno_freq_table,1);
            matrix1[hap_map[hp_father.second]] += hap_freq_table[hp_father.first]*geno_freq_table[hp_mother]*transProb(hp_offspring,hp_father,hp_mother,hap_freq_table,geno_freq_table,1);   
          }         
        }
        
        const vector<double> hap_num = missing_pat_trio_genotype[i]->get_hap_num();
        for(int j=0;j<dim;j++){
          for(int k=0;k<row;k++)
            result[j][k]+=hap_num[j]*matrix1[k];        
        }        
      }  // trio hap config is possible
    }  // complete i
  }
  else if(typed_par1==1 && typed_par2==0)
  {
    for(int i=0;i<missing_pat_trio_genotype.size();i++)
    {
      memset(matrix1,0,row*sizeof(double));
      const vector<inheritance> iv = missing_pat_trio_genotype[i]->get_inh_vec();
      if(iv.size()>0)
      {
        for(int j=0;j<iv.size();j++)
        {
          pair<int, int> hp_offspring = iv[j].hp_offspring;
          pair<int, int> hp_father = iv[j].hp_par1;

          if(hp_father.first==hp_father.second)
            matrix1[hap_map[hp_father.first]] += hap_freq_table[hp_father.first]*transProbonep(hp_offspring,hp_father,hap_freq_table,geno_freq_table,1);
          else{
            matrix1[hap_map[hp_father.first]] += hap_freq_table[hp_father.second]*transProbonep(hp_offspring,hp_father,hap_freq_table,geno_freq_table,1);
            matrix1[hap_map[hp_father.second]] += hap_freq_table[hp_father.first]*transProbonep(hp_offspring,hp_father,hap_freq_table,geno_freq_table,1);   
          }
        }
        
        const vector<double> hap_num = missing_pat_trio_genotype[i]->get_hap_num();
        for(int j=0;j<dim;j++){
          for(int k=0;k<row;k++)
            result[j][k]+=hap_num[j]*matrix1[k];        
        }
      }  // trio hap config is possible
    }  // complete i
  }

  delete[] matrix1;
}


void Missingtrio::cond_par1_two(double **result, map<int, double> & hap_freq_table, map<pair<int, int>, double> & geno_freq_table, map<pair<int, int>, int> & hap_pair_map){
// result is of dim*g_size
  int col = geno_freq_table.size();
  
  double *matrix1 = new double[col];
  memset(matrix1,0,col*sizeof(double));
  
  if(typed_par1==1 && typed_par2==1)
  {
    for(int i=0;i<missing_pat_trio_genotype.size();i++)
    {
      memset(matrix1,0,col*sizeof(double));
      const vector<inheritance> iv = missing_pat_trio_genotype[i]->get_inh_vec();
      if(iv.size()>0)
      {
        for(int j=0;j<iv.size();j++)
        {
          pair<int, int> hp_offspring = iv[j].hp_offspring;
          pair<int, int> hp_father = iv[j].hp_par1;
          pair<int, int> hp_mother = iv[j].hp_par2;
          
          matrix1[hap_pair_map[hp_father]] += geno_freq_table[hp_mother]*transProb(hp_offspring,hp_father,hp_mother,hap_freq_table,geno_freq_table,1);         
        }
        
        const vector<double> hap_num = missing_pat_trio_genotype[i]->get_hap_num();
        for(int j=0;j<dim;j++){
          for(int k=0;k<col;k++)
            result[j][k]+=hap_num[j]*matrix1[k];        
        }        
      }  // trio hap config is possible
    }  // complete i
  }
  else if(typed_par1==1 && typed_par2==0)
  {
    for(int i=0;i<missing_pat_trio_genotype.size();i++)
    {
      memset(matrix1,0,col*sizeof(double));
      const vector<inheritance> iv = missing_pat_trio_genotype[i]->get_inh_vec();
      if(iv.size()>0)
      {
        for(int j=0;j<iv.size();j++)
        {
          pair<int, int> hp_offspring = iv[j].hp_offspring;
          pair<int, int> hp_father = iv[j].hp_par1;

          matrix1[hap_pair_map[hp_father]] += transProbonep(hp_offspring,hp_father,hap_freq_table,geno_freq_table,1);
        }
        
        const vector<double> hap_num = missing_pat_trio_genotype[i]->get_hap_num();
        for(int j=0;j<dim;j++){
          for(int k=0;k<col;k++)
            result[j][k]+=hap_num[j]*matrix1[k];        
        }
      }  // trio hap config is possible
    }  // complete i
  }

  delete[] matrix1;          
}


void Missingtrio::cond_par2_one(double **result, map<int, double> & hap_freq_table, map<pair<int, int>, double> & geno_freq_table, map<int, int> & hap_map){
// result is of dim*h_size
  int row = hap_freq_table.size();
  
  double *matrix1 = new double[row];
  memset(matrix1,0,row*sizeof(double));
  
  if(typed_par1==1 && typed_par2==1)
  {
    for(int i=0;i<missing_pat_trio_genotype.size();i++)
    {
      memset(matrix1,0,row*sizeof(double));
      const vector<inheritance> iv = missing_pat_trio_genotype[i]->get_inh_vec();
      if(iv.size()>0)
      {
        for(int j=0;j<iv.size();j++)
        {
          pair<int, int> hp_offspring = iv[j].hp_offspring;
          pair<int, int> hp_father = iv[j].hp_par1;
          pair<int, int> hp_mother = iv[j].hp_par2;
          
          if(hp_mother.first==hp_mother.second)
            matrix1[hap_map[hp_mother.first]] += geno_freq_table[hp_father]*hap_freq_table[hp_mother.first]*transProb(hp_offspring,hp_father,hp_mother,hap_freq_table,geno_freq_table,1);
          else{
            matrix1[hap_map[hp_mother.first]] += geno_freq_table[hp_father]*hap_freq_table[hp_mother.second]*transProb(hp_offspring,hp_father,hp_mother,hap_freq_table,geno_freq_table,1);
            matrix1[hap_map[hp_mother.second]] += geno_freq_table[hp_father]*hap_freq_table[hp_mother.first]*transProb(hp_offspring,hp_father,hp_mother,hap_freq_table,geno_freq_table,1);   
          }         
        }
        
        const vector<double> hap_num = missing_pat_trio_genotype[i]->get_hap_num();
        for(int j=0;j<dim;j++){
          for(int k=0;k<row;k++)
            result[j][k]+=hap_num[j]*matrix1[k];        
        }        
      }  // trio hap config is possible
    }  // complete i
  }
  else if(typed_par1==1 && typed_par2==0)
  {
    for(int i=0;i<missing_pat_trio_genotype.size();i++)
    {
      memset(matrix1,0,row*sizeof(double));
      const vector<inheritance> iv = missing_pat_trio_genotype[i]->get_inh_vec_complete();
      if(iv.size()>0)
      {
        for(int j=0;j<iv.size();j++)
        {
          pair<int, int> hp_offspring = iv[j].hp_offspring;
          pair<int, int> hp_father = iv[j].hp_par1;
          pair<int, int> hp_mother = iv[j].hp_par2;
          
          if(hp_mother.first==hp_mother.second)
            matrix1[hap_map[hp_mother.first]] += geno_freq_table[hp_father]*hap_freq_table[hp_mother.first]*transProb(hp_offspring,hp_father,hp_mother,hap_freq_table,geno_freq_table,1);
          else{
            matrix1[hap_map[hp_mother.first]] += geno_freq_table[hp_father]*hap_freq_table[hp_mother.second]*transProb(hp_offspring,hp_father,hp_mother,hap_freq_table,geno_freq_table,1);
            matrix1[hap_map[hp_mother.second]] += geno_freq_table[hp_father]*hap_freq_table[hp_mother.first]*transProb(hp_offspring,hp_father,hp_mother,hap_freq_table,geno_freq_table,1);   
          }
        }
        
        const vector<double> hap_num = missing_pat_trio_genotype[i]->get_hap_num();
        for(int j=0;j<dim;j++){
          for(int k=0;k<row;k++)
            result[j][k]+=hap_num[j]*matrix1[k];        
        }
      }  // trio hap config is possible
    }  // complete i
  }

  delete[] matrix1;          
}


void Missingtrio::cond_par2_two(double **result, map<int, double> & hap_freq_table, map<pair<int, int>, double> & geno_freq_table, map<pair<int, int>, int> & hap_pair_map){
// result is of dim*g_size
  int col = geno_freq_table.size();
  
  double *matrix1 = new double[col];
  memset(matrix1,0,col*sizeof(double));
  
  if(typed_par1==1 && typed_par2==1)
  {
    for(int i=0;i<missing_pat_trio_genotype.size();i++)
    {
      memset(matrix1,0,col*sizeof(double));
      const vector<inheritance> iv = missing_pat_trio_genotype[i]->get_inh_vec();
      if(iv.size()>0)
      {
        for(int j=0;j<iv.size();j++)
        {
          pair<int, int> hp_offspring = iv[j].hp_offspring;
          pair<int, int> hp_father = iv[j].hp_par1;
          pair<int, int> hp_mother = iv[j].hp_par2;
          
          matrix1[hap_pair_map[hp_mother]] += geno_freq_table[hp_father]*transProb(hp_offspring,hp_father,hp_mother,hap_freq_table,geno_freq_table,1);         
        }
        
        const vector<double> hap_num = missing_pat_trio_genotype[i]->get_hap_num();
        for(int j=0;j<dim;j++){
          for(int k=0;k<col;k++)
            result[j][k]+=hap_num[j]*matrix1[k];        
        }        
      }  // trio hap config is possible
    }  // complete i
  }
  else if(typed_par1==1 && typed_par2==0)
  {
    for(int i=0;i<missing_pat_trio_genotype.size();i++)
    {
      memset(matrix1,0,col*sizeof(double));
      const vector<inheritance> iv = missing_pat_trio_genotype[i]->get_inh_vec_complete();
      if(iv.size()>0)
      {
        for(int j=0;j<iv.size();j++)
        {
          pair<int, int> hp_offspring = iv[j].hp_offspring;
          pair<int, int> hp_father = iv[j].hp_par1;
          pair<int, int> hp_mother = iv[j].hp_par2;
          
          matrix1[hap_pair_map[hp_mother]] += geno_freq_table[hp_father]*transProb(hp_offspring,hp_father,hp_mother,hap_freq_table,geno_freq_table,1);
        }
        
        const vector<double> hap_num = missing_pat_trio_genotype[i]->get_hap_num();
        for(int j=0;j<dim;j++){
          for(int k=0;k<col;k++)
            result[j][k]+=hap_num[j]*matrix1[k];        
        }
      }  // trio hap config is possible
    }  // complete i
  }

  delete[] matrix1;
}


void Missingtrio::cond_off_one(double **result, map<int, double> & hap_freq_table, map<pair<int, int>, double> & geno_freq_table, map<int, int> & hap_map){
// result is of dim*h_size
  int row = hap_freq_table.size();
  
  double *matrix1 = new double[row];
  memset(matrix1,0,row*sizeof(double));
  
  if(typed_par1==1 && typed_par2==1)
  {
    for(int i=0;i<missing_pat_trio_genotype.size();i++)
    {
      memset(matrix1,0,row*sizeof(double));
      const vector<inheritance> iv = missing_pat_trio_genotype[i]->get_inh_vec();
      if(iv.size()>0)
      {
        for(int j=0;j<iv.size();j++)
        {
          pair<int, int> hp_offspring = iv[j].hp_offspring;
          pair<int, int> hp_father = iv[j].hp_par1;
          pair<int, int> hp_mother = iv[j].hp_par2;
          
          if(hp_offspring.first==hp_offspring.second)
            matrix1[hap_map[hp_offspring.first]] += geno_freq_table[hp_father]*geno_freq_table[hp_mother]*transProb(hp_offspring,hp_father,hp_mother,hap_freq_table,geno_freq_table,1)/hap_freq_table[hp_offspring.first];
          else{
            matrix1[hap_map[hp_offspring.first]] += geno_freq_table[hp_father]*geno_freq_table[hp_mother]*transProb(hp_offspring,hp_father,hp_mother,hap_freq_table,geno_freq_table,1)/(2.0*hap_freq_table[hp_offspring.first]);
            matrix1[hap_map[hp_offspring.second]] += geno_freq_table[hp_father]*geno_freq_table[hp_mother]*transProb(hp_offspring,hp_father,hp_mother,hap_freq_table,geno_freq_table,1)/(2.0*hap_freq_table[hp_offspring.second]);   
          }         
        }
        
        const vector<double> hap_num = missing_pat_trio_genotype[i]->get_hap_num();
        for(int j=0;j<dim;j++){
          for(int k=0;k<row;k++)
            result[j][k]+=hap_num[j]*matrix1[k];        
        }        
      }  // trio hap config is possible
    }  // complete i
  }
  else if(typed_par1==1 && typed_par2==0)
  {
    for(int i=0;i<missing_pat_trio_genotype.size();i++)
    {
      memset(matrix1,0,row*sizeof(double));
      const vector<inheritance> iv = missing_pat_trio_genotype[i]->get_inh_vec();
      if(iv.size()>0)
      {
        for(int j=0;j<iv.size();j++)
        {
          pair<int, int> hp_offspring = iv[j].hp_offspring;
          pair<int, int> hp_father = iv[j].hp_par1;

          if(hp_offspring.first==hp_offspring.second)
            matrix1[hap_map[hp_offspring.first]] += geno_freq_table[hp_father]*transProbonep(hp_offspring,hp_father,hap_freq_table,geno_freq_table,1)/hap_freq_table[hp_offspring.first];
          else{
            matrix1[hap_map[hp_offspring.first]] += geno_freq_table[hp_father]*transProbonep(hp_offspring,hp_father,hap_freq_table,geno_freq_table,1)/(2.0*hap_freq_table[hp_offspring.first]);
            matrix1[hap_map[hp_offspring.second]] += geno_freq_table[hp_father]*transProbonep(hp_offspring,hp_father,hap_freq_table,geno_freq_table,1)/(2.0*hap_freq_table[hp_offspring.second]);   
          }
        }
        
        const vector<double> hap_num = missing_pat_trio_genotype[i]->get_hap_num();
        for(int j=0;j<dim;j++){
          for(int k=0;k<row;k++)
            result[j][k]+=hap_num[j]*matrix1[k];        
        }
      }  // trio hap config is possible
    }  // complete i
  }
  else if(typed_par1==0 && typed_par2==0)
  {
    for(int i=0;i<missing_pat_trio_genotype.size();i++)
    {
      memset(matrix1,0,row*sizeof(double));
      const vector<inheritance> iv = missing_pat_trio_genotype[i]->get_inh_vec();
      if(iv.size()>0)
      {
        for(int j=0;j<iv.size();j++)
        {
          pair<int, int> hp_offspring = iv[j].hp_offspring;

          if(hp_offspring.first==hp_offspring.second)
            matrix1[hap_map[hp_offspring.first]] += hap_freq_table[hp_offspring.first];
          else{
            matrix1[hap_map[hp_offspring.first]] += hap_freq_table[hp_offspring.second];
            matrix1[hap_map[hp_offspring.second]] += hap_freq_table[hp_offspring.first];
          }
        }
        
        const vector<double> hap_num = missing_pat_trio_genotype[i]->get_hap_num();
        for(int j=0;j<dim;j++){
          for(int k=0;k<row;k++)
            result[j][k]+=hap_num[j]*matrix1[k];        
        }
      }  // trio hap config is possible
    }  // complete i
  }

  delete[] matrix1;
}


void Missingtrio::cond_off_two(double **result, map<int, double> & hap_freq_table, map<pair<int, int>, double> & geno_freq_table, map<pair<int, int>, int> & hap_pair_map){
// result is of dim*g_size
  int col = geno_freq_table.size();
  
  double *matrix1 = new double[col];
  memset(matrix1,0,col*sizeof(double));
  
  if(typed_par1==1 && typed_par2==1)
  {
    for(int i=0;i<missing_pat_trio_genotype.size();i++)
    {
      memset(matrix1,0,col*sizeof(double));
      const vector<inheritance> iv = missing_pat_trio_genotype[i]->get_inh_vec();
      if(iv.size()>0)
      {
        for(int j=0;j<iv.size();j++)
        {
          pair<int, int> hp_offspring = iv[j].hp_offspring;
          pair<int, int> hp_father = iv[j].hp_par1;
          pair<int, int> hp_mother = iv[j].hp_par2;
          
          matrix1[hap_pair_map[hp_offspring]] += geno_freq_table[hp_father]*geno_freq_table[hp_mother]*transProb(hp_offspring,hp_father,hp_mother,hap_freq_table,geno_freq_table,1)/geno_freq_table[hp_offspring];         
        }
        
        const vector<double> hap_num = missing_pat_trio_genotype[i]->get_hap_num();
        for(int j=0;j<dim;j++){
          for(int k=0;k<col;k++)
            result[j][k]+=hap_num[j]*matrix1[k];        
        }        
      }  // trio hap config is possible
    }  // complete i
  }
  else if(typed_par1==1 && typed_par2==0)
  {
    for(int i=0;i<missing_pat_trio_genotype.size();i++)
    {
      memset(matrix1,0,col*sizeof(double));
      const vector<inheritance> iv = missing_pat_trio_genotype[i]->get_inh_vec();
      if(iv.size()>0)
      {
        for(int j=0;j<iv.size();j++)
        {
          pair<int, int> hp_offspring = iv[j].hp_offspring;
          pair<int, int> hp_father = iv[j].hp_par1;

          matrix1[hap_pair_map[hp_offspring]] += geno_freq_table[hp_father]*transProbonep(hp_offspring,hp_father,hap_freq_table,geno_freq_table,1)/geno_freq_table[hp_offspring];
        }
        
        const vector<double> hap_num = missing_pat_trio_genotype[i]->get_hap_num();
        for(int j=0;j<dim;j++){
          for(int k=0;k<col;k++)
            result[j][k]+=hap_num[j]*matrix1[k];        
        }
      }  // trio hap config is possible
    }  // complete i
  }
  else if(typed_par1==0 && typed_par2==0)
  {
    for(int i=0;i<missing_pat_trio_genotype.size();i++)
    {
      memset(matrix1,0,col*sizeof(double));
      const vector<inheritance> iv = missing_pat_trio_genotype[i]->get_inh_vec();
      if(iv.size()>0)
      {
        for(int j=0;j<iv.size();j++)
        {
          pair<int, int> hp_offspring = iv[j].hp_offspring;

          matrix1[hap_pair_map[hp_offspring]] += 1;
        }
        
        const vector<double> hap_num = missing_pat_trio_genotype[i]->get_hap_num();
        for(int j=0;j<dim;j++){
          for(int k=0;k<col;k++)
            result[j][k]+=hap_num[j]*matrix1[k];        
        }
      }  // trio hap config is possible
    }  // complete i
  }

  delete[] matrix1;
}


void Missingtrio::cond_parents_one(double ***result, map<int, double> & hap_freq_table, map<pair<int, int>, double> & geno_freq_table, map<int, int> & hap_map){
// result is of dim*h_size*h_size
  int row = hap_freq_table.size();
  
  double **matrix1 = new double*[row];
  for(int i=0;i<row;i++){
    matrix1[i] = new double[row];
    memset(matrix1[i],0,row*sizeof(double));
  }
  
  if(typed_par1==1 && typed_par2==1)
  {
    for(int i=0;i<missing_pat_trio_genotype.size();i++)
    {
      for(int j=0;j<row;j++)
        memset(matrix1[j],0,row*sizeof(double));

      const vector<inheritance> iv = missing_pat_trio_genotype[i]->get_inh_vec();
      if(iv.size()>0)
      {
        for(int j=0;j<iv.size();j++)
        {
          pair<int, int> hp_offspring = iv[j].hp_offspring;
          pair<int, int> hp_father = iv[j].hp_par1;
          pair<int, int> hp_mother = iv[j].hp_par2;
          
          if(hp_father.first==hp_father.second && hp_mother.first==hp_mother.second){
            matrix1[hap_map[hp_father.first]][hap_map[hp_mother.first]] += hap_freq_table[hp_father.first]*hap_freq_table[hp_mother.first]*transProb(hp_offspring,hp_father,hp_mother,hap_freq_table,geno_freq_table,1);                                     
          }
          else if(hp_father.first==hp_father.second && hp_mother.first!=hp_mother.second){
            matrix1[hap_map[hp_father.first]][hap_map[hp_mother.first]] += hap_freq_table[hp_father.first]*hap_freq_table[hp_mother.second]*transProb(hp_offspring,hp_father,hp_mother,hap_freq_table,geno_freq_table,1);                                     
            matrix1[hap_map[hp_father.first]][hap_map[hp_mother.second]] += hap_freq_table[hp_father.first]*hap_freq_table[hp_mother.first]*transProb(hp_offspring,hp_father,hp_mother,hap_freq_table,geno_freq_table,1);                                               
          }
          else if(hp_father.first!=hp_father.second && hp_mother.first==hp_mother.second){
            matrix1[hap_map[hp_father.first]][hap_map[hp_mother.first]] += hap_freq_table[hp_father.second]*hap_freq_table[hp_mother.first]*transProb(hp_offspring,hp_father,hp_mother,hap_freq_table,geno_freq_table,1);                                     
            matrix1[hap_map[hp_father.second]][hap_map[hp_mother.first]] += hap_freq_table[hp_father.first]*hap_freq_table[hp_mother.first]*transProb(hp_offspring,hp_father,hp_mother,hap_freq_table,geno_freq_table,1);                                               
          }
          else if(hp_father.first!=hp_father.second && hp_mother.first!=hp_mother.second){
            matrix1[hap_map[hp_father.first]][hap_map[hp_mother.first]] += hap_freq_table[hp_father.second]*hap_freq_table[hp_mother.second]*transProb(hp_offspring,hp_father,hp_mother,hap_freq_table,geno_freq_table,1);                                     
            matrix1[hap_map[hp_father.first]][hap_map[hp_mother.second]] += hap_freq_table[hp_father.second]*hap_freq_table[hp_mother.first]*transProb(hp_offspring,hp_father,hp_mother,hap_freq_table,geno_freq_table,1);                                               
            matrix1[hap_map[hp_father.second]][hap_map[hp_mother.first]] += hap_freq_table[hp_father.first]*hap_freq_table[hp_mother.second]*transProb(hp_offspring,hp_father,hp_mother,hap_freq_table,geno_freq_table,1);                                     
            matrix1[hap_map[hp_father.second]][hap_map[hp_mother.second]] += hap_freq_table[hp_father.first]*hap_freq_table[hp_mother.first]*transProb(hp_offspring,hp_father,hp_mother,hap_freq_table,geno_freq_table,1);                                              
          }
        }
        
        const vector<double> hap_num = missing_pat_trio_genotype[i]->get_hap_num();
        for(int j=0;j<dim;j++){
          for(int k=0;k<row;k++)
            for(int l=0;l<row;l++)
              result[j][k][l]+=hap_num[j]*matrix1[k][l];
        }        
      }  // trio hap config is possible
    }  // complete i
  }
  else if(typed_par1==1 && typed_par2==0)
  {
    for(int i=0;i<missing_pat_trio_genotype.size();i++)
    {
      for(int j=0;j<row;j++)
        memset(matrix1[j],0,row*sizeof(double));

      const vector<inheritance> iv = missing_pat_trio_genotype[i]->get_inh_vec_complete();
      if(iv.size()>0)
      {
        for(int j=0;j<iv.size();j++)
        {
          pair<int, int> hp_offspring = iv[j].hp_offspring;
          pair<int, int> hp_father = iv[j].hp_par1;
          pair<int, int> hp_mother = iv[j].hp_par2;

          if(hp_father.first==hp_father.second && hp_mother.first==hp_mother.second){
            matrix1[hap_map[hp_father.first]][hap_map[hp_mother.first]] += hap_freq_table[hp_father.first]*hap_freq_table[hp_mother.first]*transProb(hp_offspring,hp_father,hp_mother,hap_freq_table,geno_freq_table,1);                                     
          }
          else if(hp_father.first==hp_father.second && hp_mother.first!=hp_mother.second){
            matrix1[hap_map[hp_father.first]][hap_map[hp_mother.first]] += hap_freq_table[hp_father.first]*hap_freq_table[hp_mother.second]*transProb(hp_offspring,hp_father,hp_mother,hap_freq_table,geno_freq_table,1);                                     
            matrix1[hap_map[hp_father.first]][hap_map[hp_mother.second]] += hap_freq_table[hp_father.first]*hap_freq_table[hp_mother.first]*transProb(hp_offspring,hp_father,hp_mother,hap_freq_table,geno_freq_table,1);                                               
          }
          else if(hp_father.first!=hp_father.second && hp_mother.first==hp_mother.second){
            matrix1[hap_map[hp_father.first]][hap_map[hp_mother.first]] += hap_freq_table[hp_father.second]*hap_freq_table[hp_mother.first]*transProb(hp_offspring,hp_father,hp_mother,hap_freq_table,geno_freq_table,1);                                     
            matrix1[hap_map[hp_father.second]][hap_map[hp_mother.first]] += hap_freq_table[hp_father.first]*hap_freq_table[hp_mother.first]*transProb(hp_offspring,hp_father,hp_mother,hap_freq_table,geno_freq_table,1);                                               
          }
          else if(hp_father.first!=hp_father.second && hp_mother.first!=hp_mother.second){
            matrix1[hap_map[hp_father.first]][hap_map[hp_mother.first]] += hap_freq_table[hp_father.second]*hap_freq_table[hp_mother.second]*transProb(hp_offspring,hp_father,hp_mother,hap_freq_table,geno_freq_table,1);                                     
            matrix1[hap_map[hp_father.first]][hap_map[hp_mother.second]] += hap_freq_table[hp_father.second]*hap_freq_table[hp_mother.first]*transProb(hp_offspring,hp_father,hp_mother,hap_freq_table,geno_freq_table,1);                                               
            matrix1[hap_map[hp_father.second]][hap_map[hp_mother.first]] += hap_freq_table[hp_father.first]*hap_freq_table[hp_mother.second]*transProb(hp_offspring,hp_father,hp_mother,hap_freq_table,geno_freq_table,1);                                     
            matrix1[hap_map[hp_father.second]][hap_map[hp_mother.second]] += hap_freq_table[hp_father.first]*hap_freq_table[hp_mother.first]*transProb(hp_offspring,hp_father,hp_mother,hap_freq_table,geno_freq_table,1);                                              
          }
        }
        
        const vector<double> hap_num = missing_pat_trio_genotype[i]->get_hap_num();
        for(int j=0;j<dim;j++){
          for(int k=0;k<row;k++)
            for(int l=0;l<row;l++)
              result[j][k][l]+=hap_num[j]*matrix1[k][l];
        }
      }  // trio hap config is possible
    }  // complete i
  }

  for(int i=0;i<row;i++)
    delete[] matrix1[i];
  delete[] matrix1;
}


void Missingtrio::cond_parents_two(double ***result, map<int, double> & hap_freq_table, map<pair<int, int>, double> & geno_freq_table, map<pair<int, int>, int> & hap_pair_map){  // for sibs
// result is of dim*g_size*g_size
  int col = geno_freq_table.size();
  
  double **matrix1 = new double*[col];
  for(int i=0;i<col;i++){
    matrix1[i] = new double[col];
    memset(matrix1[i],0,col*sizeof(double));
  }
  
  if(typed_par1==1 && typed_par2==1)
  {
    for(int i=0;i<missing_pat_trio_genotype.size();i++)
    {
      for(int j=0;j<col;j++)
        memset(matrix1[j],0,col*sizeof(double));

      const vector<inheritance> iv = missing_pat_trio_genotype[i]->get_inh_vec();
      if(iv.size()>0)
      {
        for(int j=0;j<iv.size();j++)
        {
          pair<int, int> hp_offspring = iv[j].hp_offspring;
          pair<int, int> hp_father = iv[j].hp_par1;
          pair<int, int> hp_mother = iv[j].hp_par2;
          
          matrix1[hap_pair_map[hp_father]][hap_pair_map[hp_mother]] += transProb(hp_offspring,hp_father,hp_mother,hap_freq_table,geno_freq_table,1);
        }
        
        const vector<double> hap_num = missing_pat_trio_genotype[i]->get_hap_num();
        for(int j=0;j<dim;j++){
          for(int k=0;k<col;k++)
            for(int l=0;l<col;l++)
              result[j][k][l]+=hap_num[j]*matrix1[k][l];
        }        
      }  // trio hap config is possible
    }  // complete i
  }
  else if(typed_par1==1 && typed_par2==0)
  {
    for(int i=0;i<missing_pat_trio_genotype.size();i++)
    {
      for(int j=0;j<col;j++)
        memset(matrix1[j],0,col*sizeof(double));

      const vector<inheritance> iv = missing_pat_trio_genotype[i]->get_inh_vec_complete();
      if(iv.size()>0)
      {
        for(int j=0;j<iv.size();j++)
        {
          pair<int, int> hp_offspring = iv[j].hp_offspring;
          pair<int, int> hp_father = iv[j].hp_par1;
          pair<int, int> hp_mother = iv[j].hp_par2;

          matrix1[hap_pair_map[hp_father]][hap_pair_map[hp_mother]] += transProb(hp_offspring,hp_father,hp_mother,hap_freq_table,geno_freq_table,1);
        }
        
        const vector<double> hap_num = missing_pat_trio_genotype[i]->get_hap_num();
        for(int j=0;j<dim;j++){
          for(int k=0;k<col;k++)
            for(int l=0;l<col;l++)
              result[j][k][l]+=hap_num[j]*matrix1[k][l];
        }
      }  // trio hap config is possible
    }  // complete i
  }

  for(int i=0;i<col;i++)
    delete[] matrix1[i];
  delete[] matrix1;
}


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

  for(int i=0;i<s.size();i++){
    if(s[i]==m)
      return i;
  }
  return -1;
}