#include "includes.h"
#include "grid.h"
#include "vector.h"

#include "time.h"

#include <fstream.h>
extern ofstream errorfs;

#include "declare.h"


/* This file contains source code for estimating two-locus haplotype
   frequencies jointly using an HMM-- this complication is necessary
   to get "consistent" estimates across overlapping pairs of loci
   when there are missing alleles, ie., alleles not typed */

/* Following this, there is source code for adjustments to these
   frequencies for the following cases:
   1. Allele appears in affected genotypes/haps. but not in controls 
   (assigning freq. of 0 seems to harsh and will lead to loglik of -inf)
   2. Haplotype appears in affecteds (or is compatible with genotypes)
   but either does not appear in controls or the estimate of freq.
   is vanishingly small (once again, seems too harsh a penalty)

   This is not trivial, because changing one joint frequency could
   change freq's at neighboring marker-pairs, etc., if not done correctly.
*/

void makehap3fr_hmm_HconBAYES(dvectors &h3_freq_lists, 
			      const ivectors &a_lists, 
			      const igrid &controls, const int &bayes,
			      const int &debug)
{

  //estimates three-marker haplotype frequencies, using HMM, assuming
  //haplotype controls

  int l,h,i,j,k,iter,conv;
  double diff,sum;
  const int numloc=controls.dimy();
  const int numhap=controls.dimx();
  //double llik=0.0,lik=0.0;

  //make vector indicating haplotypes with missing values
  int is_missing=0;
  ivector missing; missing.init(numhap);
  make_missing_Hfr(missing, is_missing, controls, 0);

  //delete me
  for (i=0;i<missing.dim();i++) missing(i)=0;

  //declare and intialize and clear lists of two-locus hap freqs 
  //corresponding to a_lists
  dvectors h2_freq_lists; 
  h2_freq_lists.init(numloc-1);
  for (l=0;l<numloc-1;l++) 
    { h2_freq_lists[l].init(a_lists[l].dim()*a_lists[l+1].dim()); 
    h2_freq_lists[l].clear(0); }
  
  //declare and initialize and clear hmm objects
  //states are two locus haplotypes
  dvectorss b(numhap);
  dvectorss alpha(numhap),beta(numhap),gamma(numhap),xi(numhap);
  for (h=0;h<numhap;h++) {
    if (missing[h]==0) {
      b[h].init(numloc-1);
      alpha[h].init(numloc-1);
      beta[h].init(numloc-1);
      gamma[h].init(numloc-1);
      for (l=0;l<numloc-1;l++) {
	b[h][l].init(a_lists[l].dim()*a_lists[l+1].dim()); b[h][l].clear(0);
	alpha[h][l].init(b[h][l].dim()); alpha[h][l].clear(0);
	beta[h][l].init(b[h][l].dim()); beta[h][l].clear(0);
      }
    }
  }
  //initialize xi vectors for all haps
  for (h=0;h<numhap;h++) {
    xi[h].init(numloc-2);
    for (l=0;l<numloc-2;l++) {
      xi[h][l].init(a_lists[l].dim()*a_lists[l+1].dim()*
		    a_lists[l+1].dim()*a_lists[l+2].dim());
      xi[h][l].clear(0); }
  }
  
  dvector init;
  init.init( a_lists[0].dim()*a_lists[1].dim() ); 
  init.clear(0.0);

  //make copy to determine convergence
  dvector initC; initC.init( init.dim() ); initC.clear(0.0);
  

  dvectors trans(numloc-2);
  for (l=0;l<numloc-2;l++) {
    trans[l].init(
		  a_lists[l].dim()*a_lists[l+1].dim()*
		  a_lists[l+2].dim());
    trans[l].clear(0); }
  
  //make copy of trans to determine convergence
  dvectors transC(numloc-2);
  for (l=0;l<numloc-2;l++) {
    transC[l].init(trans[l].dim()); transC[l].clear(0); }

  //assign initial values to parameters
  for (i=0;i<init.dim();i++) 
    init(i)=1.0/ ( init.dim() );
  for (l=0;l<numloc-2;l++) {
    for (i=0;i<a_lists[l].dim();i++) {
      for (j=0;j<a_lists[l+1].dim();j++) {
	for (k=0;k<a_lists[l+2].dim();k++) {
	  trans[l][ 
		   i*a_lists[l+1].dim()*a_lists[l+2].dim()+
		   j*a_lists[l+2].dim()+k ]=
	    1.0/( a_lists[l+2].dim() ); } } } }

  //make h3_freq_lists from init,trans
  //makeh3freqlistsmkv2(h3_freq_lists,h2_freq_lists,init,trans,a_lists,0);
  //makeh2freqlistsmkv2(h2_freq_lists,h3_freq_lists,a_lists,0);

  //fill in xi for haplotypes with no missing alleles
  makeximkv2_nomissing_Hfr(xi, missing, a_lists, controls, 0);

  //make observation probabilities [b]
  makebmkv2_Hfr(b, missing, a_lists, controls, 0);

  //iterate hmm
  int maxiter=200;
  for (iter=1;iter<=maxiter;iter++) {

    //copy init to use later in checking convergence
    initC=init;

    //copy hapl. freq. to use later in checking convergence
    for (l=0;l<transC.dim();l++) 
      transC[l]=trans[l];

    //make trans (from h3_freq_lists)
    //maketransmkv2_Hfr(trans, h3_freq_lists, h2_freq_lists, a_lists,0);

    //make alpha
    makealphamkv2_HfrBAYES(alpha, missing, a_lists, init, trans, b, 0);

    //make beta
    makebetamkv2_HfrBAYES(beta, missing, a_lists, trans, b, 0);

    //make xi
    makeximkv2_HfrBAYES(xi, missing, a_lists, alpha, beta, trans, b, 0);

    //compute likelihood from HMM
    //cout << "iter " << iter << " ";
    //llik=0; 
    //for (h=0;h<numhap;h++) {
    //lik=0.0;
    //for (i=0;i<alpha[h][0].dim();i++) lik+=alpha[h][0][i]*beta[h][0][i];
    //llik+=log(lik);
    //}
    //cout << "llik " << llik << endl;

    //update h3_freq_lists
    //makeh3freqlistsHfrxi(h3_freq_lists, xi, a_lists, 0);
    //makeh2freqlistsmkv2(h2_freq_lists,h3_freq_lists,a_lists,0);
    
    //give expected counts of 3-locus haplotypes
    makeh3countsHfrxi(h3_freq_lists, xi, a_lists, 0);
    //from this make expected counts of two-locus haplotypes
    makeh2countsGmkv2(h2_freq_lists,numhap/2,h3_freq_lists,a_lists,0);

    //Bayesian adjustment/update to init
    init=h2_freq_lists[0];
    for (i=0;i<init.dim();i++) { 
      init(i)+=1.0*bayes; 
      init(i)/=numhap+1.0*init.dim()*bayes; }

    //check that init sums to 1
    sum=0;
    for (i=0;i<init.dim();i++) sum+=init(i);
    if ( !(0.9999<sum && sum<1.0001) ) {
      errorfs << "ERROR: Bayesian adjustment to init in hap. freq. est.; "
	      << "sums to " << sum << endl;
      exit(1); }

    //    //Bayesian adjustment/update to trans
    //for (l=0;l<numloc-2;l++) {
    //for (i=0;i<a_lists[l  ].dim();i++) {
    //for (j=0;j<a_lists[l+1].dim();j++) {
    //  for (k=0;k<a_lists[l+2].dim();k++) {
    //    trans[l][i*a_lists[l+1].dim()*a_lists[l+2].dim()+
    //	    j*a_lists[l+2].dim()+k]=
    //      ( h3_freq_lists[l][i*a_lists[l+1].dim()*a_lists[l+2].dim()+
    //			j*a_lists[l+2].dim()+k]    +1 )/
    //      ( h2_freq_lists[l][i*a_lists[l+1].dim()+j] + a_lists[l+2].dim());
    //  }
    //}
    //}
    //}

    //Bayesian adjustment/update to trans
    for (l=0;l<numloc-2;l++) {
      for (i=0;i<a_lists[l  ].dim();i++) {
	for (j=0;j<a_lists[l+1].dim();j++) {
	  for (k=0;k<a_lists[l+2].dim();k++) {
	    if (h3_freq_lists[l][i*a_lists[l+1].dim()*a_lists[l+2].dim()+
				j*a_lists[l+2].dim()+k]  +1.0*bayes <1e-100
		&&
		h2_freq_lists[l][i*a_lists[l+1].dim()+j]+ 
		1.0*a_lists[l+2].dim()*bayes<1e-100)
	      trans[l][i*a_lists[l+1].dim()*a_lists[l+2].dim()+
		      j*a_lists[l+2].dim()+k]=0.0;
	    else
	      trans[l][i*a_lists[l+1].dim()*a_lists[l+2].dim()+
		      j*a_lists[l+2].dim()+k]=
		( h3_freq_lists[l][i*a_lists[l+1].dim()*a_lists[l+2].dim()+
				  j*a_lists[l+2].dim()+k]  +1.0*bayes )/
		( h2_freq_lists[l][i*a_lists[l+1].dim()+j]+ 
		  1.0*a_lists[l+2].dim()*bayes);
	  }
	  //if they're all 0, set first to 1
	  conv=0;
	  for (k=0;k<a_lists[l+2].dim();k++) 
	    if (trans[l][i*a_lists[l+1].dim()*a_lists[l+2].dim()+
			j*a_lists[l+2].dim()+k]>0) { conv=1; break; }
	  if (conv==0) trans[l][i*a_lists[l+1].dim()*a_lists[l+2].dim()+
			j*a_lists[l+2].dim()  ]=1;
	}
      }
    }//l

    //check trans
    for (l=0;l<numloc-2;l++) {
      for (i=0;i<a_lists[l  ].dim();i++) {
	for (j=0;j<a_lists[l+1].dim();j++) {
	  sum=0;
	  for (k=0;k<a_lists[l+2].dim();k++) {
	    sum+=trans[l][
			  i*a_lists[l+1].dim()*a_lists[l+2].dim()+
			  j*a_lists[l+2].dim()+k]; }
	  if ( !(0.9999<sum && sum<1.0001) ) {
	    errorfs << "ERROR in Bayes adj. of trans. Doesn't sum to 1 for "
		    << "i " << i << " j " << j << " k " << k 
		    << " at locus " << l << " ( " << sum << " ) " << endl;
	    exit(1); }}}}


    //check convergence
    conv=0;
    for (i=0;i<init.dim();i++) {
      diff=init(i)-initC(i);
      if ( diff<-1e-8 || 1e-8<diff ) { conv=1; break; } }
    if (conv==0) {
      for (l=0;l<numloc-2;l++) {
	for (i=0;i<trans[l].dim();i++) {
	  diff=trans[l][i]-transC[l][i];
	  if ( diff<-1e-8 || 1e-8<diff ) { conv=1; break; } }
	if (conv==1) break; } }
    if (conv==0) break; 

  }//iter

  //make h3_freq_lists from trans, init
  makeh3freqlistsGmkv2(h3_freq_lists,h2_freq_lists,init,trans,a_lists,0);

  if (debug) {
    cout << "iter " << iter << endl;
    cout << "init " << endl;
    for (i=0;i<init.dim();i++) cout << init(i) << " "; cout << endl;
    cout << "trans " << endl;
    for (l=0;l<numloc-2;l++) {
      for (i=0;i<trans[l].dim();i++) {
	cout << trans[l][i] << " "; }
      cout << endl; }

    cout << "h3_freq_lists" << endl;
    for (l=0;l<numloc-2;l++) {
      for (i=0;i<h3_freq_lists[l].dim();i++) {	
	cout << h3_freq_lists[l][i] << " "; }
      cout << endl;
    }

    //make likelihood
    //make alpha
    makealphamkv2_HfrBAYES(alpha, missing, a_lists, init, trans, b, 0);
    double lik=0,loglik=0;
    for (h=0;h<numhap;h++) {
      lik=0;
      for (i=0;i<alpha[h][numloc-2].dim();i++) {
	lik+=alpha[h][numloc-2][i]; }
      if ( !(0<=lik && lik<=1.0) ) {
	errorfs << "ERROR: lik. for hap " << h << "+1 is " << lik << endl;
	exit(1); }
      loglik+=log(lik); }
    cout << "lik " << loglik << endl; 
  }
}

