#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;

  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];
      memset(deriv1,0,dim*sizeof(double));

      deriv2 = new double[dim];
      memset(deriv2,0,dim*sizeof(double));

      deriv3 = new double[dim];
      memset(deriv3,0,dim*sizeof(double));
    }
    else if(typed_par1==0 && typed_par2==0)
    {
      deriv_r = new double[dim];
      memset(deriv_r,0,dim*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;

  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));
    }

    deriv_r = new double[dim];
    memset(deriv_r,0,dim*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, map<int, double> & hap_weight, int **count)
{
  memset(deriv1,0,dim*sizeof(double));
  memset(deriv2,0,dim*sizeof(double));
  memset(deriv3,0,dim*sizeof(double));

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

  double matrix1 = 0;
  double *matrix2 = new double[dim];
  memset(matrix2,0,dim*sizeof(double));
  double matrix3 = 0;
  double *matrix4 = new double[dim];
  memset(matrix4,0,dim*sizeof(double));
  double matrix5 = 0;
  double *matrix6 = new double[dim];
  memset(matrix6,0,dim*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)
        {
          matrix1 = 0;
          memset(matrix2,0,dim*sizeof(double));
          matrix3 = 0;
          memset(matrix4,0,dim*sizeof(double));
          matrix5 = 0;
          memset(matrix6,0,dim*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 = geno_freq_table[hp_father]*geno_freq_table[hp_mother]*transProb(hp_offspring,hp_father,hp_mother,hap_freq_table,geno_freq_table,1);
  
            value1 = value*(hap_weight[hp_offspring.first]+hap_weight[hp_offspring.second]);
            value2 = value*(hap_weight[hp_father.first]+hap_weight[hp_father.second]);
            value3 = value*(hap_weight[hp_mother.first]+hap_weight[hp_mother.second]);
            matrix1+=value1;
            matrix3+=value2;
            matrix5+=value3;

            for(int k=0;k<dim;k++)
            {
              if(count[sites[k]][loc]>0)
              {
                matrix2[k]+=count[sites[k]][loc]*value1;
                matrix4[k]+=count[sites[k]][loc]*value2;
                matrix6[k]+=count[sites[k]][loc]*value3;
              }
            }
          }

          for(int j=0;j<dim;j++)
          {
            deriv1[j]-=(matrix2[j]/2.0-weight[j][i]*matrix1);
            deriv2[j]-=(matrix4[j]/2.0-weight[j][i]*matrix3);
            deriv3[j]-=(matrix6[j]/2.0-weight[j][i]*matrix5);
          }
        }
      }
    }
  }
  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)
        {
          matrix1 = 0;
          memset(matrix2,0,dim*sizeof(double));
          matrix3 = 0;
          memset(matrix4,0,dim*sizeof(double));
          matrix5 = 0;
          memset(matrix6,0,dim*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 = geno_freq_table[hp_father]*transProbonep(hp_offspring,hp_father,hap_freq_table,geno_freq_table,1);
  
            value1 = value*(hap_weight[hp_offspring.first]+hap_weight[hp_offspring.second]);
            value2 = value*(hap_weight[hp_father.first]+hap_weight[hp_father.second]);
            int trans_hap;
            if(hp_father.first==hp_father.second)
            {
              if(hp_offspring.first==hp_father.first)
                trans_hap = hp_offspring.second;
              else if(hp_offspring.second==hp_father.first)
                trans_hap = hp_offspring.first;
            }
            else
            {
              if(hp_offspring.first==hp_offspring.second)
              {
                if(hp_offspring.first==hp_father.first || hp_offspring.first==hp_father.second)
                  trans_hap = hp_offspring.second;
              }
              else
              {
                if(hp_offspring.first==hp_father.first && hp_offspring.second==hp_father.second)
                  trans_hap = -1;
                else if(hp_offspring.first==hp_father.first || hp_offspring.first==hp_father.second)
                  trans_hap = hp_offspring.second;
                else if(hp_offspring.second==hp_father.first || hp_offspring.second==hp_father.second)
                  trans_hap = hp_offspring.first;
              }
            }
            if(trans_hap==-1)
              value3 = value*(hap_freq_table[hp_offspring.first]*hap_weight[hp_offspring.first]+hap_freq_table[hp_offspring.second]*hap_weight[hp_offspring.second])/(hap_freq_table[hp_offspring.first]+hap_freq_table[hp_offspring.second]);
            else
              value3 = value*hap_weight[trans_hap];
            matrix1+=value1;
            matrix3+=value2;
            matrix5+=value3;

            for(int k=0;k<dim;k++)
            {
              if(count[sites[k]][loc]>0)
              {
                matrix2[k]+=count[sites[k]][loc]*value1;
                matrix4[k]+=count[sites[k]][loc]*value2;
                matrix6[k]+=count[sites[k]][loc]*value3;
              }
            }
          }

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

  delete[] matrix2;
  delete[] matrix4;
  delete[] matrix6;
}


void Missingtrio::compute_deriv_r(double *freq, map<pair<int, int>, double> & geno_freq_table, map<pair<int, int>, int> & hap_pair_map, map<int, double> & hap_weight, vector<double> & h_weight, int **count)
{
  memset(deriv_r,0,dim*sizeof(double));

  int loc;
  double value;

  double matrix1 = 0;
  double *matrix2 = new double[dim];
  memset(matrix2,0,dim*sizeof(double));

  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()>1)
    {
      matrix1 = 0;
      memset(matrix2,0,dim*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];
        value = geno_freq_table[hp_offspring]*(hap_weight[hp_offspring.first]+hap_weight[hp_offspring.second]);
        matrix1+=value;

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

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

  for(int i=0;i<dim;i++)
    deriv_r[i]+=freq[sites[i]]*h_weight[sites[i]];

  delete[] matrix2;
}


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]);
  }
}


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)
    {
      delete[] deriv1;
      delete[] deriv2;
      delete[] deriv3;
    }
    else if(typed_par1==0 && typed_par2==0)
      delete[] deriv_r;
  }
}


// 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;
}