void makehap3fr_hmm_Hcon(dvectors &h3_freq_lists, 
			 const ivectors &a_lists, 
			 const igrid &controls, 
			 const int &debug)
{

  //estimates three-marker haplotype frequencies, using HMM, assuming
  //haplotype controls

  int l,h,i,j,k,iter,conv;
  double diff;
  const int numloc=controls.dimy();
  const int numhap=controls.dimx();
  //double llik=0.0,lik=0.0;

  //make vector indicating haplotypes with missing values
  int is_missing=0;
  ivector missing; missing.init(numhap);
  make_missing_Hfr(missing, is_missing, controls, 0);

  //delete me
  //for (i=0;i<missing.dim();i++) missing(i)=0;

  //declare and intialize and clear lists of two-locus hap freqs 
  //corresponding to a_lists
  dvectors h2_freq_lists; 
  h2_freq_lists.init(numloc-1);
  for (l=0;l<numloc-1;l++) 
    { h2_freq_lists[l].init(a_lists[l].dim()*a_lists[l+1].dim()); 
    h2_freq_lists[l].clear(0); }
  
  //declare and initialize and clear copy to be used to determine convergence
  dvectors h3_freq_lists_copy(h3_freq_lists.dim());
  for (l=0;l<numloc-2;l++) 
    { 
      h3_freq_lists_copy[l].init(h3_freq_lists[l].dim()); 
      h3_freq_lists_copy[l].clear(0); 
    }
  
  //declare and initialize and clear hmm objects
  //states are two locus haplotypes
  dvectorss b(numhap);
  dvectorss alpha(numhap),beta(numhap),gamma(numhap),xi(numhap);
  for (h=0;h<numhap;h++) {
    if (missing[h]==0) {
      b[h].init(numloc-1);
      alpha[h].init(numloc-1);
      beta[h].init(numloc-1);
      gamma[h].init(numloc-1);
      for (l=0;l<numloc-1;l++) {
	b[h][l].init(a_lists[l].dim()*a_lists[l+1].dim()); b[h][l].clear(0);
	alpha[h][l].init(b[h][l].dim()); alpha[h][l].clear(0);
	beta[h][l].init(b[h][l].dim()); beta[h][l].clear(0);
      }
    }
  }
  //initialize xi vectors for all haps
  for (h=0;h<numhap;h++) {
    xi[h].init(numloc-2);
    for (l=0;l<numloc-2;l++) {
      xi[h][l].init(a_lists[l].dim()*a_lists[l+1].dim()*
		    a_lists[l+1].dim()*a_lists[l+2].dim());
      xi[h][l].clear(0); }
  }
  
  dvector init;
  init.init( a_lists[0].dim()*a_lists[1].dim() ); 
  init.clear(0.0);

  dvectors trans(numloc-2);
  for (l=0;l<numloc-2;l++) {
    trans[l].init(
		  a_lists[l].dim()*a_lists[l+1].dim()*
		  a_lists[l+1].dim()*a_lists[l+2].dim());
    trans[l].clear(0); }
  
  //assign initial values to parameters
  for (i=0;i<init.dim();i++) 
    init(i)=1.0/ ( a_lists[0].dim()*a_lists[1].dim() );
  for (l=0;l<numloc-2;l++) {
    for (i=0;i<a_lists[l].dim();i++) {
      for (j=0;j<a_lists[l+1].dim();j++) {
	for (k=0;k<a_lists[l+2].dim();k++) {
	  trans[l][ 
		   (i*a_lists[l+1].dim()+j)* 
		   a_lists[l+1].dim()*a_lists[l+2].dim()
		   +
		   (j*a_lists[l+2].dim()+k) ]=
	    1.0/( a_lists[l+2].dim() ); } } } }
  
  //make h3_freq_lists from init,trans
  makeh3freqlistsmkv2(h3_freq_lists,h2_freq_lists,init,trans,a_lists,0);
  makeh2freqlistsmkv2(h2_freq_lists,h3_freq_lists,a_lists,0);

  //fill in xi for haplotypes with no missing alleles
  makeximkv2_nomissing_Hfr(xi, missing, a_lists, controls, 0);

  //make observation probabilities [b]
  makebmkv2_Hfr(b, missing, a_lists, controls, 0);

  //iterate hmm
  int maxiter=200;
  for (iter=1;iter<=maxiter;iter++) {

    //copy hapl. freq. to use later in checking convergence
    for (l=0;l<h3_freq_lists.dim();l++) 
      h3_freq_lists_copy[l]=h3_freq_lists[l];
    
    //make trans (from h3_freq_lists)
    maketransmkv2_Hfr(trans, h3_freq_lists, h2_freq_lists, a_lists,0);

    //make alpha
    makealphamkv2_Hfr(alpha, missing, a_lists, init, trans, b, 0);

    //make beta
    makebetamkv2_Hfr(beta, missing, trans, b, 0);

    //make xi
    makeximkv2_Hfr(xi, missing, alpha, beta, trans, b, 0);

    //compute likelihood from HMM
    //cout << "iter " << iter << " ";
    //llik=0; 
    //for (h=0;h<numhap;h++) {
    //lik=0.0;
    //for (i=0;i<alpha[h][0].dim();i++) lik+=alpha[h][0][i]*beta[h][0][i];
    //llik+=log(lik);
    //}
    //cout << "llik " << llik << endl;

    //update h3_freq_lists
    makeh3freqlistsHfrxi(h3_freq_lists, xi, a_lists, 0);
    makeh2freqlistsmkv2(h2_freq_lists,h3_freq_lists,a_lists,0);
    
    //check convergence
    conv=0;
    for (l=0;l<numloc-2;l++) {
      for (i=0;i<h3_freq_lists[l].dim();i++) {
	diff=h3_freq_lists[l][i]-h3_freq_lists_copy[l][i];
	//cout << "diff " << is_missing << " " 
	//   << l << " " << i << " " << diff << endl;
	if ( diff<-1e-8 || 1e-8<diff ) { conv=1; break; } }
      if (conv==1) break; }// || is_missing==1) break; }
    if (conv==0) break;

  }//iter

  if (debug) {
    cout << "h3_freq_lists" << endl;
    for (l=0;l<numloc-2;l++) {
      for (i=0;i<h3_freq_lists[l].dim();i++) {	
	cout << h3_freq_lists[l][i] << " "; }
      cout << endl;
    }
  }
}

void makeximkv2_nomissing_Hfr(dvectorss &xi, const ivector &missing, 
			      const ivectors &a_lists, 
			      const igrid &controls, const int &debug)
{
  //makes xi, for control haplotypes without a missing allele

  int i,l,h;
  const int numhap=controls.dimx();
  const int numloc=controls.dimy();
  int alp0=0,alp1=0,alp2=0;

  for (h=0;h<numhap;h++) {

    if (missing[h]==1) {
      for (l=0;l<numloc-2;l++) {

	alp0= getposinlist( controls(h,l  ), a_lists[l] );
	alp1= getposinlist( controls(h,l+1), a_lists[l+1] );
	alp2= getposinlist( controls(h,l+2), a_lists[l+2] );

	xi[h][l][ 
		 ( alp0*a_lists[l+1].dim()+alp1 )
		 *a_lists[l+1].dim()*a_lists[l+2].dim()+
		 ( alp1*a_lists[l+2].dim()+alp2 )
	]=1.0;
      }
    }
  }

  if (debug) {
    cout << "xi [in makeximkv2_nomissing_Hfr] " << endl;
    for (h=0;h<numhap;h++) {
      if (missing[h]==1) {
	cout << "hap " << h << endl;
	for (l=0;l<numloc-2;l++) {
	  for (i=0;i<xi[h][l].dim();i++) {	
	    cout << xi[h][l][i] << " "; }
	  cout << endl;
	}
      }
    }
  }

}

void makebmkv2_Hfr(dvectorss &b, const ivector &missing, 
		   const ivectors &a_lists,
		   const igrid &controls, const int &debug)
{
  //makes observation probabilities, assuming haplotype data

  int i,h,l;

  const int numhap=controls.dimx();
  const int numloc=controls.dimy();

  for (h=0;h<numhap;h++) {

    if (missing[h]==0) {
      
      for (l=0;l<numloc-1;l++) {
	
	if ( controls(h,l)!=0) {
	  if (controls(h,l+1)!=0) {
	    b[h][l][ getposinlist( controls(h,l), a_lists[l] )*
		   a_lists[l+1].dim()+
		   getposinlist( controls(h,l+1), a_lists[l+1] ) ]
	      =1.0;
	  }
	  else {//controls(h,l+1)==0
	    for (i=0;i<a_lists[l+1].dim();i++)
	      b[h][l][ getposinlist( controls(h,l), a_lists[l] )*
		     a_lists[l+1].dim()+ i ]
		=1.0;
	  }
	}
	else {//controls(h,l)==0
	  if (controls(h,l+1)!=0) {
	    for (i=0;i<a_lists[l].dim();i++)
	      b[h][l][ 
		       i*a_lists[l+1].dim()+ 
		       getposinlist( controls(h,l+1), a_lists[l+1] ) ]
		=1.0;	    
	  }
	  else {//controls(h,l+1)==0
	    for (i=0;i<b[h][l].dim();i++)
	      b[h][l][i]=1.0;
	  }
	}

      }
    }
  }

  if (debug) {
    cout << "b [obs. probs.] " << endl;
    for (h=0;h<numhap;h++) {
      for (l=0;l<numloc-1;l++) {
	for (i=0;i<b[h][l].dim();i++) {
	  cout << b[h][l][i] << " "; }
	cout << endl;
      }
    }
  }
  
}

void maketransmkv2_Hfr(dvectors &trans, const dvectors &h3_freq_lists,
		       const dvectors &h2_freq_lists,
		       const ivectors &a_lists, const int &debug)
{
  //makes transition probabilities for haplotype HMM

  int i,j,k,l;
  const int numloc=a_lists.dim();
  double value;

  //clear trans
  for (l=0;l<numloc-2;l++) {
    trans[l].clear(0.0); }
  
  for (l=0;l<numloc-2;l++) {

    for (i=0;i<a_lists[l  ].dim();i++) {
      for (j=0;j<a_lists[l+1].dim();j++) {
	for (k=0;k<a_lists[l+2].dim();k++) {
	
	if (h2_freq_lists[l][i*a_lists[l+1].dim()+j]!=0.0)
	  trans[l][ 
		   (i*a_lists[l+1].dim()+j)*
		   h2_freq_lists[l+1].dim()+
		   (j*a_lists[l+2].dim()+k) ]=	    
	    h3_freq_lists[l][
			     i*a_lists[l+1].dim()*a_lists[l+2].dim()+
			     j*a_lists[l+2].dim()+
			     k ]/
	    h2_freq_lists[l][ i*a_lists[l+1].dim()+j ];
	}
      }
    }
  }

  if (debug) {
    cout << "trans [in maketransmkv2_Hfr] " << endl;
    for (l=0;l<numloc-2;l++) {
      for (i=0;i<trans[l].dim();i++) {	
	cout << trans[l][i] << " "; }
      cout << endl;
    }
  }

  //check that P(omega)=1
  for (l=0;l<numloc-2;l++) {
    
    for (i=0;i<h2_freq_lists[l].dim();i++) {
      value=0.0;
      for (j=0;j<h2_freq_lists[l+1].dim();j++) {
	value+=trans[l][i*h2_freq_lists[l+1].dim()+j];
	//cout << l << " " << i << " " << j << " " 
	//   << trans[l][i*h2_freq_lists[l+1].dim()+j] << endl;
      }
      if ( !( (0.9999<value && value<1.0001) || value<0.01) ) {
	errorfs << "ERROR: in maketransmkv2_Hfr [" << l 
		<< "], probs do not sum to 1 ["
		<< value << "] " << endl;
	exit(1); } } }
  

}

void makealphamkv2_Hfr(dvectorss &alpha, const ivector &missing,
		       const ivectors &a_lists,
		       const dvector &init,
		       const dvectors &trans, const dvectorss &b,
		       const int &debug)
{
  //makes alpha object, assuming haplotype controls

  int h,i,j,l;
  const int numhap=alpha.dim();
  const int numloc=a_lists.dim();

  //clear alpha
  for (h=0;h<numhap;h++) {
    if (missing[h]==0) {
      for (l=0;l<numloc-1;l++) {
	alpha[h][l].clear(0); } } }

  for (h=0;h<numhap;h++) {

    if (missing[h]==0) {
      //make first column
      for (i=0;i<alpha[h][0].dim();i++) {

	if (b[h][0][i] != 0) 
	  alpha[h][0][i]=
	    init(i)*b[h][0][i];
	
      }
    
    
      //make rest of alpha
      for (l=0;l<numloc-2;l++) {
      
	for (j=0;j<alpha[h][l+1].dim();j++) {
	
	  if (b[h][l+1][j] != 0) { 
	    for (i=0;i<alpha[h][l].dim();i++) {
	    
	      if (alpha[h][l][i] != 0) {
		alpha[h][l+1][j]+=
		alpha[h][l][i]*
		  trans[l][i*alpha[h][l+1].dim()+j]*
		  b[h][l+1][j];

	      }
	    }
	  }
	}
	
      }
    }
  }

  if (debug) {
    cout << "alpha [in maketransmkv2_Hfr] " << endl;
    for (h=0;h<numhap;h++) {
      for (l=0;l<numloc-1;l++) {
	for (i=0;i<alpha[h][l].dim();i++) {	
	  cout << alpha[h][l][i] << " "; }
	cout << endl;
      }
    }
  }

  //check that all entries are between 0 and 1
  for (h=0;h<numhap;h++) {
    if (missing[h]==0) 
      for (l=0;l<numloc-1;l++) {
	for (i=0;i<alpha[h][l].dim();i++) {
	  if ( !(0.0<=alpha[h][l][i] && alpha[h][l][i]<1.0001) ) {
	    errorfs << "ERROR (hap est): alpha[" << h << "][" << l << "][" << i
		 << "] is not between 0 and 1 ["
		 << alpha[h][l][i] << "]" << endl;
	    exit(1); } } } }
  
  
}

void makealphamkv2_HfrBAYES(dvectorss &alpha, const ivector &missing,
			    const ivectors &a_lists,
			    const dvector &init,
			    const dvectors &trans, const dvectorss &b,
			    const int &debug)
{
  //makes alpha object, assuming haplotype controls

  int h,i,j,k,l;
  const int numhap=alpha.dim();
  const int numloc=a_lists.dim();

  //clear alpha
  for (h=0;h<numhap;h++) {
    if (missing[h]==0) {
      for (l=0;l<numloc-1;l++) {
	alpha[h][l].clear(0); } } }

  for (h=0;h<numhap;h++) {

    if (missing[h]==0) {
      //make first column
      for (i=0;i<alpha[h][0].dim();i++) {

	if (b[h][0][i] != 0) 
	  alpha[h][0][i]=
	    init(i)*b[h][0][i];
	
      }
    
    
      //make rest of alpha
      for (l=0;l<numloc-2;l++) {
      
	for (k=0;k<a_lists[l+2].dim();k++) {
	  for (j=0;j<a_lists[l+1].dim();j++) {
	    for (i=0;i<a_lists[l  ].dim();i++) {

	      alpha[h][l+1][j*a_lists[l+2].dim()+k]+=
		alpha[h][l  ][i*a_lists[l+1].dim()+j]*
		trans[l][
			 i*a_lists[l+1].dim()*a_lists[l+2].dim()+
			 j*a_lists[l+2].dim()+k]*
		b[h][l+1][j*a_lists[l+2].dim()+k];

	    }
	  }
	}
      }
    }
  }
  
  if (debug) {
    cout << "alpha [in maketransmkv2_HfrBAYES] " << endl;
    for (h=0;h<numhap;h++) {
      for (l=0;l<numloc-1;l++) {
	for (i=0;i<alpha[h][l].dim();i++) {	
	  cout << alpha[h][l][i] << " "; }
	cout << endl;
      }
    }
  }

  //check that all entries are between 0 and 1
  for (h=0;h<numhap;h++) {
    if (missing[h]==0) 
      for (l=0;l<numloc-1;l++) {
	for (i=0;i<alpha[h][l].dim();i++) {
	  if ( !(0.0<=alpha[h][l][i] && alpha[h][l][i]<1.0001) ) {
	    errorfs << "ERROR (hap est): alpha[" << h << "][" << l << "][" << i
		 << "] is not between 0 and 1 ["
		 << alpha[h][l][i] << "]" << endl;
	    exit(1); } } } }
  
  
}

void makebetamkv2_Hfr(dvectorss &beta, const ivector &missing, 
		      const dvectors &trans, const dvectorss &b,
		      const int &debug)
{
  //makes beta object, assuming haplotype controls

  int h,i,j,l;
  const int numhap=beta.dim();
  const int numloc=trans.dim()+2;

  //clear beta
  for (h=0;h<numhap;h++) {
    if (missing[h]==0) {
      for (l=0;l<numloc-1;l++) {
	beta[h][l].clear(0); } } }

  for (h=0;h<numhap;h++) {

    if (missing[h]==0) {
      //make last column
      for (i=0;i<beta[h][numloc-2].dim();i++) {
	
	beta[h][numloc-2][i]=1.0;
	
      }
      
      //make rest of beta
      for (l=numloc-3;l>=0;l--) {
	
	for (i=0;i<beta[h][l].dim();i++) {
	  
	  for (j=0;j<beta[h][l+1].dim();j++) {
	    if (b[h][l+1][j] != 0) { 
	      
	      beta[h][l][i]+=
		beta[h][l+1][j]*
		trans[l][i*beta[h][l+1].dim()+j]*
		b[h][l+1][j];
	    }
	  }
	}
	
      }
    }
  }

  if (debug) {
    cout << "beta [in makebetamkv2_Hfr] " << endl;
    for (h=0;h<numhap;h++) {
      for (l=0;l<numloc;l++) {
	for (i=0;i<beta[h][l].dim();i++) {	
	  cout << beta[h][l][i] << " "; }
	cout << endl;
      }
    }
  }

  //check that all entries are between 0 and 1
  for (h=0;h<numhap;h++) {
    if (missing[h]==0)
      for (l=0;l<numloc-1;l++) {
	for (i=0;i<beta[h][l].dim();i++) {
	  if ( !(0.0<=beta[h][l][i] && beta[h][l][i]<1.0001) ) {
	    errorfs << "ERROR (hap est): beta[" << h << "][" << l << "][" << i
		 << "] is not between 0 and 1 ["
		 << beta[h][l][i] << "]" << endl;
	    exit(1); } } } }
  
  
}

void makebetamkv2_HfrBAYES(dvectorss &beta, const ivector &missing, 
			   const ivectors &a_lists,
			   const dvectors &trans, const dvectorss &b,
			   const int &debug)
{
  //makes beta object, assuming haplotype controls

  int h,i,j,k,l;
  const int numhap=beta.dim();
  const int numloc=trans.dim()+2;

  //clear beta
  for (h=0;h<numhap;h++) {
    if (missing[h]==0) {
      for (l=0;l<numloc-1;l++) {
	beta[h][l].clear(0); } } }

  for (h=0;h<numhap;h++) {

    if (missing[h]==0) {
      //make last column
      for (i=0;i<beta[h][numloc-2].dim();i++) {
	
	beta[h][numloc-2][i]=1.0;
	
      }
      
      //make rest of beta
      for (l=numloc-3;l>=0;l--) {
	
	for (i=0;i<a_lists[l  ].dim();i++) {
	  for (j=0;j<a_lists[l+1].dim();j++) {
	    for (k=0;k<a_lists[l+2].dim();k++) {
	  
	      beta[h][l][i*a_lists[l+1].dim()+j]+=
		beta[h][l+1][j*a_lists[l+2].dim()+k]*
		trans[l][
			 i*a_lists[l+1].dim()*a_lists[l+2].dim()+
			 j*a_lists[l+2].dim()+k]*
		b[h][l+1][j*a_lists[l+2].dim()+k];
	    }
	  }
	}
	
      }
    }
  }

  if (debug) {
    cout << "beta [in makebetamkv2_HfrBAYES] " << endl;
    for (h=0;h<numhap;h++) {
      for (l=0;l<numloc;l++) {
	for (i=0;i<beta[h][l].dim();i++) {	
	  cout << beta[h][l][i] << " "; }
	cout << endl;
      }
    }
  }

  //check that all entries are between 0 and 1
  for (h=0;h<numhap;h++) {
    if (missing[h]==0)
      for (l=0;l<numloc-1;l++) {
	for (i=0;i<beta[h][l].dim();i++) {
	  if ( !(0.0<=beta[h][l][i] && beta[h][l][i]<1.0001) ) {
	    errorfs << "ERROR (hap est): beta[" << h << "][" << l << "][" << i
		 << "] is not between 0 and 1 ["
		 << beta[h][l][i] << "]" << endl;
	    exit(1); } } } }
  
  
}

void makeximkv2_Hfr(dvectorss &xi, const ivector &missing, 
		    const dvectorss &alpha, const dvectorss &beta, 
		    const dvectors &trans, const dvectorss &b, 
		    const int &debug)
{
  //makes xi object assuming haplotype controls

  int h,i,j,l;
  double denom;
  const int numhap=alpha.dim();
  const int numloc=trans.dim()+2;

  //clear xi
  for (h=0;h<numhap;h++) {
    if (missing[h]==0) {
      for (l=0;l<numloc-2;l++) {
	xi[h][l].clear(0); } } }

  //make xi
  for (h=0;h<numhap;h++) {

    if (missing[h]==0) {
      for (l=0;l<numloc-2;l++) {  
	for (i=0;i<alpha[h][l].dim();i++) {

	  for (j=0;j<alpha[h][l+1].dim();j++) {
	
	    if (b[h][l+1][j] != 0)
	      xi[h][l][i*alpha[h][l+1].dim()+j]=
		alpha[h][l][i]*
		trans[l][i*alpha[h][l+1].dim()+j]*
		b[h][l+1][j]*
		beta[h][l+1][j];
	    
	  }
	}
      }
      //make denominator
      denom=0.0;
      for (i=0;i<xi[h][0].dim();i++) denom+=xi[h][0][i];
      if (denom<=0) {
	errorfs << "ERROR: in makeximkv2_Hfr,denominator for xi at haplotype " 
	     << h << "(+1) <=0 [" << denom << "] " << endl;
	exit(1); }
      //divide all entries by this number
      for (l=0;l<numloc-2;l++) {  
	for (i=0;i<xi[h][l].dim();i++) { 
	  xi[h][l][i]/=denom;
	}
      }

    }
  }

  if (debug) {
    cout << "xi [in makeximkv2_Hfr] " << endl;
    for (h=0;h<numhap;h++) {
      for (l=0;l<numloc-2;l++) {
	for (i=0;i<xi[h][l].dim();i++) {	
	  cout << xi[h][l][i] << " "; }
	cout << endl;
      }
    }
  }

  //check that all entries sum to 1 at each locus for each haplotype
  for (h=0;h<numhap;h++) {
    for (l=0;l<numloc-2;l++) {  
      denom=0.0;
      for (i=0;i<xi[h][l].dim();i++) denom+=xi[h][l][i];
      if ( !(0.9999<denom && denom<1.0001) ) {
	errorfs << "ERROR: xi[" << h << "][" << l 
		<< "][] does not sum to 1.0 ["
		<< denom << "]" << endl;
      }
    }
  }

  
}

void makeximkv2_HfrBAYES(dvectorss &xi, const ivector &missing, 
			 const ivectors &a_lists,
			 const dvectorss &alpha, const dvectorss &beta, 
			 const dvectors &trans, const dvectorss &b, 
			 const int &debug)
{
  //makes xi object assuming haplotype controls

  int h,i,j,k,l;
  double denom;
  const int numhap=alpha.dim();
  const int numloc=trans.dim()+2;

  //clear xi
  for (h=0;h<numhap;h++) {
    if (missing[h]==0) {
      for (l=0;l<numloc-2;l++) {
	xi[h][l].clear(0); } } }

  //make xi
  for (h=0;h<numhap;h++) {

    if (missing[h]==0) {
      for (l=0;l<numloc-2;l++) {  

	for (i=0;i<a_lists[l  ].dim();i++) {
	  for (j=0;j<a_lists[l+1].dim();j++) {
	    for (k=0;k<a_lists[l+2].dim();k++) {
	      
	      if (b[h][l+1][j*a_lists[l+2].dim()+k]>0)
		xi[h][l][
			 (i*a_lists[l+1].dim()+j)*alpha[h][l+1].dim()+
			 (j*a_lists[l+2].dim()+k)]=
		  alpha[h][l][i*a_lists[l+1].dim()+j]*
		  trans[l][
			   i*a_lists[l+1].dim()*a_lists[l+2].dim()+
			   j*a_lists[l+2].dim()+k]*
		  b[h][l+1][j*a_lists[l+2].dim()+k]*
		  beta[h][l+1][j*a_lists[l+2].dim()+k];
	      
	    }
	  }
	}
      }
      //make denominator
      denom=0.0;
      for (i=0;i<xi[h][0].dim();i++) denom+=xi[h][0][i];
      if (denom<=0) {
	errorfs << "ERROR: in makeximkv2_Hfr,denominator for xi at haplotype " 
	     << h << "(+1) <=0 [" << denom << "] " << endl;
	exit(1); }
      //divide all entries by this number
      for (l=0;l<numloc-2;l++) {  
	for (i=0;i<xi[h][l].dim();i++) { 
	  xi[h][l][i]/=denom;
	}
      }

    }
  }

  if (debug) {
    cout << "xi [in makeximkv2_HfrBAYES] " << endl;
    for (h=0;h<numhap;h++) {
      for (l=0;l<numloc-2;l++) {
	for (i=0;i<xi[h][l].dim();i++) {	
	  cout << xi[h][l][i] << " "; }
	cout << endl;
      }
    }
  }

  //check that all entries sum to 1 at each locus for each haplotype
  for (h=0;h<numhap;h++) {
    for (l=0;l<numloc-2;l++) {  
      denom=0.0;
      for (i=0;i<xi[h][l].dim();i++) denom+=xi[h][l][i];
      if ( !(0.9999<denom && denom<1.0001) ) {
	errorfs << "ERROR: xi[" << h << "][" << l 
		<< "][] does not sum to 1.0 ["
		<< denom << "]" << endl;
      }
    }
  }

  
}

void makeh3freqlistsHfrxi(dvectors &h3_freq_lists, 
			  const dvectorss &xi, 
			  const ivectors &a_lists, const int &debug)
{
  //makes haplotype frequencies from xi, assuming haplotype controls

  int h,i,j,k,l;
  double value;

  const int numhap=xi.dim();
  const int numloc=h3_freq_lists.dim()+2;

  //clear h3_freq_lists
  for (l=0;l<numloc-2;l++) {
    h3_freq_lists[l].clear(0); } 

  for (h=0;h<numhap;h++) {
    for (l=0;l<numloc-2;l++) {
      
      for (i=0;i<a_lists[l  ].dim();i++) {
	for (j=0;j<a_lists[l+1].dim();j++) {
	  for (k=0;k<a_lists[l+2].dim();k++) {

	  h3_freq_lists[l][
			  i*a_lists[l+1].dim()*a_lists[l+2].dim()+
			  j*a_lists[l+2].dim()+k]

	    +=xi[h][l][
		       (i*a_lists[l+1].dim()+j)*
		       a_lists[l+1].dim()*a_lists[l+2].dim()+
		       (j*a_lists[l+2].dim()+k) ];
	      
	  }
	}
      }
    }
  }

  //divide by numhap
  for (l=0;l<numloc-2;l++) {
    for (i=0;i<h3_freq_lists[l].dim();i++) {
      h3_freq_lists[l][i]/=numhap; } }
  
  if (debug) {
    cout << "h3_freq_lists [in makeximkv2_Hfr] " << endl;
    for (l=0;l<numloc-2;l++) {
      for (i=0;i<h3_freq_lists[l].dim();i++) {	
	cout << h3_freq_lists[l][i] << " "; }
      cout << endl;
    }
  }

  //check that frequencies sum to 1 for each pair
  for (l=0;l<numloc-2;l++) {
    value=0.0;
    for (i=0;i<h3_freq_lists[l].dim();i++) value+=h3_freq_lists[l][i];
    if ( !(0.9999<value && value<1.0001) ) {
      errorfs << "ERROR in makeh3freqlistsHfrxi: h3_freq_lists[" 
	      << l << "] does not sum to 1 ["
	      << value << "]" << endl;
      exit(1); } }

  
}

void makeh3countsHfrxi(dvectors &h3_freq_lists, 
		       const dvectorss &xi, 
		       const ivectors &a_lists, const int &debug)
{
  //makes haplotype frequencies from xi, assuming haplotype controls

  int h,i,j,k,l;
  double value;

  const int numhap=xi.dim();
  const int numloc=h3_freq_lists.dim()+2;

  //clear h3_freq_lists
  for (l=0;l<numloc-2;l++) {
    h3_freq_lists[l].clear(0); } 

  for (h=0;h<numhap;h++) {
    for (l=0;l<numloc-2;l++) {
      
      for (i=0;i<a_lists[l  ].dim();i++) {
	for (j=0;j<a_lists[l+1].dim();j++) {
	  for (k=0;k<a_lists[l+2].dim();k++) {

	  h3_freq_lists[l][
			  i*a_lists[l+1].dim()*a_lists[l+2].dim()+
			  j*a_lists[l+2].dim()+k]

	    +=xi[h][l][
		       (i*a_lists[l+1].dim()+j)*
		       a_lists[l+1].dim()*a_lists[l+2].dim()+
		       (j*a_lists[l+2].dim()+k) ];
	      
	  }
	}
      }
    }
  }

  //divide by numhap
  //for (l=0;l<numloc-2;l++) {
  //for (i=0;i<h3_freq_lists[l].dim();i++) {
  //h3_freq_lists[l][i]/=numhap; } }
  
  if (debug) {
    cout << "h3_freq_lists [in makeximkv2_Hfr] " << endl;
    for (l=0;l<numloc-2;l++) {
      for (i=0;i<h3_freq_lists[l].dim();i++) {	
	cout << h3_freq_lists[l][i] << " "; }
      cout << endl;
    }
  }

  //check that frequencies sum to numhap for each pair
  for (l=0;l<numloc-2;l++) {
    value=0.0;
    for (i=0;i<h3_freq_lists[l].dim();i++) value+=h3_freq_lists[l][i];
    if ( !(numhap-0.0001<value && value<numhap+0.0001) ) {
      errorfs << "ERROR in makeh3countsHfrxi: h3_freq_lists[" 
	      << l << "] does not sum to numhap ["
	      << value << "]" << endl;
      exit(1); } }

  
}

void makeh3freqlistsmkv2(dvectors &h3_freq_lists, dvectors &h2_freq_lists,
			 const dvector &init,
			 const dvectors &trans,const ivectors &a_lists,
			 const int &debug)
{
  int l,i,j,k;
  const int numloc=a_lists.dim();
  double value=0.0;
  
  for (l=0;l<numloc-2;l++) {
    h2_freq_lists[l].clear(0.0);
    h3_freq_lists[l].clear(0.0); }
  h2_freq_lists[numloc-2].clear(0.0);

  h2_freq_lists[0]=init;
  //cout << "init" << " ";
  //for (i=0;i<h2_freq_lists[0].dim();i++) 
  //cout << init(i) << " "; cout << endl;

  for (l=0;l<numloc-2;l++) {

    for (i=0;i<a_lists[l  ].dim();i++) {
      for (j=0;j<a_lists[l+1].dim();j++) {
	for (k=0;k<a_lists[l+2].dim();k++) {

	  h3_freq_lists[l][
			   i*a_lists[l+1].dim()*a_lists[l+2].dim()+
			   j*a_lists[l+2].dim()+
			   k ]=
	    h2_freq_lists[l][ i*a_lists[l+1].dim() +j]*
	    trans[l][
		     (i*a_lists[l+1].dim() +j)*
		     a_lists[l+1].dim()*a_lists[l+2].dim()+
		     (j*a_lists[l+2].dim() +k)];

	  //update h2_freq_lists[l+1]
	  h2_freq_lists[l+1][ j*a_lists[l+2].dim() +k]+=
	    h3_freq_lists[l][
			     i*a_lists[l+1].dim()*a_lists[l+2].dim()+
			     j*a_lists[l+2].dim()+
			     k ];
	}
      }
    }
  }
  if (debug) {
    cout << "h3_freq_lists [in makeh3freqlistsmkv2] " << endl;
    for (l=0;l<numloc-2;l++) {
      for (i=0;i<h3_freq_lists[l].dim();i++) {	
	cout << h3_freq_lists[l][i] << " "; }
      cout << endl;
    }
  }

  //check that frequencies sum to 1 for each pair
  for (l=0;l<numloc-2;l++) {
    value=0.0;
    for (i=0;i<h3_freq_lists[l].dim();i++) value+=h3_freq_lists[l][i];
    if ( !(0.9999<value && value<1.0001) ) {
      errorfs << "ERROR in makeh3freqlistsmkv2: h3_freq_lists[" 
	      << l << "] does not sum to 1 ["
	      << value << "]" << endl;
      exit(1); } }


}

void makeh2freqlistsmkv2(dvectors &h2_freq_lists,
			 const dvectors &h3_freq_lists, 
			 const ivectors &a_lists, const int &debug)
{
  int i,j,k,l;
  double value=0.0;
  const int numloc=a_lists.dim();

  for (l=0;l<numloc-1;l++) {
    h2_freq_lists[l].clear(0.0);
  }

  for (l=0;l<numloc-2;l++) {

    for (i=0;i<a_lists[l  ].dim();i++) {
      for (j=0;j<a_lists[l+1].dim();j++) {
	for (k=0;k<a_lists[l+2].dim();k++) {

	  if (l==0)
	    h2_freq_lists[l][ i*a_lists[l+1].dim() +j]+=
	      h3_freq_lists[l][
			       i*a_lists[l+1].dim()*a_lists[l+2].dim()+
			       j*a_lists[l+2].dim()+
			       k ];
	  
	  h2_freq_lists[l+1][ j*a_lists[l+2].dim() +k]+=
	    h3_freq_lists[l][
			     i*a_lists[l+1].dim()*a_lists[l+2].dim()+
			     j*a_lists[l+2].dim()+
			     k ];

	}
      }
    }
  }
    
  if (debug) {
    cout << "h2_freq_lists " << endl;
    for (l=0;l<numloc-1;l++) {
      for (i=0;i<h2_freq_lists[l].dim();i++) {	
	cout << h2_freq_lists[l][i] << " "; }
      cout << endl;
    }
  }

  //check that frequencies sum to 1 for each pair
  for (l=0;l<numloc-1;l++) {
    value=0.0;
    for (i=0;i<h2_freq_lists[l].dim();i++) value+=h2_freq_lists[l][i];
    if ( !(0.9999<value && value<1.0001) ) {
      errorfs << "ERROR: h2_freq_lists[" << l << "] does not sum to 1 ["
	   << value << "]" << endl;
      exit(1); } }
}

//adjustment function

void adjustmkv2(dvectors &h3_freq_lists_adj, const dvectors& h3_freq_lists,
		dvectors &h2_freq_lists, const double &e,
		const ivectors &a_lists, const int &debug)
{

  //this function adjust for joint frequencies with estimated freq. of 0

  int i,j,k,l;
  const int numloc=a_lists.dim();
  int place;
  int ct=0;
  double den=0.0;

  //initialize h3_freq_lists_adj
  h3_freq_lists_adj.init(numloc-2);
  for (l=0;l<h3_freq_lists_adj.dim();l++) {
    h3_freq_lists_adj[l].init(a_lists[l].dim()*
			      a_lists[l+1].dim()*a_lists[l+2].dim());
    h3_freq_lists_adj[l].clear(0); }

  //copy h3_freq_lists into h3_freq_lists_adj
  for (l=0;l<h3_freq_lists_adj.dim();l++) {
    h3_freq_lists_adj[l]=h3_freq_lists[l]; }

  //translate into hmm objects

  //declare and initialize objects
  dvector init(a_lists[0].dim()); init.clear(0);
  dvector trans(a_lists[0].dim()*a_lists[1].dim()); trans.clear(0.0);  
  dvector init2(trans.dim()); init2.clear(0.0); 
  dvectors trans2(numloc-2);
  for (l=0;l<numloc-2;l++) {
    trans2[l].init(a_lists[l].dim()*a_lists[l+1].dim()*
		   a_lists[l+1].dim()*a_lists[l+2].dim());
    trans2[l].clear(0); }

  //fill in init (assume no need to adjust initial frequencies)
  for (i=0;i<a_lists[0].dim();i++) {
    for (j=0;j<a_lists[1].dim();j++) {
      init(i)+=h2_freq_lists[0][i*a_lists[1+1].dim()+j]; 
} }

  //fill in trans
  for (i=0;i<a_lists[0].dim();i++) {
    for (j=0;j<a_lists[1].dim();j++) {
      trans(i*a_lists[1].dim()+j)=h2_freq_lists[0][i*a_lists[l].dim()+j]/
	init(i); } }
  //cout << "trans (pre) ";
  //for (i=0;i<trans.dim();i++) cout << trans(i) << " "; cout << endl;
  //check that trans sums to 1 for each i

  //fill in trans2
  maketransmkv2_Hfr(trans2, h3_freq_lists,h2_freq_lists,a_lists,0);

  //adjust trans
  for (i=0;i<a_lists[0].dim();i++) {
    ct=0;
    den=0.0;
    for (j=0;j<a_lists[1].dim();j++) {
      if (trans(i*a_lists[1].dim()+j)<e) ct+=1; 
      else den+=trans(i*a_lists[1].dim()+j);
    }
    if (ct==a_lists[1].dim()) 
      for (j=0;j<a_lists[1].dim();j++) trans(i*a_lists[1].dim()+j)=1.0/ct;
    else { 
      for (j=0;j<a_lists[1].dim();j++) {
	if (trans(i*a_lists[1].dim()+j)<e) 
	  trans(i*a_lists[1].dim()+j)=e/(den+ct*e);
	else trans(i*a_lists[1].dim()+j)/=(den+ct*e);
      } }
  }
  //cout << "trans (post) ";
  //for (i=0;i<trans.dim();i++) cout << trans(i) << " "; cout << endl;
  //check that trans still sums to 1
  //cout << trans(0)+trans(1) << " " << trans(2)+trans(3) << endl;
  //exit(1);

  //update trans2
  for (l=0;l<numloc-2;l++) {
    for (i=0;i<a_lists[l].dim();i++) {
      for (j=0;j<a_lists[l+1].dim();j++) {
	ct=0;
	den=0.0;
	for (k=0;k<a_lists[l+2].dim();k++) {
	  place=
	    (i*a_lists[l+1].dim()+j)*h2_freq_lists[l+1].dim()+
	    (j*a_lists[l+2].dim()+k);
	  if (trans2[l][place]<e) ct+=1;
	  else den+=trans2[l][place];
	}
	if (ct==h2_freq_lists[l+1].dim())
	  for (k=0;k<a_lists[l+2].dim();k++) {
	    place=
	      (i*a_lists[l+1].dim()+j)*h2_freq_lists[l+1].dim()+
	      (j*a_lists[l+2].dim()+k);
	  trans2[l][place]=1.0/ct; }
      else {
	for (k=0;k<a_lists[l+2].dim();k++) {
	  place=
	    (i*a_lists[l+1].dim()+j)*h2_freq_lists[l+1].dim()+
	    (j*a_lists[l+2].dim()+k);
	  if (trans2[l][place]<e)
	    trans2[l][place]=e/(den+ct*e);
	  else trans2[l][place]/=(den+ct*e);
	} }
      }
    }
  }

  //update two-locus hap freq.
  for (i=0;i<a_lists[0].dim();i++) {
    for (j=0;j<a_lists[1].dim();j++) {
      init2(i*a_lists[1].dim()+j)=
	init(i)*trans(i*a_lists[1+1].dim()+j); } }
  //cout << "init2 ";
  //for (i=0;i<init2.dim();i++) cout << init2(i) << " "; cout << endl;
  
  makeh3freqlistsmkv2(h3_freq_lists_adj,h2_freq_lists,init2,trans2,a_lists,0);
  makeh2freqlistsmkv2(h2_freq_lists,h3_freq_lists_adj,a_lists,0);

  if (debug) {
    cout << "h3_freq_lists_adj" << endl;
    for (l=0;l<numloc-2;l++) {
      for (i=0;i<h3_freq_lists_adj[l].dim();i++) 
	cout << h3_freq_lists_adj[l][i] << " ";
      cout << endl;
    }
  }

}

//genotype version
void makehap3fr_hmm_GconBAYES(dvectors &h3_freq_lists, 
			      const ivectors &a_lists, 
			      const igrid &controls, const int &bayes, 
			      const int &debug)
{

  //estimates three-marker haplotype frequencies, using HMM, assuming
  //genotype controls

  //note that bayes=0 makes this the EM
  int l,g,i,j,k,iter,conv;
  double diff,sum;
  const int numloc=controls.dimy();
  const int numgen=controls.dimx()/2;
  //double llik=0.0,lik=0.0;

  //declare and intialize and clear lists of two-locus hap freqs 
  //corresponding to a_lists
  dvectors h2_freq_lists; 
  h2_freq_lists.init(numloc-1);
  for (l=0;l<numloc-1;l++) 
    { h2_freq_lists[l].init(a_lists[l].dim()*a_lists[l+1].dim()); 
    h2_freq_lists[l].clear(0); }
  
  //declare and initialize and clear hmm objects
  //states are two locus haplotypes
  dvectorss b(numgen);
  dvectorss alpha(numgen),beta(numgen),gamma(numgen),xi(numgen);
  for (g=0;g<numgen;g++) {
    b[g].init(numloc-1);
    alpha[g].init(numloc-1);
    beta[g].init(numloc-1);
    gamma[g].init(numloc-1);
    for (l=0;l<numloc-1;l++) {
      b[g][l].init(a_lists[l].dim()*a_lists[l+1].dim()*
		   a_lists[l].dim()*a_lists[l+1].dim()); 
      b[g][l].clear(0);
      alpha[g][l].init(b[g][l].dim()); alpha[g][l].clear(0);
      beta[g][l].init(b[g][l].dim()); beta[g][l].clear(0);
    }
  }
  //initialize xi vectors for all genotypes
  for (g=0;g<numgen;g++) {
    xi[g].init(numloc-2);
    for (l=0;l<numloc-2;l++) {
      xi[g][l].init(a_lists[l].dim()*a_lists[l+1].dim()*
		    a_lists[l+2].dim()*
		    a_lists[l].dim()*a_lists[l+1].dim()*
		    a_lists[l+2].dim());
      xi[g][l].clear(0.0); }
  }
  
  dvector init;
  init.init( a_lists[0].dim()*a_lists[1].dim() );
  init.clear(0.0);

  //make copy to determine convergence
  dvector initC; initC.init( init.dim() ); initC.clear(0.0);

  dvectors trans(numloc-2);
  for (l=0;l<numloc-2;l++) {
    trans[l].init(a_lists[l].dim()*a_lists[l+1].dim()*
		  a_lists[l+2].dim());
    trans[l].clear(0); }
  
  //make copy of trans to determine convergence
  dvectors transC(numloc-2);
  for (l=0;l<numloc-2;l++) {
    transC[l].init(trans[l].dim()); transC[l].clear(0); }
  
  //assign initial values to parameters
  for (i=0;i<init.dim();i++) 
    init(i)=1.0/ ( init.dim() );
  for (l=0;l<numloc-2;l++) {
    for (i=0;i<a_lists[l].dim();i++) {
      for (j=0;j<a_lists[l+1].dim();j++) {
	for (k=0;k<a_lists[l+2].dim();k++) {
	  trans[l][ 
		   i*a_lists[l+1].dim()*a_lists[l+2].dim()+
		   j*a_lists[l+2].dim()+k ]=
	    1.0/( a_lists[l+2].dim() ); } } } }
  
  //make h3_freq_lists from init,trans
  //makeh3freqlistsmkv2(h3_freq_lists,h2_freq_lists,init,trans,a_lists,0);
  //makeh2freqlistsmkv2(h2_freq_lists,h3_freq_lists,a_lists,0);

  //make observation probabilities [b]
  makebmkv2_Gfr(b, a_lists, controls, 0);

  //delete me
  //cout << "b ";
  //for (i=0;i<b[75][0].dim();i++) cout << b[75][0][i] << " "; cout << endl;

  //iterate hmm
  int maxiter=500;
  for (iter=1;iter<=maxiter;iter++) {

    //copy init to use later in checking convergence
    initC=init;

    //copy hapl. freq. to use later in checking convergence
    for (l=0;l<transC.dim();l++) 
      transC[l]=trans[l];
    
    //make trans (from h3_freq_lists)
    //maketransmkv2_Hfr(trans, h3_freq_lists, h2_freq_lists, a_lists,0);

    //make alpha
    makealphamkv2_Gfr(alpha, a_lists, init, trans, b, 0);

    //make beta
    makebetamkv2_Gfr(beta, a_lists, trans, b, 0);

    //make xi
    makeximkv2_Gfr(xi, a_lists, alpha, beta, trans, b, 0);

    //compute likelihood from HMM
    //cout << "iter " << iter << " ";
    //llik=0; 
    //for (h=0;h<numhap;h++) {
    //lik=0.0;
    //for (i=0;i<alpha[h][0].dim();i++) lik+=alpha[h][0][i]*beta[h][0][i];
    //llik+=log(lik);
    //}
    //cout << "llik " << llik << endl;

    //give expected counts of 3-locus haplotypes
    makeh3countsGfrxi(h3_freq_lists, xi, a_lists, 0);
    //from this make expected counts of two-locus haplotypes
    makeh2countsGmkv2(h2_freq_lists,numgen,h3_freq_lists,a_lists,0);

    //cout << "iter " << iter << endl;
      // << h3_freq_lists[70][4] << " " << h3_freq_lists[70][5] << " "
      // << h2_freq_lists[70][2] << " " 
      // << h3_freq_lists[70][4]/h2_freq_lists[70][2] << " "
      // << h3_freq_lists[70][5]/h2_freq_lists[70][2] << endl;


    //Bayesian adjustment/update to init
    init=h2_freq_lists[0];
    for (i=0;i<init.dim();i++) { 
      init(i)+=1.0*bayes;
      init(i)/=2.0*numgen+1.0*init.dim()*bayes; }

    //check that init sums to 1
    sum=0;
    for (i=0;i<init.dim();i++) sum+=init(i);
    if ( !(0.9999<sum && sum<1.0001) ) {
      errorfs << "ERROR: Bayesian adjustment to init in gen. freq. est.; "
	      << "sums to " << sum << endl;
      exit(1); }

    //delete me
    //cout << "init ";
    //for (i=0;i<init.dim();i++) cout << init(i) << " "; cout << endl;

    //Bayesian adjustment/update to trans
    for (l=0;l<numloc-2;l++) {
      for (i=0;i<a_lists[l  ].dim();i++) {
	for (j=0;j<a_lists[l+1].dim();j++) {
	  for (k=0;k<a_lists[l+2].dim();k++) {
	    if (h3_freq_lists[l][i*a_lists[l+1].dim()*a_lists[l+2].dim()+
				j*a_lists[l+2].dim()+k]  +1.0*bayes <1e-100
		&&
		h2_freq_lists[l][i*a_lists[l+1].dim()+j]+ 
		1.0*a_lists[l+2].dim()*bayes<1e-100)
	      trans[l][i*a_lists[l+1].dim()*a_lists[l+2].dim()+
		      j*a_lists[l+2].dim()+k]=0.0;
	    else
	      trans[l][i*a_lists[l+1].dim()*a_lists[l+2].dim()+
		      j*a_lists[l+2].dim()+k]=
		( h3_freq_lists[l][i*a_lists[l+1].dim()*a_lists[l+2].dim()+
				  j*a_lists[l+2].dim()+k]  +1.0*bayes )/
		( h2_freq_lists[l][i*a_lists[l+1].dim()+j]+ 
		  1.0*a_lists[l+2].dim()*bayes);
	  }
	  //if they're all 0, set first to 1
	  conv=0;
	  for (k=0;k<a_lists[l+2].dim();k++) 
	    if (trans[l][i*a_lists[l+1].dim()*a_lists[l+2].dim()+
			j*a_lists[l+2].dim()+k]>0) { conv=1; break; }
	  if (conv==0) trans[l][i*a_lists[l+1].dim()*a_lists[l+2].dim()+
			j*a_lists[l+2].dim()  ]=1;
	}
      }
    }//l

    //delete me
    //  cout << "iter " << iter << " "
    //< trans[70][4] << " " << trans[70][5] << endl;

     //check trans
    for (l=0;l<numloc-2;l++) {
      for (i=0;i<a_lists[l  ].dim();i++) {
	for (j=0;j<a_lists[l+1].dim();j++) {
	  sum=0;
	  for (k=0;k<a_lists[l+2].dim();k++) {
	    sum+=trans[l][
			  i*a_lists[l+1].dim()*a_lists[l+2].dim()+
			  j*a_lists[l+2].dim()+k]; }
	  if ( !(0.9999<sum && sum<1.0001) ) {

	    for (k=0;k<a_lists[l+2].dim();k++) 
	      errorfs << trans[l][
				  i*a_lists[l+1].dim()*a_lists[l+2].dim()+
				  j*a_lists[l+2].dim()+k] << " ";
	    errorfs << endl;
	    errorfs << "ERROR in Bayes adj. of trans. Doesn't sum to 1 for "
		    << "i " << i << " j " << j << " pos "
		    << i*a_lists[l+1].dim()*a_lists[l+2].dim()+
	      j*a_lists[l+2].dim() 
		    << " at locus " << l << " ( " << sum << " ), iter " 
		    << iter << endl;
	    
	    errorfs << trans[70][4] << " " << trans[70][5] << endl;
	    exit(1); }}}}
    
    //delete me
    //cout << "trans[0] ";
    //for (i=0;i<trans[0].dim();i++) cout << trans[0][i] << " "; cout << endl;
    //cout << "init*trans[0] "
    // << init(0)*trans[0][0] +init(2)*trans[0][4] << " "
    // << init(0)*trans[0][1] +init(2)*trans[0][5] << " "
    // << init(1)*trans[0][2] +init(3)*trans[0][6] << " "
    // << init(1)*trans[0][3] +init(3)*trans[0][7] << endl;

    //check convergence
    conv=0;
    for (i=0;i<init.dim();i++) {
      diff=init(i)-initC(i);
      if ( diff<-1e-8 || 1e-8<diff ) { conv=1; break; } }
    if (conv==0) {
      for (l=0;l<numloc-2;l++) {
	for (i=0;i<trans[l].dim();i++) {
	  diff=trans[l][i]-transC[l][i];
	  if ( diff<-1e-8 || 1e-8<diff ) { conv=1; break; } }
	if (conv==1) break; } }
    if (conv==0) break; 

  }//iter

  //make h3_freq_lists from trans, init
  makeh3freqlistsGmkv2(h3_freq_lists,h2_freq_lists,init,trans,a_lists,0);

  if (debug) {
    cout << "init " << endl;
    for (i=0;i<init.dim();i++) cout << init(i) << " "; cout << endl;
    cout << "trans " << endl;
    for (l=0;l<numloc-2;l++) {
      for (i=0;i<trans[l].dim();i++) {
	cout << trans[l][i] << " "; }
      cout << endl; }

    cout << "h3_freq_lists [" << iter << " iterations]" << endl;
    for (l=0;l<numloc-2;l++) {
      for (i=0;i<h3_freq_lists[l].dim();i++) {	
	cout << h3_freq_lists[l][i] << " "; }
      cout << endl;
    }
  }
}

void makebmkv2_Gfr(dvectorss &b, 
		   const ivectors &a_lists,
		   const igrid &controls, const int &debug)
{
  //makes observation probabilities, assuming haplotype data

  int i,j,g,k,m,l,ind=0;

  const int numgen=controls.dimx()/2;
  const int numloc=controls.dimy();

  for (g=0;g<numgen;g++) {

    for (l=0;l<numloc-1;l++) {
      
      for (i=0;i<a_lists[l].dim();i++) {
	for (j=0;j<a_lists[l].dim();j++) { 
	  for (k=0;k<a_lists[l+1].dim();k++) {
	    for (m=0;m<a_lists[l+1].dim();m++) { 
	      
	      ind=0;
	      if ( controls(2*g  ,l)==a_lists[l][i] || 
		   controls(2*g  ,l)==0 )
		if ( controls(2*g+1,l)==a_lists[l][j] || 
		     controls(2*g+1,l)==0 ) {
		  
		  if ( controls(2*g  ,l+1)==a_lists[l+1][k] || 
		       controls(2*g  ,l+1)==0 )
		    if ( controls(2*g+1,l+1)==a_lists[l+1][m] || 
			 controls(2*g+1,l+1)==0 )
		      ind=1;
		  
		  if ( controls(2*g  ,l+1)==a_lists[l+1][m] || 
		       controls(2*g  ,l+1)==0 )
		    if ( controls(2*g+1,l+1)==a_lists[l+1][k] || 
			 controls(2*g+1,l+1)==0 )
		      ind=1; }
		  
	      if ( controls(2*g  ,l)==a_lists[l][j] || 
		   controls(2*g  ,l)==0 )
		if ( controls(2*g+1,l)==a_lists[l][i] || 
		     controls(2*g+1,l)==0 ) {
		  
		  if ( controls(2*g  ,l+1)==a_lists[l+1][k] || 
		       controls(2*g  ,l+1)==0 )
		    if ( controls(2*g+1,l+1)==a_lists[l+1][m] || 
			 controls(2*g+1,l+1)==0 )
		      ind=1;
		  
		  if ( controls(2*g  ,l+1)==a_lists[l+1][m] || 
		       controls(2*g  ,l+1)==0 )
		    if ( controls(2*g+1,l+1)==a_lists[l+1][k] || 
			 controls(2*g+1,l+1)==0 )
		      ind=1; }

	      if (ind==1)
		b[g][l][
			i*
			a_lists[l].dim()*
			a_lists[l+1].dim()*
			a_lists[l+1].dim()+
			j*
			a_lists[l+1].dim()*
			a_lists[l+1].dim()+
			m*
			a_lists[l+1].dim()+
			k
		]=1.0;

	    }
	  }
	}
      }
    }
  }

  if (debug) {
    g=0;
    cout << "b [obs. probs. for genotype " << g << "]" << endl;
    for (l=0;l<numloc-1;l++) {
      for (i=0;i<b[g][l].dim();i++) {
	cout << b[g][l][i] << " "; }
      cout << endl;
    }
    
  }
  
}

void makealphamkv2_Gfr(dvectorss &alpha, 
		       const ivectors &a_lists,
		       const dvector &init,
		       const dvectors &trans, const dvectorss &b,
		       const int &debug)
{
  //makes alpha object, assuming haplotype controls
  
  int g,i,j,k,l,m,n,o,place;
  const int numgen=alpha.dim();
  const int numloc=a_lists.dim();
  
  //clear alpha
  for (g=0;g<numgen;g++) {
    for (l=0;l<numloc-1;l++) {
      alpha[g][l].clear(0); } }

  for (g=0;g<numgen;g++) {

    //make first column
    for (i=0;i<a_lists[0].dim();i++) {
      for (j=0;j<a_lists[0].dim();j++) {
	for (k=0;k<a_lists[1].dim();k++) {
	  for (m=0;m<a_lists[1].dim();m++) {
	    
	    place=
	      i*a_lists[0].dim()*a_lists[1].dim()*a_lists[1].dim()+
	      j*a_lists[1].dim()*a_lists[1].dim()+
	      k*a_lists[1].dim()+
	      m;
	    alpha[g][0][place]=
	      init(i*a_lists[1].dim()+k)*init(j*a_lists[1].dim()+m)*
	      b[g][0][place]; 
	  } } } }
    
    //make rest of alpha
    for (l=0;l<numloc-2;l++) {
      
      for (k=0;k<a_lists[l+1].dim();k++) {
	for (m=0;m<a_lists[l+1].dim();m++) {
	  for (n=0;n<a_lists[l+2].dim();n++) {
	    for (o=0;o<a_lists[l+2].dim();o++) {

	      for (i=0;i<a_lists[l].dim();i++) {
		for (j=0;j<a_lists[l].dim();j++) {

		  place=
		    k*
		    a_lists[l+1].dim()*
		    a_lists[l+2].dim()*
		    a_lists[l+2].dim()+
		    m*
		    a_lists[l+2].dim()*
		    a_lists[l+2].dim()+
		    n*
		    a_lists[l+2].dim()+
		    o;
				  
		  if (b[g][l+1][place]>0)
		    alpha[g][l+1][place]+=
		      alpha[g][l][
				  i*
				  a_lists[l  ].dim()*
				  a_lists[l+1].dim()*
				  a_lists[l+1].dim()+
				  j*
				  a_lists[l+1].dim()*
				  a_lists[l+1].dim()+
				  k*
				  a_lists[l+1].dim()+
				  m
		      ]*
		      trans[l][
			       i*a_lists[l+1].dim()*a_lists[l+2].dim()+
			       k*a_lists[l+2].dim()+
			       n]*
		      trans[l][
			       j*a_lists[l+1].dim()*a_lists[l+2].dim()+
			       m*a_lists[l+2].dim()+
			       o]*
		      b[g][l+1][place];
		}
	      }
	    }
	  }
	}
      }
    }
  }

  if (debug) {
    g=0;
    cout << "alpha [for genotype " << g << "]" << endl;
    for (l=0;l<numloc-1;l++) {
      for (i=0;i<alpha[g][l].dim();i++) {
	cout << alpha[g][l][i] << " "; }
      cout << endl;
    }
  }

  //check that all entries are between 0 and 1
  for (g=0;g<numgen;g++) {
    for (l=0;l<numloc-1;l++) {
      for (i=0;i<alpha[g][l].dim();i++) {
	if ( !(0.0<=alpha[g][l][i] && alpha[g][l][i]<1.0001) ) {
	    errorfs << "ERROR (hap est): alpha[" << g << "][" << l << "][" << i
		    << "] is not between 0 and 1 ["
		    << alpha[g][l][i] << "]" << endl;
	    exit(1); } } } }

}
  
  
void makebetamkv2_Gfr(dvectorss &beta,
		      const ivectors &a_lists,
		      const dvectors &trans, const dvectorss &b,
		      const int &debug)
{
  //makes beta object, assuming haplotype controls
  
  int g,i,j,k,l,m,n,o,place;
  const int numgen=beta.dim();
  const int numloc=a_lists.dim();
  
  //clear beta
  for (g=0;g<numgen;g++) {
    for (l=0;l<numloc-1;l++) {
      beta[g][l].clear(0.0); } }

  for (g=0;g<numgen;g++) {

    //make last column of beta
    for (i=0;i<beta[g][numloc-2].dim();i++) {
      beta[g][numloc-2][i]=1.0; }

    //make rest of beta
    for (l=numloc-3;l>=0;l--) {
      
      for (i=0;i<a_lists[l].dim();i++) {
	for (j=0;j<a_lists[l].dim();j++) {
	  for (k=0;k<a_lists[l+1].dim();k++) {
	    for (m=0;m<a_lists[l+1].dim();m++) {
	      for (n=0;n<a_lists[l+2].dim();n++) {
		for (o=0;o<a_lists[l+2].dim();o++) {


		  place=
		    k*
		    a_lists[l+1].dim()*
		    a_lists[l+2].dim()*
		    a_lists[l+2].dim()+
		    m*
		    a_lists[l+2].dim()*
		    a_lists[l+2].dim()+
		    n*
		    a_lists[l+2].dim()+
		    o;


		  if (b[g][l+1][place]>0)
		    beta[g][l][
			       i*
			       a_lists[l  ].dim()*
			       a_lists[l+1].dim()*
			       a_lists[l+1].dim()+
			       j*
			       a_lists[l+1].dim()*
			       a_lists[l+1].dim()+
			       k*
			       a_lists[l+1].dim()+
			       m
		    ]+=
		      beta[g][l+1][place]*
		      trans[l][
			       i*a_lists[l+1].dim()*a_lists[l+2].dim()+
			       k*a_lists[l+2].dim()+
			       n]*
		      trans[l][
			       j*a_lists[l+1].dim()*a_lists[l+2].dim()+
			       m*a_lists[l+2].dim()+
			       o]*
		      b[g][l+1][place];
		}
	      }
	    }
	  }
	}
      }
    }
  }

  if (debug) {
    g=0;
    cout << "beta [for genotype " << g << "]" << endl;
    for (l=0;l<numloc-1;l++) {
      for (i=0;i<beta[g][l].dim();i++) {
	cout << beta[g][l][i] << " "; }
      cout << endl;
    }
  }

  //check that all entries are between 0 and 1
  for (g=0;g<numgen;g++) {
    for (l=0;l<numloc-1;l++) {
      for (i=0;i<beta[g][l].dim();i++) {
	if ( !(0.0<=beta[g][l][i] && beta[g][l][i]<1.0001) ) {
	    errorfs << "ERROR (hap est): beta[" << g << "][" << l << "][" << i
		    << "] is not between 0 and 1 ["
		    << beta[g][l][i] << "]" << endl;
	    exit(1); } } } }

}

void makeximkv2_Gfr(dvectorss &xi, const ivectors &a_lists,
		    const dvectorss &alpha, const dvectorss &beta, 
		    const dvectors &trans, const dvectorss &b, 
		    const int &debug)
{
  //makes xi object assuming genotype controls

  int g,i,j,k,l,m,n,o,placel,placelp1;
  double denom;
  const int numgen=alpha.dim();
  const int numloc=trans.dim()+2;

  //clear xi
  for (g=0;g<numgen;g++) {
    for (l=0;l<numloc-2;l++) {
      xi[g][l].clear(0.0); } }

  //make xi
  for (g=0;g<numgen;g++) {
    
    for (l=0;l<numloc-2;l++) {  
      
      for (i=0;i<a_lists[l].dim();i++) {
	for (j=0;j<a_lists[l].dim();j++) {
	  for (k=0;k<a_lists[l+1].dim();k++) {
	    for (m=0;m<a_lists[l+1].dim();m++) {
	      for (n=0;n<a_lists[l+2].dim();n++) {
		for (o=0;o<a_lists[l+2].dim();o++) {
		  
		  placel=
		    i*
		    a_lists[l  ].dim()*
		    a_lists[l+1].dim()*
		    a_lists[l+1].dim()+
		    j*
		    a_lists[l+1].dim()*
		    a_lists[l+1].dim()+
		    k*
		    a_lists[l+1].dim()+
		    m;
		  
		  
		  placelp1=
		    k*
		    a_lists[l+1].dim()*
		    a_lists[l+2].dim()*
		    a_lists[l+2].dim()+
		    m*
		    a_lists[l+2].dim()*
		    a_lists[l+2].dim()+
		    n*
		    a_lists[l+2].dim()+
		    o;
		  
		  
		  if (b[g][l+1][placelp1]>0)
		    xi[g][l][
			     i*a_lists[l].dim()*
			     a_lists[l+1].dim()*a_lists[l+1].dim()*
			     a_lists[l+2].dim()*a_lists[l+2].dim()+
			     j*
			     a_lists[l+1].dim()*a_lists[l+1].dim()*
			     a_lists[l+2].dim()*a_lists[l+2].dim()+
			     k*
			     a_lists[l+1].dim()*
			     a_lists[l+2].dim()*a_lists[l+2].dim()+
			     m*
			     a_lists[l+2].dim()*a_lists[l+2].dim()+
			     n*
			     a_lists[l+2].dim()+
			     o
		    ]=
		      alpha[g][l][placel]*
		      trans[l][
			       i*a_lists[l+1].dim()*a_lists[l+2].dim()+
			       k*a_lists[l+2].dim()+
			       n]*
		      trans[l][
			       j*a_lists[l+1].dim()*a_lists[l+2].dim()+
			       m*a_lists[l+2].dim()+
			       o]*
		      b[g][l+1][placelp1]*
		      beta[g][l+1][placelp1];

		  //delete me
		  //if (g==75 && l==0 && 
		  //  i==0 && j==0 && k==0 && m==0 && n==0 && o==0)
		  //cout << "te " << alpha[g][l][placel] << " "
		  // << beta[g][l+1][placelp1] << " " 
		  //<< b[g][l+1][placelp1] << endl;
		}
	      }
	    }
	  }
	}
      }
      
    }
    //make denominator
    denom=0.0;
    for (i=0;i<xi[g][0].dim();i++) denom+=xi[g][0][i];
      if (denom<=0) {
	errorfs << "ERROR: in makeximkv2_Gfr,denominator for xi at genotype " 
	     << g << "(+1) <=0 [" << denom << "] " << endl;
	exit(1); }
      //divide all entries by this number
      for (l=0;l<numloc-2;l++) {  
	for (i=0;i<xi[g][l].dim();i++) { 
	  xi[g][l][i]/=denom;
	}
      }
      
  }


  if (debug) {
    g=109;
    cout << "xi [in makeximkv2_Gfr] for genotype " << g << endl;
    for (l=0;l<1;l++) {//delete me
      cout << xi[g][l].dim() << endl;
      for (i=0;i<xi[g][l].dim();i++) {	
	cout << i << " " << xi[g][l][i] << endl; }
      cout << endl;
    }
  }
  

  //check that all entries sum to 1 at each locus for each haplotype
  for (g=0;g<numgen;g++) {
    for (l=0;l<numloc-2;l++) {  
      denom=0.0;
      for (i=0;i<xi[g][l].dim();i++) denom+=xi[g][l][i];
      if ( !(0.9999<denom && denom<1.0001) ) {
	errorfs << "ERROR: xi[" << g << "][" << l 
		<< "][] does not sum to 1.0 ["
		<< denom << "]" << endl;
      }
    }
  }

  
}

void makeh3countsGfrxi(dvectors &h3_freq_lists, 
		       const dvectorss &xi, 
		       const ivectors &a_lists, const int &debug)
{
  //makes haplotype frequencies from xi, assuming haplotype controls
  
  int g,i,j,k,l,m,n,o;
  double value;
  
  const int numgen=xi.dim();
  const int numloc=h3_freq_lists.dim()+2;
  
  //clear h3_freq_lists
  for (l=0;l<numloc-2;l++) {
    h3_freq_lists[l].clear(0.0); } 
  
  for (g=0;g<numgen;g++) {
    for (l=0;l<numloc-2;l++) {
      
      for (i=0;i<a_lists[l  ].dim();i++) {
	for (k=0;k<a_lists[l+1].dim();k++) {
	  for (n=0;n<a_lists[l+2].dim();n++) {
	    for (j=0;j<a_lists[l  ].dim();j++) {
	      for (m=0;m<a_lists[l+1].dim();m++) {
		for (o=0;o<a_lists[l+2].dim();o++) {
		  
		  value=xi[g][l][
				 i*
				 a_lists[l].dim()*
				 a_lists[l+1].dim()*a_lists[l+1].dim()*
				 a_lists[l+2].dim()*a_lists[l+2].dim()+
				 j*
				 a_lists[l+1].dim()*a_lists[l+1].dim()*
				 a_lists[l+2].dim()*a_lists[l+2].dim()+
				 k*
				 a_lists[l+1].dim()*
				 a_lists[l+2].dim()*a_lists[l+2].dim()+
				 m*
				 a_lists[l+2].dim()*a_lists[l+2].dim()+
				 n*
				 a_lists[l+2].dim()+
				 o
		  ];
		  
		  h3_freq_lists[l][
				   i*a_lists[l+1].dim()*a_lists[l+2].dim()+
				   k*a_lists[l+2].dim()+n]+=value;
		  h3_freq_lists[l][
				   j*a_lists[l+1].dim()*a_lists[l+2].dim()+
				   m*a_lists[l+2].dim()+o]+=value;
		  
		}
	      }
	    }
	  }
	}
      }
    }
  }

  //check that frequencies sum to 2*numgen
  for (l=0;l<numloc-2;l++) {
    value=0.0;
    for (i=0;i<h3_freq_lists[l].dim();i++) value+=h3_freq_lists[l][i];
    if ( !(2.0*numgen-0.001<value && value<2.0*numgen+0.001) ) {
      errorfs << "ERROR in makeh3freqlistsHfrxi: h3_freq_lists[" 
	      << l << "] does not sum to " << 2*numgen << " ["
	      << value << "]" << endl;
      exit(1); } }

  if (debug) {
    cout << "h3counts " << endl;
    for (i=0;i<1;i++) { 
      for (j=0;j<h3_freq_lists[i].dim();j++)
	cout << h3_freq_lists[i][j] << " ";
      cout << endl; }
  }
}

void makeh2countsGmkv2(dvectors &h2_freq_lists, const int &numgen,
			 const dvectors &h3_freq_lists, 
			 const ivectors &a_lists, const int &debug)
{
  int i,j,k,l;
  double value=0.0;
  const int numloc=a_lists.dim();

  for (l=0;l<numloc-1;l++) {
    h2_freq_lists[l].clear(0.0);
  }

  for (l=0;l<numloc-2;l++) {

    for (i=0;i<a_lists[l  ].dim();i++) {
      for (j=0;j<a_lists[l+1].dim();j++) {
	for (k=0;k<a_lists[l+2].dim();k++) {

	  if (l==0)
	    h2_freq_lists[l][ i*a_lists[l+1].dim() +j]+=
	      h3_freq_lists[l][
			       i*a_lists[l+1].dim()*a_lists[l+2].dim()+
			       j*a_lists[l+2].dim()+
			       k ];
	  
	  h2_freq_lists[l+1][ j*a_lists[l+2].dim() +k]+=
	    h3_freq_lists[l][
			     i*a_lists[l+1].dim()*a_lists[l+2].dim()+
			     j*a_lists[l+2].dim()+
			     k ];

	}
      }
    }
  }
    
  if (debug) {
    cout << "h2counts " << endl;
    for (l=0;l<2;l++) {
      for (i=0;i<h2_freq_lists[l].dim();i++) {	
	cout << h2_freq_lists[l][i] << " "; }
      cout << endl;
    }
  }

  //check that frequencies sum to 2*numgen for each pair
  for (l=0;l<numloc-1;l++) {
    value=0.0;
    for (i=0;i<h2_freq_lists[l].dim();i++) value+=h2_freq_lists[l][i];
    if ( !(2.0*numgen-0.001<value && value<2.0*numgen+0.001) ) {
      errorfs << "ERROR: h2_freq_lists[" << l << "] does not sum to "
	      << 2*numgen << " ["
	      << value << "]" << endl;
      exit(1); } }
}

void makeh3freqlistsGmkv2(dvectors &h3_freq_lists, dvectors &h2_freq_lists,
			  const dvector &init,
			  const dvectors &trans,const ivectors &a_lists,
			  const int &debug)
{
  int l,i,j,k;
  const int numloc=a_lists.dim();
  double value=0.0;
  
  //delete me
  //cout << init(0)+init(1)+init(2)+init(3) << endl;
  //cout << trans[0][0]+trans[0][1] << " "
  //   << trans[0][2]+trans[0][3] << " "
  //   << trans[0][4]+trans[0][5] << " "
  //   << trans[0][6]+trans[0][7] << endl;

  for (l=0;l<numloc-2;l++) {
    h2_freq_lists[l].clear(0.0);
    h3_freq_lists[l].clear(0.0); }
  h2_freq_lists[numloc-2].clear(0.0);

  h2_freq_lists[0]=init;
  //cout << "init" << " ";
  //for (i=0;i<h2_freq_lists[0].dim();i++) 
  //cout << init(i) << " "; cout << endl;

  for (l=0;l<numloc-2;l++) {

    for (i=0;i<a_lists[l  ].dim();i++) {
      for (j=0;j<a_lists[l+1].dim();j++) {
	for (k=0;k<a_lists[l+2].dim();k++) {

	  h3_freq_lists[l][
			   i*a_lists[l+1].dim()*a_lists[l+2].dim()+
			   j*a_lists[l+2].dim()+
			   k ]=
	    h2_freq_lists[l][ i*a_lists[l+1].dim() +j]*
	    trans[l][
		     i*a_lists[l+1].dim()*a_lists[l+2].dim()+
		     j*a_lists[l+2].dim()+
		     k];

	  //update h2_freq_lists[l+1]
	  h2_freq_lists[l+1][ j*a_lists[l+2].dim() +k]+=
	    h3_freq_lists[l][
			     i*a_lists[l+1].dim()*a_lists[l+2].dim()+
			     j*a_lists[l+2].dim()+
			     k ];
	}
      }
    }
  }
  if (debug) {
    cout << "h3_freq_lists [in makeh3freqlistsGmkv2] " << endl;
    for (l=0;l<numloc-2;l++) {
      for (i=0;i<h3_freq_lists[l].dim();i++) {	
	cout << h3_freq_lists[l][i] << " "; }
      cout << endl;
    }
  }

  //check that frequencies sum to 1 for each pair
  for (l=0;l<numloc-2;l++) {
    value=0.0;
    for (i=0;i<h3_freq_lists[l].dim();i++) value+=h3_freq_lists[l][i];
    if ( !(0.9999<value && value<1.0001) ) {
      errorfs << "ERROR in makeh3freqlistsGmkv2: h3_freq_lists[" 
	      << l << "] does not sum to 1 ["
	      << value << "]" << endl;
      exit(1); } }


}

//bayesian hap freq est. for mkv1

void makehap2fr_hmm_HconBAYES(dvectors &h2_freq_lists, const int &bayes,
			      const ivectors &a_lists, 
			      const igrid &controls, 
			      const int &debug)
{

  //estimates two-marker haplotype frequencies, using HMM, assuming
  //haplotype controls

  int l,h,i,j,iter,conv;
  double diff;
  const int numloc=controls.dimy();
  const int numhap=controls.dimx();
  double sum=0.0;

  //make vector indicating haplotypes with missing values
  int is_missing=0;
  ivector missing; missing.init(numhap);
  make_missing_Hfr(missing, is_missing, controls, 0);

  //delete me
  for (i=0;i<missing.dim();i++) missing(i)=0;  

  //declare and intialize and clear lists of allele freqs 
  //corresponding to a_lists
  dvectors a_freq_lists; 
  a_freq_lists.init(numloc);
  for (l=0;l<numloc;l++) 
    { a_freq_lists[l].init(a_lists[l].dim()); a_freq_lists[l].clear(0); }
  
  //declare and initialize and clear copy to be used to determine convergence
  //dvectors h2_freq_lists_copy(h2_freq_lists.dim());
  //for (l=0;l<numloc-1;l++) 
  //{ 
  //  h2_freq_lists_copy[l].init(h2_freq_lists[l].dim()); 
  //  h2_freq_lists_copy[l].clear(0); 
  //}
  
  //declare and initialize and clear hmm objects
  dvectorss b(numhap);
  dvectorss alpha(numhap),beta(numhap),gamma(numhap),xi(numhap);
  for (h=0;h<numhap;h++) {
    if (missing[h]==0) {
      b[h].init(numloc);
      alpha[h].init(numloc);
      beta[h].init(numloc);
      gamma[h].init(numloc);
      for (l=0;l<numloc;l++) {
	b[h][l].init(a_lists[l].dim()); b[h][l].clear(0);
	alpha[h][l].init(a_lists[l].dim()); alpha[h][l].clear(0);
	beta[h][l].init(a_lists[l].dim()); beta[h][l].clear(0);
	gamma[h][l].init(a_lists[l].dim()); gamma[h][l].clear(0);
      }
    }
  }
  //initialize xi vectors for all haps
  for (h=0;h<numhap;h++) {
    xi[h].init(numloc-1);
    for (l=0;l<numloc-1;l++) {
      xi[h][l].init(a_lists[l].dim()*a_lists[l+1].dim());
      xi[h][l].clear(0); }
  }
  
  //make init  
  dvector init;
  init.init( a_lists[0].dim() ); 
  init.clear(0.0);

  //make copy to determine convergence
  dvector initC; initC.init( init.dim() ); initC.clear(0.0);

  //trans
  dvectors trans(numloc-1);
  for (l=0;l<numloc-1;l++) {
    trans[l].init(a_lists[l].dim()*a_lists[l+1].dim());
    trans[l].clear(0); }

  //copy to determine convergence
  dvectors transC(numloc-1);
  for (l=0;l<numloc-1;l++) {
    transC[l].init(trans[l].dim()); transC[l].clear(0); }

  //assign initial values to parameters
  for (i=0;i<init.dim();i++) init(i)=1.0/init.dim();
  for (l=0;l<numloc-1;l++) {
    for (i=0;i<a_lists[l].dim();i++) {
      for (j=0;j<a_lists[l+1].dim();j++) {
	trans[l][i*a_lists[l+1].dim()+j]=
	  1.0/( a_lists[l+1].dim() ); } } }
 
  //fill in xi for haplotypes with no missing alleles
  makexi_nomissing_Hfr(xi, missing, a_lists, controls, 0);

  //make observation probabilities [b]
  makeb_Hfr(b, missing, a_lists, controls, 0);

  //make initial allele freq's
  //make_a_freq_lists_Gfr(a_freq_lists, h2_freq_lists, 0);

  //iterate hmm
  int maxiter=200;
  for (iter=1;iter<=maxiter;iter++) {

    //copy init to use later in checking convergence
    initC=init;

    //copy hapl. freq. to use later in checking convergence
    for (l=0;l<transC.dim();l++) 
      transC[l]=trans[l];

    //copy hapl. freq. to use later in checking convergence
    //for (l=0;l<h2_freq_lists.dim();l++) h2_freq_lists_copy[l]=h2_freq_lists[l];

    //make trans
    //maketrans_Hfr(trans, a_freq_lists, h2_freq_lists, 0);

    //make alpha
    makealpha_HfrBAYES(alpha, missing, a_lists, init, trans, b, 0);

    //make beta
    makebeta_HfrBAYES(beta, missing, a_lists, trans, b, 0);

    //make xi
    makexi_Hfr(xi, missing, alpha, beta, trans, b, 0);

    //update h2_freq_lists
    //make_h2_freq_lists_Hfr(h2_freq_lists, xi, 0);
    
    //make allele freq's
    //make_a_freq_lists_Gfr(a_freq_lists, h2_freq_lists, 0);

    //give expected counts of 2-locus haplotypes
    makeh2countsHfrxi(h2_freq_lists, xi, a_lists, 0);
    make_a_freq_lists_HfrBAYES(a_freq_lists, numhap, h2_freq_lists, 0);

    //Bayesian adjustment/update to init
    init=a_freq_lists[0];
    for (i=0;i<init.dim();i++) { 
      init(i)+=1.0*bayes; 
      init(i)/=numhap+1.0*init.dim()*bayes; }

    //check that init sums to 1
    sum=0;
    for (i=0;i<init.dim();i++) sum+=init(i);
    if ( !(0.9999<sum && sum<1.0001) ) {
      errorfs << "ERROR: Bayesian adjustment to init in hap. freq. est.; "
	      << "sums to " << sum << endl;
      exit(1); }

    //Bayesian adjustment/update to trans
    for (l=0;l<numloc-1;l++) {
      for (i=0;i<a_lists[l  ].dim();i++) {
	for (j=0;j<a_lists[l+1].dim();j++) {
	  if (h2_freq_lists[l][i*a_lists[l+1].dim()+j]+1.0*bayes <1e-100
	      && a_freq_lists[l][i]+ 1.0*a_lists[l+1].dim()*bayes<1e-100)
	    trans[l][i*a_lists[l+1].dim()+j]=0.0;
	  else
	    trans[l][i*a_lists[l+1].dim()+j]=
	      ( h2_freq_lists[l][i*a_lists[l+1].dim()+j]+1.0*bayes )/
	      ( a_freq_lists[l][i]+ 1.0*a_lists[l+1].dim()*bayes);
	}
	//if they're all 0, set first to 1
	conv=0;
	for (j=0;j<a_lists[l+1].dim();j++) 
	  if (trans[l][i*a_lists[l+1].dim()+j]>0) { conv=1; break; }
	if (conv==0) trans[l][i*a_lists[l+1].dim()  ]=1.0;
      }
    }//l

    //check trans
    for (l=0;l<numloc-1;l++) {
      for (i=0;i<a_lists[l  ].dim();i++) {
	sum=0;
	for (j=0;j<a_lists[l+1].dim();j++) {
	  sum+=trans[l][i*a_lists[l+1].dim()+j]; }
	if ( !(0.9999<sum && sum<1.0001) ) {
	  errorfs << "ERROR in Bayes adj. of trans. Doesn't sum to 1 for "
		  << "i " << i << " j " << j 
		  << " at locus " << l << " ( " << sum << " ) " << endl;
	  exit(1); }} }


    //check convergence
    conv=0;
    for (i=0;i<init.dim();i++) {
      diff=init(i)-initC(i);
      if ( diff<-1e-8 || 1e-8<diff ) { conv=1; break; } }
    if (conv==0) {
      for (l=0;l<numloc-1;l++) {
	for (i=0;i<trans[l].dim();i++) {
	  diff=trans[l][i]-transC[l][i];
	  if ( diff<-1e-8 || 1e-8<diff ) { conv=1; break; } }
	if (conv==1) break; } }
    if (conv==0) break; 

  }//iter

  //check that two ways of getting a_freq_lists are identical
  //check_a_freq_lists(a_freq_lists, h2_freq_lists);

  //make h2_freq_lists from trans, init
  makeh2freqlistsH(h2_freq_lists,a_freq_lists,init,trans,a_lists,0);

  if (debug) {
    cout << "iter " << iter << endl;
    cout << "init " << endl;
    for (i=0;i<init.dim();i++) cout << init(i) << " "; cout << endl;
    cout << "trans " << endl;
    for (l=0;l<numloc-1;l++) {
      for (i=0;i<trans[l].dim();i++) {
	cout << trans[l][i] << " "; }
      cout << endl; }

    cout << "h2_freq_lists" << endl;
    for (l=0;l<numloc-1;l++) {
      for (i=0;i<h2_freq_lists[l].dim();i++) {	
	cout << h2_freq_lists[l][i] << " "; }
      cout << endl;
    }

    //make likelihood
    //make alpha
    makealpha_HfrBAYES(alpha, missing, a_lists, init, trans, b, 0);
    double lik=0,loglik=0;
    for (h=0;h<numhap;h++) {
      lik=0;
      for (i=0;i<alpha[h][numloc-1].dim();i++) {
	lik+=alpha[h][numloc-1][i]; }
      if ( !(0<=lik && lik<=1.0) ) {
	errorfs << "ERROR: lik. for hap " << h << "+1 is " << lik << endl;
	exit(1); }
      loglik+=log(lik); }
    cout << "lik " << loglik << endl; 

  }
}


void makealpha_HfrBAYES(dvectorss &alpha, const ivector &missing,
			const ivectors &a_lists,
			const dvector &init,
			const dvectors &trans, const dvectorss &b,
			const int &debug)
{
  //makes alpha object, assuming haplotype controls

  int h,i,j,l;
  const int numhap=alpha.dim();
  const int numloc=a_lists.dim();

  //clear alpha
  for (h=0;h<numhap;h++) {
    if (missing[h]==0) {
      for (l=0;l<numloc;l++) {
	alpha[h][l].clear(0); } } }

  for (h=0;h<numhap;h++) {

    if (missing[h]==0) {
      //make first column
      for (i=0;i<alpha[h][0].dim();i++) {

	if (b[h][0][i] != 0) 
	  alpha[h][0][i]=
	    init(i)*b[h][0][i];
	
      }
    
    
      //make rest of alpha
      for (l=0;l<numloc-1;l++) {
      
	for (j=0;j<a_lists[l+1].dim();j++) {
	  for (i=0;i<a_lists[l  ].dim();i++) {

	    alpha[h][l+1][j]+=
	      alpha[h][l  ][i]*
	      trans[l][i*a_lists[l+1].dim()+j]*
	      b[h][l+1][j];

	    
	  }
	}
      }
    }
  }
  
  if (debug) {
    cout << "alpha [in maketrans_HfrBAYES] " << endl;
    for (h=0;h<numhap;h++) {
      if (missing[h]==0) {
	for (l=0;l<numloc;l++) {
	  for (i=0;i<alpha[h][l].dim();i++) {	
	    cout << alpha[h][l][i] << " "; }
	  cout << endl;
	}
      }
    }
  }
  //check that all entries are between 0 and 1
  for (h=0;h<numhap;h++) {
    if (missing[h]==0) 
      for (l=0;l<numloc;l++) {
	for (i=0;i<alpha[h][l].dim();i++) {
	  if ( !(0.0<=alpha[h][l][i] && alpha[h][l][i]<1.0001) ) {
	    errorfs << "ERROR (hap est): alpha[" << h << "][" << l << "][" << i
		 << "] is not between 0 and 1 ["
		    << alpha[h][l][i] << "]" << endl;
	    exit(1); } } } }
  
  
}

void makebeta_HfrBAYES(dvectorss &beta, const ivector &missing, 
		       const ivectors &a_lists,
		       const dvectors &trans, const dvectorss &b,
		       const int &debug)
{
  //makes beta object, assuming haplotype controls

  int h,i,j,l;
  const int numhap=beta.dim();
  const int numloc=trans.dim()+1;

  //clear beta
  for (h=0;h<numhap;h++) {
    if (missing[h]==0) {
      for (l=0;l<numloc;l++) {
	beta[h][l].clear(0); } } }

  for (h=0;h<numhap;h++) {

    if (missing[h]==0) {
      //make last column
      for (i=0;i<beta[h][numloc-1].dim();i++) {
	
	beta[h][numloc-1][i]=1.0;
	
      }
      
      //make rest of beta
      for (l=numloc-2;l>=0;l--) {
	
	for (i=0;i<a_lists[l  ].dim();i++) {
	  for (j=0;j<a_lists[l+1].dim();j++) {
	  
	    beta[h][l][i]+=
	      beta[h][l+1][j]*
	      trans[l][i*a_lists[l+1].dim()+j]*
	      b[h][l+1][j];

	    //delete me
	    //if (h==0 && l==0)
	    //cout << "beta " << i << " " << j << " "
	    //   << beta[h][l+1][j] << " "
	    //   << trans[l][i*a_lists[l+1].dim()+j] << " "
	    //   << b[h][l+1][j] << " "
	    //   << endl;

	  }
	}
      }
    }
  }

  if (debug) {
    cout << "beta [in makebeta_HfrBAYES] " << endl;
    for (h=0;h<numhap;h++) {
      if (missing[h]==0) {
	for (l=0;l<numloc;l++) {
	  for (i=0;i<beta[h][l].dim();i++) {	
	    cout << beta[h][l][i] << " "; }
	  cout << endl;
	}
      }
    }
  }

  //check that all entries are between 0 and 1
  for (h=0;h<numhap;h++) {
    if (missing[h]==0)
      for (l=0;l<numloc;l++) {
	for (i=0;i<beta[h][l].dim();i++) {
	  if ( !(0.0<=beta[h][l][i] && beta[h][l][i]<1.0001) ) {
	    errorfs << "ERROR (hap est): beta[" << h << "][" << l << "][" << i
		    << "] is not between 0 and 1 ["
		    << beta[h][l][i] << "]" << endl;
	    exit(1); } } } }
  
  
}

void makeh2countsHfrxi(dvectors &h2_freq_lists, 
		       const dvectorss &xi, 
		       const ivectors &a_lists, const int &debug)
{
  //makes haplotype frequencies from xi, assuming haplotype controls

  int h,i,j,l;
  double value;

  const int numhap=xi.dim();
  const int numloc=h2_freq_lists.dim()+1;

  //clear h2_freq_lists
  for (l=0;l<numloc-1;l++) {
    h2_freq_lists[l].clear(0); } 

  for (h=0;h<numhap;h++) {
    for (l=0;l<numloc-1;l++) {
      
      for (i=0;i<a_lists[l  ].dim();i++) {
	for (j=0;j<a_lists[l+1].dim();j++) {

	  h2_freq_lists[l][i*a_lists[l+1].dim()+j]
	    +=xi[h][l][i*a_lists[l+1].dim()+j];
	}
      }
    }
  }
  
  //divide by numhap
  //for (l=0;l<numloc-2;l++) {
  //for (i=0;i<h3_freq_lists[l].dim();i++) {
  //h3_freq_lists[l][i]/=numhap; } }
  
  if (debug) {
    cout << "h2_freq_lists [in makexi_Hfr] " << endl;
    for (l=0;l<numloc-1;l++) {
      for (i=0;i<h2_freq_lists[l].dim();i++) {	
	cout << h2_freq_lists[l][i] << " "; }
      cout << endl;
    }
  }

  //check that frequencies sum to numhap for each pair
  for (l=0;l<numloc-1;l++) {
    value=0.0;
    for (i=0;i<h2_freq_lists[l].dim();i++) value+=h2_freq_lists[l][i];
    if ( !(numhap-0.0001<value && value<numhap+0.0001) ) {
      errorfs << "ERROR in makeh2countsHfrxi: h2_freq_lists[" 
	      << l << "] does not sum to numhap ["
	      << value << "]" << endl;
      exit(1); } }

  
}

void make_a_freq_lists_HfrBAYES(dvectors &a_freq_lists, const int &numhap, 
				const dvectors &h2_freq_lists,
				const int &debug)
{
  //this function generates the list of allele frequencies at the loci
  //by summing over the joint frequencies

  int i,j,l;
  const int numloc=a_freq_lists.dim();
  //double diff;

  //check that get same estimates by summing first and second elements of pair
  //dvectors scratch; scratch.init(a_freq_lists.dim());
  //for (l=0;l<numloc;l++) {
  //scratch[l].init(a_freq_lists[l].dim());
  // scratch[l].clear(0); }

  //zero-out lists
  for (l=0;l<numloc;l++) {
    a_freq_lists[l].clear(0); 
  }

  for (l=0;l<numloc-1;l++) {
    for (i=0;i<a_freq_lists[l].dim();i++) {
      for (j=0;j<a_freq_lists[l+1].dim();j++) {

	a_freq_lists[l][i]+=h2_freq_lists[l][i*a_freq_lists[l+1].dim()+j];
	if(l==numloc-2) 
	  a_freq_lists[l+1][j]+=h2_freq_lists[l][i*a_freq_lists[l+1].dim()+j];

	//scratch[l+1][j]+=h2_freq_lists[l][i*a_freq_lists[l+1].dim()+j];
      }
    }
  }

  //check that columns sum to numhap
  double value;
  for (l=0;l<numloc;l++) {
    value=0.0;
    for (i=0;i<a_freq_lists[l].dim();i++) {
      value+=a_freq_lists[l][i];
    }
    if ( !(numhap-0.0001<value && value <numhap+0.0001) ) {
      errorfs << "ERROR: all. f.'s of makehap2fr_hmm_Gcon do not sum to "
	      << "numhap at locus "
	      << l << "(+1) [" << value << "] " 
	      << endl;
      //assert( 0.9999<value && value <1.0001 );
      exit(1); }
  }
  
  
  //check that scratch==a_freq_lists
  //for (l=1;l<numloc;l++) {
  //for (i=0;i<a_freq_lists[l].dim();i++) {
  //  diff=scratch[l][i]-a_freq_lists[l][i];
  //  if ( !(-0.0001<diff && diff<0.0001) ) {
  //errorfs << "ERROR in make_a_freq_lists_Gfr; " << endl
  //     << "two estimates for a_freq_lists_Gfr[" << l 
  //     << "][" << i << "] are not equal" << endl;
  //exit(1);
  //  }
  //}
  //}

  if (debug) {
    cout << "a_freq_lists" << endl;
    for (l=0;l<numloc;l++) {
      for (i=0;i<a_freq_lists[l].dim();i++) cout << a_freq_lists[l][i] << " "; 
      cout << endl; 
    }
  }
  
}

void makeh2freqlistsH(dvectors &h2_freq_lists,
			 dvectors &a_freq_lists,
			 const dvector &init,
			 const dvectors &trans,const ivectors &a_lists,
			 const int &debug)
{
  int l,i,j;
  const int numloc=a_lists.dim();
  double value=0.0;
  
  for (l=0;l<numloc-1;l++) {
    a_freq_lists[l].clear(0.0);
    h2_freq_lists[l].clear(0.0); }
  a_freq_lists[numloc-1].clear(0.0);
  
  a_freq_lists[0]=init;
  //cout << "init" << " ";
  //for (i=0;i<h2_freq_lists[0].dim();i++) 
  //cout << init(i) << " "; cout << endl;

  for (l=0;l<numloc-1;l++) {

    for (i=0;i<a_lists[l  ].dim();i++) {
      for (j=0;j<a_lists[l+1].dim();j++) {

	h2_freq_lists[l][i*a_lists[l+1].dim()+j]=
	  a_freq_lists[l][i]*
	  trans[l][i*a_lists[l+1].dim() +j];

	//update a_freq_lists[l+1]
	a_freq_lists[l+1][j]+=
	  h2_freq_lists[l][i*a_lists[l+1].dim()+j];
      }
    }
  }
  if (debug) {
    cout << "h2_freq_lists [in makeh2freqlistsH] " << endl;
    for (l=0;l<numloc-1;l++) {
      for (i=0;i<h2_freq_lists[l].dim();i++) {	
	cout << h2_freq_lists[l][i] << " "; }
      cout << endl;
    }
  }

  //check that frequencies sum to 1 for each pair
  for (l=0;l<numloc-1;l++) {
    value=0.0;
    for (i=0;i<h2_freq_lists[l].dim();i++) value+=h2_freq_lists[l][i];
    if ( !(0.9999<value && value<1.0001) ) {
      errorfs << "ERROR in makeh2freqlistsH: h2_freq_lists[" 
	      << l << "] does not sum to 1 ["
	      << value << "]" << endl;
      exit(1); } }

  
}