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

#include <fstream.h>
extern ofstream errorfs;

//declarations
#include "declare.h"

/*
  This file contains source code for functions involved in estimation
  of 1/tau and p given the location of the gene and the ancestral
  hap., assuming genotypes from affecteds-- much of code parallels that
  in estfunctions.C
 */

void makenullsharegen(dgrid &nullprobs, const igrid &data,
		      const igrid &freqallele, const dgrid &freqfreq,
		      const dgrid &condfreq2l, const dgrid &condfreq2r,
		      const dgrid &condfreq2C,
		      const ivector &allelecounts,
		      const int &LE, const int &RE, const int &C)
{
  //this function computes parts of observation probabilities that do 
  //not involve the parameters-- it stores these pieces in nullprobs
  
  int i,j;
  const int numgen=data.dimx()/2;
  double f0,f1,cf00,cf10,cf01,cf11;
  int LHSle, LHSre, RHSle, RHSre, Cadj=-1;

  if (LE<C && C<RE) {
    LHSle=LE+1; LHSre=C-1;
    RHSle=C+1;  RHSre=RE-1; }
  else {
    if (C<LE) {
      Cadj=LE;
      LHSle=0,LHSre=-1;
      RHSle=LE;  RHSre=RE-1; }
    else {//RE<C
      Cadj=RE;
      LHSle=LE+1; LHSre=RE;
      RHSle=0,RHSre=-1; }
  }

  for (i=0;i<numgen;i++) {
    
    if (LE<C)
      {
	//make first column
	if (data(2*i  ,LE)==0) {
	  if (data(2*i+1,LE)==0) {
	    f0=1;
	    f1=1; }
	  else {
	    f0=1;
	    f1=getfreq(freqallele, freqfreq, data(2*i+1,LE), LE); } }
	else {
	  if (data(2*i+1,LE)==0) {
	    f0=getfreq(freqallele, freqfreq, data(2*i  ,LE), LE);
	    f1=1; }
	  else {
	    f0=getfreq(freqallele, freqfreq, data(2*i  ,LE), LE);
	    f1=getfreq(freqallele, freqfreq, data(2*i+1,LE), LE); } }
	
	nullprobs(7*i  ,LE)=1;    //AA
	nullprobs(7*i+1,LE)=f1;   //AN1
	nullprobs(7*i+2,LE)=f0;   //NA1
	nullprobs(7*i+3,LE)=f0*f1;//NN1
	nullprobs(7*i+4,LE)=0;    //AN0
	nullprobs(7*i+5,LE)=0;    //NA0
	nullprobs(7*i+6,LE)=0;    //NN0
      }


    //make columns LE+1:C-1
    for (j=LHSle;j<=LHSre;j++) {
	  
      if (data(2*i  ,j)==0) {
	if (data(2*i+1,j)==0) {
	  
	  cf00=1;
	  cf10=1;
	  cf01=1;
	  cf11=1; }
	
	else {//data(2*i  ,j)==0),data(2*i+1,j)!=0)
	  
	  if (data(2*i  ,j-1)==0) {
	    if (data(2*i+1,j-1)==0) {
	      cf00=1;
	      cf10=1;
	      cf01=getfreq(freqallele, freqfreq, data(2*i+1,j), j);
	      cf11=getfreq(freqallele, freqfreq, data(2*i+1,j), j); 
	    } 
	    else {
	      cf00=1;
	      cf10=1;
	      cf01=getfreq(freqallele, freqfreq, data(2*i+1,j), j);
	      cf11=getcondfreq2l(freqallele, allelecounts, condfreq2l,
				 j-1, data(2*i+1,j-1), data(2*i+1,j)); 
	    } }
	  else {
	    if (data(2*i+1,j-1)==0) {
	      cf00=1;
	      cf10=1;
	      cf01=getcondfreq2l(freqallele, allelecounts, condfreq2l,
				 j-1, data(2*i  ,j-1), data(2*i+1,j));
	      cf11=getfreq(freqallele, freqfreq, data(2*i+1,j), j); 
	    } 
	    else {
	      cf00=1;
	      cf10=1;
	      cf01=getcondfreq2l(freqallele, allelecounts, condfreq2l,
				 j-1, data(2*i  ,j-1), data(2*i+1,j));
	      cf11=getcondfreq2l(freqallele, allelecounts, condfreq2l,
				 j-1, data(2*i+1,j-1), data(2*i+1,j)); 
	    } } 
	  
	}
      }
      
      else {//(data(2*i  ,j)!=0)
	if (data(2*i+1,j)==0) {
	  
	  if (data(2*i  ,j-1)==0) {
	    if (data(2*i+1,j-1)==0) {
	      cf00=getfreq(freqallele, freqfreq, data(2*i  ,j), j);
	      cf10=getfreq(freqallele, freqfreq, data(2*i  ,j), j);
	      cf01=1;
	      cf11=1;
	    } 
	    else {
	      cf00=getfreq(freqallele, freqfreq, data(2*i  ,j), j);;
	      cf10=getcondfreq2l(freqallele, allelecounts, condfreq2l,
				 j-1, data(2*i+1,j-1), data(2*i  ,j));
	      cf01=1;
	      cf11=1;
	    } }
	  else {
	    if (data(2*i+1,j-1)==0) {
	      cf00=getcondfreq2l(freqallele, allelecounts, condfreq2l,
				 j-1, data(2*i  ,j-1), data(2*i  ,j));
	      cf10=getfreq(freqallele, freqfreq, data(2*i  ,j), j);
	      cf01=1;
	      cf11=1;
	    } 
	    else {
	      cf00=getcondfreq2l(freqallele, allelecounts, condfreq2l,
				 j-1, data(2*i  ,j-1), data(2*i  ,j));
	      cf10=getcondfreq2l(freqallele, allelecounts, condfreq2l,
				 j-1, data(2*i+1,j-1), data(2*i  ,j));
	      cf01=1;
	      cf11=1;
	    } 
	  }
	      
	}
	    
	else {//data(2*i  ,j)!=0),data(2*i+1,j)!=0)
	      
	  if (data(2*i  ,j-1)==0) {
	    if (data(2*i+1,j-1)==0) {
	      cf00=getfreq(freqallele, freqfreq, data(2*i  ,j), j);
	      cf10=getfreq(freqallele, freqfreq, data(2*i  ,j), j);
	      cf01=getfreq(freqallele, freqfreq, data(2*i+1,j), j);
	      cf11=getfreq(freqallele, freqfreq, data(2*i+1,j), j); 
	    } 
	    else {
	      cf00=getfreq(freqallele, freqfreq, data(2*i  ,j), j);
	      cf10=getcondfreq2l(freqallele, allelecounts, condfreq2l,
				 j-1, data(2*i+1,j-1), data(2*i  ,j));
	      cf01=getfreq(freqallele, freqfreq, data(2*i+1,j), j);
	      cf11=getcondfreq2l(freqallele, allelecounts, condfreq2l,
				 j-1, data(2*i+1,j-1), data(2*i+1,j)); 
	    } }
	  else {
	    if (data(2*i+1,j-1)==0) {
	      cf00=getcondfreq2l(freqallele, allelecounts, condfreq2l,
				 j-1, data(2*i  ,j-1), data(2*i  ,j));
	      cf10=getfreq(freqallele, freqfreq, data(2*i  ,j), j); 
	      cf01=getcondfreq2l(freqallele, allelecounts, condfreq2l,
				 j-1, data(2*i  ,j-1), data(2*i+1,j));
	      cf11=getfreq(freqallele, freqfreq, data(2*i+1,j), j); 
	    } 
	    else {
	      cf00=getcondfreq2l(freqallele, allelecounts, condfreq2l,
				 j-1, data(2*i  ,j-1), data(2*i  ,j));
	      cf10=getcondfreq2l(freqallele, allelecounts, condfreq2l,
				 j-1, data(2*i+1,j-1), data(2*i  ,j));
	      cf01=getcondfreq2l(freqallele, allelecounts, condfreq2l,
				 j-1, data(2*i  ,j-1), data(2*i+1,j));
	      cf11=getcondfreq2l(freqallele, allelecounts, condfreq2l,
				 j-1, data(2*i+1,j-1), data(2*i+1,j)); 
	    } } 
	  
	}
      }
	  
      nullprobs(7*i  ,j)=1;        //AA
      nullprobs(7*i+1,j)=cf11;     //AN1
      nullprobs(7*i+2,j)=cf00;     //NA1
      nullprobs(7*i+3,j)=cf00*cf11;//NN1
      nullprobs(7*i+4,j)=cf01;     //AN0
      nullprobs(7*i+5,j)=cf10;     //NA0
      nullprobs(7*i+6,j)=cf10*cf01;//NN0
      
    }


    //insert Cstuff here
    if (C<LE || RE<C) {
      
      //make column C
      if (data(2*i  ,C)==0) {
	if (data(2*i+1,C)==0) {
	  
	  cf00=1;
	  cf10=1;
	  cf01=1;
	  cf11=1; }
	
	else {//data(2*i  ,C)==0),data(2*i+1,C)!=0)
	  
	  if (data(2*i  ,Cadj)==0) {
	    if (data(2*i+1,Cadj)==0) {
	      cf00=1;
	      cf10=1;
	      cf01=getfreq(freqallele, freqfreq, data(2*i+1,C), C);
	      cf11=getfreq(freqallele, freqfreq, data(2*i+1,C), C); 
	    } 
	    else {
	      cf00=1;
	      cf10=1;
	      cf01=getfreq(freqallele, freqfreq, data(2*i+1,C), C);
	      cf11=getjointfreq2C(freqallele, condfreq2C, 
				  data(2*i+1,Cadj), Cadj, data(2*i+1,C),C);
	    } }
	  else {
	    if (data(2*i+1,Cadj)==0) {
	      cf00=1;
	      cf10=1;
	      cf01=getjointfreq2C(freqallele, condfreq2C, 
				  data(2*i  ,Cadj), Cadj, data(2*i+1,C),C);
	      cf11=getfreq(freqallele, freqfreq, data(2*i+1,C), C); 
	    } 
	    else {
	      cf00=1;
	      cf10=1;
	      cf01=getjointfreq2C(freqallele, condfreq2C, 
				  data(2*i  ,Cadj), Cadj, data(2*i+1,C),C);
	      cf11=getjointfreq2C(freqallele, condfreq2C, 
				  data(2*i+1,Cadj), Cadj, data(2*i+1,C),C);
	    } } 
	  
	}
      }
      
      else {//(data(2*i  ,C)!=0)
	if (data(2*i+1,C)==0) {
	  
	  if (data(2*i  ,Cadj)==0) {
	    if (data(2*i+1,Cadj)==0) {
	      cf00=getfreq(freqallele, freqfreq, data(2*i  ,C), C);
	      cf10=getfreq(freqallele, freqfreq, data(2*i  ,C), C);
	      cf01=1;
	      cf11=1;
	    } 
	    else {
	      cf00=getfreq(freqallele, freqfreq, data(2*i  ,C), C);;
	      cf10=getjointfreq2C(freqallele, condfreq2C, 
				  data(2*i+1,Cadj), Cadj, data(2*i  ,C),C);
	      cf01=1;
	      cf11=1;
	    } }
	  else {
	    if (data(2*i+1,Cadj)==0) {
	      cf00=getjointfreq2C(freqallele, condfreq2C, 
				  data(2*i  ,Cadj), Cadj, data(2*i  ,C),C);
	      cf10=getfreq(freqallele, freqfreq, data(2*i  ,C), C);
	      cf01=1;
	      cf11=1;
	    } 
	    else {
	      cf00=getjointfreq2C(freqallele, condfreq2C, 
				  data(2*i  ,Cadj), Cadj, data(2*i  ,C),C);
	      cf10=getjointfreq2C(freqallele, condfreq2C, 
				  data(2*i+1,Cadj), Cadj, data(2*i  ,C),C);
	      cf01=1;
	      cf11=1;
	    } 
	  }
	  
	}
	
	else {//data(2*i  ,C)!=0),data(2*i+1,C)!=0)
	  
	  if (data(2*i  ,Cadj)==0) {
	    if (data(2*i+1,Cadj)==0) {
	      cf00=getfreq(freqallele, freqfreq, data(2*i  ,C), C);
	      cf10=getfreq(freqallele, freqfreq, data(2*i  ,C), C);
	      cf01=getfreq(freqallele, freqfreq, data(2*i+1,C), C);
	      cf11=getfreq(freqallele, freqfreq, data(2*i+1,C), C); 
	    } 
	    else {
	      cf00=getfreq(freqallele, freqfreq, data(2*i  ,C), C);
	      cf10=getjointfreq2C(freqallele, condfreq2C, 
				  data(2*i+1,Cadj), Cadj, data(2*i  ,C),C);
	      cf01=getfreq(freqallele, freqfreq, data(2*i+1,C), C);
	      cf11=getjointfreq2C(freqallele, condfreq2C, 
				  data(2*i+1,Cadj), Cadj, data(2*i+1,C),C);
	    } }
	  else {
	    if (data(2*i+1,Cadj)==0) {
	      cf00=getjointfreq2C(freqallele, condfreq2C, 
				  data(2*i  ,Cadj), Cadj, data(2*i  ,C),C);
	      cf10=getfreq(freqallele, freqfreq, data(2*i  ,C), C); 
	      cf01=getjointfreq2C(freqallele, condfreq2C, 
				  data(2*i  ,Cadj), Cadj, data(2*i+1,C),C);
	      cf11=getfreq(freqallele, freqfreq, data(2*i+1,C), C); 
	    } 
	    else {
	      cf00=getjointfreq2C(freqallele, condfreq2C, 
				  data(2*i  ,Cadj), Cadj, data(2*i  ,C),C);
	      cf10=getjointfreq2C(freqallele, condfreq2C, 
				  data(2*i+1,Cadj), Cadj, data(2*i  ,C),C);
	      cf01=getjointfreq2C(freqallele, condfreq2C, 
				  data(2*i  ,Cadj), Cadj, data(2*i+1,C),C);
	      cf11=getjointfreq2C(freqallele, condfreq2C, 
				  data(2*i+1,Cadj), Cadj, data(2*i+1,C),C);
	    } } 
	  
	}
      }
      
      nullprobs(7*i  ,C)=1;        //AA
      nullprobs(7*i+1,C)=cf11;     //AN1
      nullprobs(7*i+2,C)=cf00;     //NA1
      nullprobs(7*i+3,C)=cf00*cf11;//NN1
      nullprobs(7*i+4,C)=cf01;     //AN0
      nullprobs(7*i+5,C)=cf10;     //NA0
      nullprobs(7*i+6,C)=cf10*cf01;//NN0
      
    }

    //make columns C+1:RE-1
    for (j=RHSle;j<=RHSre;j++) {
      
      if (data(2*i  ,j)==0) {
	if (data(2*i+1,j)==0) {
	  
	  cf00=1;
	  cf10=1;
	  cf01=1;
	  cf11=1; }
	
	else {//data(2*i  ,j)==0),data(2*i+1,j)!=0)
	  
	  if (data(2*i  ,j+1)==0) {
	    if (data(2*i+1,j+1)==0) {
	      cf00=1;
	      cf10=1;
	      cf01=getfreq(freqallele, freqfreq, data(2*i+1,j), j);
	      cf11=getfreq(freqallele, freqfreq, data(2*i+1,j), j); 
	    } 
	    else {
	      cf00=1;
	      cf10=1;
	      cf01=getfreq(freqallele, freqfreq, data(2*i+1,j), j);
	      cf11=getcondfreq2r(freqallele, allelecounts, condfreq2r,
				 j, data(2*i+1,j), data(2*i+1,j+1), C);
	    } }
	  else {
	    if (data(2*i+1,j+1)==0) {
	      cf00=1;
	      cf10=1;
	      cf01=getcondfreq2r(freqallele, allelecounts, condfreq2r,
				 j, data(2*i+1,j), data(2*i  ,j+1), C);
	      cf11=getfreq(freqallele, freqfreq, data(2*i+1,j), j); 
	    } 
	    else {
	      cf00=1;
	      cf10=1;
	      cf01=getcondfreq2r(freqallele, allelecounts, condfreq2r,
				 j, data(2*i+1,j), data(2*i  ,j+1), C);
	      cf11=getcondfreq2r(freqallele, allelecounts, condfreq2r,
				 j, data(2*i+1,j), data(2*i+1,j+1), C);
	    } } 
	  
	}
      }
      
      else {//(data(2*i  ,j)!=0)
	if (data(2*i+1,j)==0) {
	  
	  if (data(2*i  ,j+1)==0) {
	    if (data(2*i+1,j+1)==0) {
	      cf00=getfreq(freqallele, freqfreq, data(2*i  ,j), j);
	      cf10=getfreq(freqallele, freqfreq, data(2*i  ,j), j);
	      cf01=1;
	      cf11=1;
	    } 
	    else {
	      cf00=getfreq(freqallele, freqfreq, data(2*i  ,j), j);;
	      cf10=getcondfreq2r(freqallele, allelecounts, condfreq2r,
				 j, data(2*i  ,j), data(2*i+1,j+1), C);
	      cf01=1;
	      cf11=1;
	    } }
	  else {
	    if (data(2*i+1,j+1)==0) {
	      cf00=getcondfreq2r(freqallele, allelecounts, condfreq2r,
				 j, data(2*i  ,j), data(2*i  ,j+1), C);
	      cf10=getfreq(freqallele, freqfreq, data(2*i  ,j), j);
	      cf01=1;
	      cf11=1;
	    } 
	    else {
	      cf00=getcondfreq2r(freqallele, allelecounts, condfreq2r,
				 j, data(2*i  ,j), data(2*i  ,j+1), C);
	      cf10=getcondfreq2r(freqallele, allelecounts, condfreq2r,
				 j, data(2*i  ,j), data(2*i+1,j+1), C);
	      cf01=1;
	      cf11=1;
	    } 
	  }
	  
	}
	
	else {//data(2*i  ,j)!=0),data(2*i+1,j)!=0)
	  
	  if (data(2*i  ,j+1)==0) {
	    if (data(2*i+1,j+1)==0) {
	      cf00=getfreq(freqallele, freqfreq, data(2*i  ,j), j);
	      cf10=getfreq(freqallele, freqfreq, data(2*i  ,j), j);
	      cf01=getfreq(freqallele, freqfreq, data(2*i+1,j), j);
	      cf11=getfreq(freqallele, freqfreq, data(2*i+1,j), j); 
	    } 
	    else {
	      cf00=getfreq(freqallele, freqfreq, data(2*i  ,j), j);
	      cf10=getcondfreq2r(freqallele, allelecounts, condfreq2r,
				 j, data(2*i  ,j), data(2*i+1,j+1), C);
	      cf01=getfreq(freqallele, freqfreq, data(2*i+1,j), j);
	      cf11=getcondfreq2r(freqallele, allelecounts, condfreq2r,
				 j, data(2*i+1,j), data(2*i+1,j+1), C);
	    } }
	  else {
	    if (data(2*i+1,j+1)==0) {
	      cf00=getcondfreq2r(freqallele, allelecounts, condfreq2r,
				 j, data(2*i  ,j), data(2*i  ,j+1), C);
	      cf10=getfreq(freqallele, freqfreq, data(2*i  ,j), j); 
	      cf01=getcondfreq2r(freqallele, allelecounts, condfreq2r,
				 j, data(2*i+1,j), data(2*i  ,j+1), C);
	      cf11=getfreq(freqallele, freqfreq, data(2*i+1,j), j); 
	    } 
	    else {
	      cf00=getcondfreq2r(freqallele, allelecounts, condfreq2r,
				 j, data(2*i  ,j), data(2*i  ,j+1), C);
	      cf10=getcondfreq2r(freqallele, allelecounts, condfreq2r,
				 j, data(2*i  ,j), data(2*i+1,j+1), C);
	      cf01=getcondfreq2r(freqallele, allelecounts, condfreq2r,
				 j, data(2*i+1,j), data(2*i  ,j+1), C);
	      cf11=getcondfreq2r(freqallele, allelecounts, condfreq2r,
				 j, data(2*i+1,j), data(2*i+1,j+1), C);
	    } } 
	  
	}
      }
	

      nullprobs(7*i  ,j)=1;        //AA
      nullprobs(7*i+1,j)=cf11;     //AN1
      nullprobs(7*i+2,j)=cf00;     //NA1
      nullprobs(7*i+3,j)=cf00*cf11;//NN1
      nullprobs(7*i+4,j)=cf01;     //AN0
      nullprobs(7*i+5,j)=cf10;     //NA0
      nullprobs(7*i+6,j)=cf10*cf01;//NN0
      
    }

    if (C<RE) {
      //make last column 
      if (data(2*i  ,RE)==0) {
	if (data(2*i+1,RE)==0) {
	  f0=1;
	  f1=1; }
	else {
	  f0=1;
	  f1=getfreq(freqallele, freqfreq, data(2*i+1,RE), RE); } }
      else {
	if (data(2*i+1,RE)==0) {
	  f0=getfreq(freqallele, freqfreq, data(2*i  ,RE), RE);
	  f1=1; }
	else {
	  f0=getfreq(freqallele, freqfreq, data(2*i  ,RE), RE);
	  f1=getfreq(freqallele, freqfreq, data(2*i+1,RE), RE); } }
      
      nullprobs(7*i  ,RE)=1;    //AA
      nullprobs(7*i+1,RE)=f1;   //AN1
      nullprobs(7*i+2,RE)=f0;   //NA1
      nullprobs(7*i+3,RE)=f0*f1;//NN1
      nullprobs(7*i+4,RE)=0;    //AN0
      nullprobs(7*i+5,RE)=0;    //NA0
      nullprobs(7*i+6,RE)=0;    //NN0
      
    }
  }

  for (i=0;i<nullprobs.dimx();i++) {
    if (C<LE || C>RE) 
      assert(-0.00001<nullprobs(i,C) && nullprobs(i,C)<1.00001);
    for (j=LE;j<=RE;j++)
      assert(-0.00001<nullprobs(i,j) && nullprobs(i,j)<1.00001); }
  
}

void makenullshareCnumgen(dgrid &nullprobsCnum, const igrid &data,
			  const igrid &freqallele, const dgrid &freqfreq,
			  const ivector &allelecounts,
			  const dgrid &jointfreq2, const dgrid &condfreq2C,
			  const dvector &condfreq3l,
			  const dvector &condfreq3r, 
			  const dgrid &condfreq3F,
			  const int &C,
			  const int &LE, const int &RE)
{

  //makes parts of observation probs. that do not depend on the parameters,
  //for numerator at alpha at C and beta at C-1

  int i,g, indl, indr, rwl,rwC,rwr,hap,t;
  const int numgen=data.dimx()/2;
  double mult, scratchalp, scratchbet;
  
  if (LE<C && C<RE) {
    
    for (g=0;g<numgen;g++) {
      
      //set all values to 1
      for (i=40*g;i<40*(g+1);i++) { 
	nullprobsCnum(i,0)=nullprobsCnum(i,1)=1; }

      //fill in first 4 rows, corresponding to j=0, i.e., "AA"
      for (hap=0;hap<=1;hap++) {

	//col 0 (for alpha)
	if (data(2*g+hap,C+1)!=0) {
	  mult=getfreq(freqallele, freqfreq, data(2*g+hap,C+1), C+1); }
	else mult=1.0;

	nullprobsCnum(40*g+ 2-hap,0)*=mult;	
	nullprobsCnum(40*g+     3,0)*=mult;

	//col 1 (for beta)
	if (data(2*g+hap,C-1)!=0) {
	  mult=getfreq(freqallele, freqfreq, data(2*g+hap,C-1), C-1); }
	else mult=1.0;

	nullprobsCnum(40*g+ 2-hap,1)*=mult;	
	nullprobsCnum(40*g+     3,1)*=mult;

      }

      //now fill in remaining 36 columns
      //here ind*=0 means order is AS GIVEN

      for (indl=0;indl<=1;indl++) {
	for (indr=0;indr<=1;indr++) {
	  
	  for (hap=0;hap<=1;hap++) {

	    //ordering based on function (ind* - hap)^2=0 if same, 1 otherw.
	    //needed as phase not known for genotype data
	    rwl=(hap-indl)*(hap-indl);
	    rwC=hap;
	    rwr=(hap-indr)*(hap-indr);

	    //multiply one-locus frequencies
	    if (data(2*g+(  rwr),C+1)!=0) {
	      scratchalp= getfreq(freqallele, freqfreq, 
				  data(2*g+(  rwr),C+1), C+1); }
	    else scratchalp=1;
	    if (data(2*g+(  rwl),C-1)!=0) {
	      scratchbet= getfreq(freqallele, freqfreq, 
				  data(2*g+(  rwl),C-1), C-1); }
	    else scratchbet=1;

	    for (i=0;i<2;i++) {
	      nullprobsCnum(40*g+ 4+18*indl+9*indr+4*(  rwC)+1+2*i,0)*=
		scratchalp;
	      nullprobsCnum(40*g+ 4+18*indl+9*indr+4*(  rwC)+2+i,1)*=
		scratchbet; }

	  

	    if (data(2*g+(  rwl),C-1)!=0) {
	      if (data(2*g+(  rwC),C  )!=0) {
		if (data(2*g+(  rwr),C+1)!=0) {
		  //d(C-1)!=0,d(C)!=0,d(C+1)!=0
		  scratchalp=getjointfreq3(freqallele, allelecounts,
					   condfreq3l, C-1,
					   data(2*g+rwl,C-1),data(2*g+rwC,C  ),
					   data(2*g+rwr,C+1) ); 
		  scratchbet=getjointfreq3(freqallele, allelecounts,
					   condfreq3r, C-1,
					   data(2*g+rwl,C-1),data(2*g+rwC,C  ),
					   data(2*g+rwr,C+1) ); }
		else {
		  //d(C-1)!=0,d(C)!=0,d(C+1)==0
		  scratchalp=getjointfreq2C(freqallele, condfreq2C, 
					    data(2*g+rwl,C-1),C-1,
					    data(2*g+rwC,C  ),C  ); 
		  scratchbet=getjointfreq2(freqallele, allelecounts,
					   jointfreq2, C-1,
					   data(2*g+rwl,C-1),
					   data(2*g+rwC,C  )); }
	      }
	      else {
		if (data(2*g+(  rwr),C+1)!=0) {
		  //d(C-1)!=0,d(C)==0,d(C+1)!=0
		  scratchalp=getcondfreq3F(freqallele, allelecounts,
					   condfreq3F, C, 
					   data(2*g+rwl,C-1), 
					   data(2*g+rwr,C+1),0);
		  scratchbet=getcondfreq3F(freqallele, allelecounts,
					   condfreq3F, C, 
					   data(2*g+rwl,C-1), 
					   data(2*g+rwr,C+1),1); }

		else {
		  //d(C-1)!=0,d(C)==0,d(C+1)==0
		  scratchalp=1;
		  scratchbet=getfreq(freqallele, freqfreq, 
				     data(2*g+rwl,C-1), C-1); }
	      }
	    }
	    else {
	      if (data(2*g+(  rwC),C  )!=0) {
		if (data(2*g+(  rwr),C+1)!=0) {
		  //d(C-1)==0,d(C)!=0,d(C+1)!=0
		  scratchalp=getjointfreq2(freqallele, allelecounts,
					   jointfreq2, C,
					   data(2*g+rwC,C),
					   data(2*g+rwr,C+1));
		  scratchbet=getjointfreq2C(freqallele, condfreq2C, 
					    data(2*g+rwr,C+1),C+1,
					    data(2*g+rwC,C  ),C  ); }
		else {
		  //d(C-1)==0,d(C)!=0,d(C+1)==0
		  scratchalp=getfreq(freqallele, freqfreq, 
				     data(2*g+rwC,C), C);
		  scratchbet=getfreq(freqallele, freqfreq, 
				     data(2*g+rwC,C), C); }
	      }
	      else {
		if (data(2*g+(  rwr),C+1)!=0) {
		  //d(C-1)==0,d(C)==0,d(C+1)!=0
		  scratchalp=getfreq(freqallele, freqfreq, 
				     data(2*g+rwr,C+1), C+1); 
		  scratchbet=1; }
		else {
		  //d(C-1)==0,d(C)==0,d(C+1)==0
		  scratchalp=1;
		  scratchbet=1; }
	      }
	    }

	    for (i=0;i<4;i++) {
	      nullprobsCnum(40*g+ 4+18*indl+9*indr+ 4*(1-rwC)+ i,0)*=
		scratchalp;
	      nullprobsCnum(40*g+ 4+18*indl+9*indr+ 4*(1-rwC)+ i,1)*=
		scratchbet; }
	    nullprobsCnum(40*g+ 4+18*indl+9*indr+ 8,0)*=scratchalp;
	    nullprobsCnum(40*g+ 4+18*indl+9*indr+ 8,1)*=scratchbet;
	  
	  }
	}
      }
    }

    for (g=0;g<nullprobsCnum.dimx();g++) {
      for (t=0;t<nullprobsCnum.dimy();t++)
	assert(-0.00001<nullprobsCnum(g,t) && nullprobsCnum(g,t)<1.00001); }
  }
}		      


void makenullshareCdengen(dgrid &nullprobsCden, const igrid &data, 
			  const igrid &freqallele,
			  const dgrid &freqfreq, const int &C, 
			  const int &LE, const int &RE)
{
  //makes parts of observation probs. that do not depend on the parameters,
  //for denominator at alpha at C and beta at C-1
  
  int g,side,E;
  const int numgen=data.dimx()/2;
  
  if ( (LE<C && C<RE) || RE<C) {
    
    for (g=0;g<numgen;g++) {
      
      for (side=0;side<=1;side++) {

	if (side==0) {
	  if (LE<C && C<RE) E=C+1; }
	else {//side==1
	  if (LE<C && C<RE) E=C-1;
	  else E=RE; }
	
	if ( (LE<C && C<RE) || side==1) {
	  
	  nullprobsCden(7*g  ,side)=1;      

	  if (data(2*g  ,E)!=0) {
	    if (data(2*g+1,E)!=0) {
	      
	      nullprobsCden(7*g+1,side)=
		getfreq(freqallele, freqfreq, data(2*g+1,E), E);
	      nullprobsCden(7*g+2,side)=
		getfreq(freqallele, freqfreq, data(2*g  ,E), E);
	      nullprobsCden(7*g+3,side)=
		getfreq(freqallele, freqfreq, data(2*g  ,E), E)*
		getfreq(freqallele, freqfreq, data(2*g+1,E), E);
	      
	    }
	    else {
	      
	      nullprobsCden(7*g+1,side)=1;
	      nullprobsCden(7*g+2,side)=nullprobsCden(7*g+3,side)=
		getfreq(freqallele, freqfreq, data(2*g  ,E), E);
	      
	    } 
	  }
	  else {
	    if (data(2*g+1,E)!=0) {
	      
	      nullprobsCden(7*g+2,side)=1;
	      nullprobsCden(7*g+1,side)=nullprobsCden(7*g+3,side)=
		getfreq(freqallele, freqfreq, data(2*g+1,E), E);
	      
	    }
	    else {
	      
	      nullprobsCden(7*g+1,side)=nullprobsCden(7*g+2,side)=
		nullprobsCden(7*g+3,side)=1;
	      
	    } 
	  }
	  
	  //switch ordering for alpha but not for beta
	  nullprobsCden(7*g+4,0)=nullprobsCden(7*g+2,0);
	  nullprobsCden(7*g+5,0)=nullprobsCden(7*g+1,0);
	  nullprobsCden(7*g+6,0)=nullprobsCden(7*g+3,0);
	
	  nullprobsCden(7*g+4,1)=nullprobsCden(7*g+1,1);
	  nullprobsCden(7*g+5,1)=nullprobsCden(7*g+2,1);
	  nullprobsCden(7*g+6,1)=nullprobsCden(7*g+3,1);
	}
      }
    }
  }

}

void makenullshareCltREgen(dvector &nullprobsCltRE, const igrid &data,
			   const igrid &freqallele, const dgrid &freqfreq,
			   const ivector &allelecounts,
			   const dgrid &jointfreq2C, 
			   const int &C,
			   const int &RE)
{
  //if (RE<C), need another observation probability-- this assembles
  //the part that doesn't depend on the parameters

  int i,g, indl, rwl,rwC,hap;
  const int numgen=data.dimx()/2;
  double mult, scratchbet;
  
  if (RE<C) {

    for (g=0;g<numgen;g++) {
      
      //set all values to 1
      for (i=14*g;i<14*(g+1);i++) { 
	nullprobsCltRE(i)=1; }

      //fill in first 4 rows, corresponding to j=0, i.e., "AA"
      for (hap=0;hap<=1;hap++) {

	if (data(2*g+hap,RE)!=0) {
	  mult=getfreq(freqallele, freqfreq, data(2*g+hap,RE), RE); }
	else mult=1;

	nullprobsCltRE(14*g+ 2-hap)*=mult;	
	nullprobsCltRE(14*g+     3)*=mult;

      }

      //now fill in remaining 10 columns
      //here ind*=0 means order is AS GIVEN
      
      for (indl=0;indl<=1;indl++) {

	for (hap=0;hap<=1;hap++) {

	  //ordering based on function (ind* - hap)^2=0 if same, 1 otherw.
	  //needed as phase not known for genotype data
	  rwl=(hap-indl)*(hap-indl);
	  rwC=hap;
	  
	  
	  if (data(2*g+(  rwl),RE)!=0) {
	    scratchbet= getfreq(freqallele, freqfreq, 
				data(2*g+(  rwl),RE), RE); }
	  else {
	    scratchbet=1; }

	  nullprobsCltRE(14*g+4+5*indl+2*(  rwC)+1)*=
	    scratchbet; 
	
	  if (data(2*g+(  rwl),RE)!=0) {
	    if (data(2*g+(  rwC),C  )!=0) {
	      //d(C-1)!=0,d(C)!=0
	      scratchbet=getjointfreq2C(freqallele, jointfreq2C, 
					data(2*g+rwl,RE),RE,
					data(2*g+rwC,C  ),C  ); }
	    else {
	      //d(C-1)!=0,d(C)==0
	      scratchbet=getfreq(freqallele, freqfreq, 
				 data(2*g+rwl,RE), RE); }
	  }
	  else {
	    if (data(2*g+(  rwC),C  )!=0) {
	      //d(C-1)==0,d(C)!=0
	      scratchbet=getfreq(freqallele, freqfreq, 
				 data(2*g+rwC,C), C); }
	    else {
	      //d(C-1)==0,d(C)==0
	      scratchbet=1; }
	  }

	  for (i=0;i<2;i++) {
	    nullprobsCltRE(14*g+4+5*indl+ 2*(1-rwC)+ i)*=
	      scratchbet; 
	    	    
	  }
	  nullprobsCltRE(14*g+4+5*indl+ 4)*=scratchbet; 

	}
      }
    }
  }
}

void makemodlsharegen(dgrid &obsprobs, const dgrid &nullprobs, 
		      const igrid &data,
		      const dvector &mutmatch, dvector &mutnomatch,
		      const ivector &anchap, 
		      const int &C, const int &LE, const int &RE)
{
  //assembles observation probabilities, including portion that depends
  //on 1/tau through m(l,tau,i,j)

  int g,t;
  const int numgen=data.dimx()/2;

  //set obsprobs equal to null probs
  obsprobs=nullprobs;
  
  //now include mutation bits
  for (g=0;g<numgen;g++) {
 
    //mutation not allowed at C
    if (C<LE || RE<C) {
      if (data(2*g  ,C)!=anchap(C)) {
	obsprobs(7*g  ,C)=obsprobs(7*g+1,C)=obsprobs(7*g+4,C)=0; }
      if (data(2*g+1,C)!=anchap(C)) {
	obsprobs(7*g  ,C)=obsprobs(7*g+2,C)=obsprobs(7*g+5,C)=0; }
    }

    for (t=LE;t<=RE;t++) {

      if (t!=C)
	{

	  if (data(2*g,t)!=0) {
	    if (data(2*g,t)==anchap(t)) {
	      obsprobs(7*g  ,t)*=mutmatch(t);
	      obsprobs(7*g+1,t)*=mutmatch(t);
	      obsprobs(7*g+4,t)*=mutmatch(t); }
	    else {
	      obsprobs(7*g  ,t)*=mutnomatch(t);
	      obsprobs(7*g+1,t)*=mutnomatch(t);
	      obsprobs(7*g+4,t)*=mutnomatch(t); }
	  }
	  
	  if (data(2*g+1,t)!=0) {
	    if (data(2*g+1,t)==anchap(t)) {
	      obsprobs(7*g  ,t)*=mutmatch(t);
	      obsprobs(7*g+2,t)*=mutmatch(t);
	      obsprobs(7*g+5,t)*=mutmatch(t); }
	    else {
	      obsprobs(7*g  ,t)*=mutnomatch(t);
	      obsprobs(7*g+2,t)*=mutnomatch(t);
	      obsprobs(7*g+5,t)*=mutnomatch(t); }
	  }
	  
	}
    }
  }
}

void makemodlshareCnumgen(dgrid &obsprobsCnum, const dgrid &nullprobsCnum,
			  const igrid &data, const ivector &anchap, 
			  const dvector &mutmatch, const dvector &mutnomatch,
			  const int &C, const int &LE, const int &RE)
{
  //assembles observation probs. for numerator of alpha and beta and entry
  //linking left and right sides; includes portion that's dependent on tau
  //through m(l,tau,i,j)

  //assumes no mutation allowed at center; we condition on this(?)
  int i,g,hap, indl, indr, rwl, rwr,t;
  const int numgen=data.dimx()/2;
  double valuer[2];
  double valuel[2];

  if (LE<C && C<RE) {
    
    for (g=0;g<numgen;g++) {
      
      //set all values to nullprobsCnum
      for (i=40*g;i<40*(g+1);i++) { 
	obsprobsCnum(i,0)=nullprobsCnum(i,0);
	obsprobsCnum(i,1)=nullprobsCnum(i,1); 

      }

      //mutation not allowed at C; if allele doesn't equal ancestral
      //allele, and state is A, prob is 0
      //NOTE: treat missing as not matching (at C)

      for (hap=0;hap<=1;hap++) {
	if (data(2*g+hap,C)!=anchap(C)) {
	  for (i=0;i<4;i++) {
	    obsprobsCnum(40*g+i,0)=obsprobsCnum(40*g+i,1)=0; 
	    for (indl=0;indl<=1;indl++) {
	      for (indr=0;indr<=1;indr++) {
		obsprobsCnum(40*g+4+ 18*indl+9*indr +4*hap+i,0)=
		  obsprobsCnum(40*g+4+ 18*indl+9*indr +4*hap+i,1)=0; } } } } }
    
      //now multiply by approprate mutation prob. for flanking markers
      //(here ind* =0 means order is AS GIVEN (* =l or =r))

      for (i=0;i<2;i++) {
	//for alpha object
	if (data(2*g+i,C+1)==anchap(C+1)) {
	  valuer[i]=mutmatch(C+1); }
	else {  
	  if (data(2*g+i,C+1)!=0) {
	    valuer[i]=mutnomatch(C+1); }
	  else valuer[i] =1; } 

	//for beta object
	if (data(2*g+i,C-1)==anchap(C-1)) {
	  valuel[i]=mutmatch(C-1); }
	else { 
	  if (data(2*g+i,C-1)!=0) {
	    valuel[i]=mutnomatch(C-1); }
	  else valuel[i] =1; } 
      }

      //fill in first four lines
      obsprobsCnum(40*g  ,0)*= valuer[0]*valuer[1];
      obsprobsCnum(40*g  ,1)*= valuel[0]*valuel[1]; 
      obsprobsCnum(40*g+1,0)*= valuer[0];
      obsprobsCnum(40*g+1,1)*= valuel[0];
      obsprobsCnum(40*g+2,0)*= valuer[1];
      obsprobsCnum(40*g+2,1)*= valuel[1];

      //now complete remaining 36 lines

      for (indl=0;indl<=1;indl++) {
	for (indr=0;indr<=1;indr++) {
	  
	  for (hap=0;hap<=1;hap++) {

	    //ordering based on function (ind* - hap)^2=0 if same, 1 otherw.
	    //needed as phase not known for genotype data
	    rwl=(hap-indl)*(hap-indl);
	    rwr=(hap-indr)*(hap-indr);

	    for (i=0;i<=1;i++) {

	      obsprobsCnum(40*g+ 4+18*indl+9*indr+4*hap+2*i,0)*=
		valuer[rwr];
	      obsprobsCnum(40*g+ 4+18*indl+9*indr+4*hap+i,1)*=
		valuel[rwl]; 

	    }
	  }
	}
      }

    }
  
    for (g=0;g<obsprobsCnum.dimx();g++) {
      for (t=0;t<obsprobsCnum.dimy();t++)
	assert(-0.00001<obsprobsCnum(g,t) && obsprobsCnum(g,t)<1.00001); }
  }
}
		
void makemodlshareCdengen(dgrid &obsprobsCden, const dgrid &nullprobsCden,
			  const igrid &data, const ivector &anchap, 
			  const dvector &mutmatch, const dvector &mutnomatch,
			  const int &C, const int &LE, const int &RE)
{
  //assembles observation probs. for denominator of alpha and beta and entry
  //linking left and right sides; includes portion that's dependent on tau
  //through m(l,tau,i,j)

  int i,g,E;
  double value;
  const int numgen=data.dimx()/2;

  if ( (LE<C && C<RE) || RE<C) {
    
    for (g=0;g<numgen;g++) {
      
      //set all values to nullprobsCden
      for (i=7*g;i<7*(g+1);i++) { 
	if (LE<C && C<RE) obsprobsCden(i,0)=nullprobsCden(i,0);
	obsprobsCden(i,1)=nullprobsCden(i,1); }

      if (LE<C && C<RE) {
	//make column for alpha
	if (data(2*g  ,C+1)==anchap(C+1)) {
	  value=mutmatch(C+1); }
	else { 
	  if (data(2*g  ,C+1)!=0) {
	    value=mutnomatch(C+1); }
	  else value = 1.0; } 
	
	obsprobsCden(7*g  ,0)*=value;
	obsprobsCden(7*g+1,0)*=value;
	obsprobsCden(7*g+5,0)*=value;
	
	if (data(2*g+1,C+1)==anchap(C+1)) {
	  value=mutmatch(C+1); }
	else { 
	  if (data(2*g+1,C+1)!=0) {
	    value=mutnomatch(C+1); }
	  else value = 1.0; } 
      
	obsprobsCden(7*g  ,0)*=value;
	obsprobsCden(7*g+2,0)*=value;
	obsprobsCden(7*g+4,0)*=value;
      }

      //now for beta
      if (LE<C && C<RE) E=C-1;
      else E=RE;

      if (data(2*g  ,E)==anchap(E)) {
	value=mutmatch(E); }
      else { 
	if (data(2*g  ,E)!=0) {
	  value=mutnomatch(E); }
	else value = 1.0; } 

      obsprobsCden(7*g  ,1)*=value;//was 0,1,4
      obsprobsCden(7*g+1,1)*=value;
      obsprobsCden(7*g+4,1)*=value;
	
	if (data(2*g+1,E)==anchap(E)) {
	  value=mutmatch(E); }
	else { 
	  if (data(2*g+1,E)!=0) {
	    value=mutnomatch(E); }
	  else value = 1.0; } 

	obsprobsCden(7*g  ,1)*=value;//was 0,2,5
	obsprobsCden(7*g+2,1)*=value;
	obsprobsCden(7*g+5,1)*=value;
	
    }
  }
}

void makemodlshareCltREgen(dvector &obsprobsCltRE, 
			   const dvector &nullprobsCltRE,
			   const igrid &data, const ivector &anchap, 
			   const dvector &mutmatch, const dvector &mutnomatch,
			   const int &C, const int &RE)
{
  //assembles extra observation probability for when RE<C

  //assumes no mutation allowed at center; we condition on this(?)
  int i,g,hap, indl,rwl;
  const int numgen=data.dimx()/2;
  double valuel[2];

  if (RE<C) {
    for (g=0;g<numgen;g++) {
      
      //set all values to nullprobsCnum
      for (i=14*g;i<14*(g+1);i++) { 
	obsprobsCltRE(i)=nullprobsCltRE(i); }

      //mutation not allowed at C; if allele doesn't equal ancestral
      //allele, and state is A, prob is 0
      //NOTE: treat missing as not matching (at C)

      for (hap=0;hap<=1;hap++) {
	if (data(2*g+hap,C)!=anchap(C)) {
	  for (i=0;i<4;i++) {
	    obsprobsCltRE(14*g+i)=0; }
	  for (i=0;i<2;i++) {
	    for (indl=0;indl<=1;indl++) {
	      obsprobsCltRE(14*g+4+ 5*indl+ 2*hap+i)=0; } } } }

      //now multiply by approprate mutation prob. for flanking markers
      //(here ind* =0 means order is AS GIVEN (* =l or =r))

      for (i=0;i<2;i++) {
	//for beta object
	if (data(2*g+i,RE)==anchap(RE)) {
	  valuel[i]=mutmatch(RE); }
	else { 
	  if (data(2*g+i,RE)!=0) {
	    valuel[i]=mutnomatch(RE); }
	  else valuel[i] =1.0; } 
      }

      //fill in first four lines
      obsprobsCltRE(14*g  )*= valuel[0]*valuel[1]; 
      obsprobsCltRE(14*g+1)*= valuel[0];
      obsprobsCltRE(14*g+2)*= valuel[1]; 

      //now complete remaining 10 lines

      for (indl=0;indl<=1;indl++) {

	for (hap=0;hap<=1;hap++) {

	  //ordering based on function (ind* - hap)^2=0 if same, 1 otherw.
	  //needed as phase not known for genotype data
	  rwl=(hap-indl)*(hap-indl);
	  
	  obsprobsCltRE(14*g+ 4+5*indl+2*hap)*=
	    valuel[rwl]; 
	}
      }
    }
  }
}

void makeinitgen(dgrid &initstore, dvector &scratchvec,
		 const dvector &x, const double &tau, 
		 const double &p, const double &freqC,
		 const int &C, const int &LE, const int &RE)
{

  int i,t;
  
  for (t=0;t<=3;t++) {

    if (t==0) {//for left end
      //initial distribution of chain on LHS
      if (C<LE) scratchvec(0)= freqC*(1-p); //A 
      else scratchvec(0)=freqC * (1-p) * exp( -tau * x[LE]);//A
      scratchvec(1)=1-scratchvec(0);//N 
    }
    else {
      if (t==1) {
	//make unconditional distribution of chain at C
	scratchvec(0)= freqC * (1-p);
	scratchvec(1)= 1-scratchvec(0); 
      }
      else {
	if (t==2) {
	  //initial distribution of chain on right end
	  if (RE<C) scratchvec(0)= freqC * (1-p); //A 
	  else scratchvec(0)=freqC * (1-p)* exp( -tau * x[RE]);//A
	  scratchvec(1)=1-scratchvec(0);//N 
	}
	else {//t==3, dist. at C-1 or RE (if C<RE)
	  if (RE<C || (LE<C &&C <RE)) {
	    if (RE<C) scratchvec(0)=freqC * (1-p)*exp( -tau * x[RE]); //A 
	    else scratchvec(0)=freqC * (1-p)* exp( -tau * x[C-1]);//A
	    scratchvec(1)=1-scratchvec(0);//N 
	  }
	}
      }
    }

    initstore(0,t)=scratchvec(0)*scratchvec(0);//AA
    initstore(1,t)=scratchvec(0)*scratchvec(1);//AN0
    initstore(2,t)=scratchvec(1)*scratchvec(0);//NA0
    initstore(3,t)=scratchvec(1)*scratchvec(1);//NN0
    
    //chain anchored at ends
    if ( (t==0 && C<LE) || (t==1 && (C<LE || RE<C) ) ||
	 (t==2 && RE<C) || (t==3 && (C-1!=LE && LE<C) ) ) {
      for (i=1;i<=3;i++) {
	initstore(i,t)*=0.5; } }
    else {
      if ( (t==1) && (LE<C && C<RE) ) {
	for (i=1;i<=3;i++) {
	  initstore(i,t)*=0.25;
	}
      }
    }
  }
}

void maketransgen(dgrid &transstore, dvector &scratchvec, const dvector &x, 
		  const double &tau, const double &p,
		  const double &freqC,
		  const int &C, const int &LE, const int &RE)
{
  //makes transition probabilities for HMM, assuming genotypes known
  //for affecteds, for transitions in the direction of C
  
  int t,i,j,k,l;
  int LHSle, RHSre;
  
  if (LE<C && C<RE) {
    LHSle=LE; 
    RHSre=RE-1; }
  else {
    if (C<LE) {
      LHSle=C;
      RHSre=RE-1; }
    else {//RE<C
      LHSle=LE; 
      RHSre=RE; }
  }

  for (t=LHSle;t<=RHSre;t++) {
    
    if (  ( ((LE<C && C<RE) && t<=C-1) || (RE<C && t<=RE-1) ) ||
	  (RE<C && t==RE) ||
	  (C<LE && t==C) ||
	  ( ((LE<C && C<RE) && C<=t) || (C<LE && LE<=t) )  ) {
      
      //LHS
      if ( ((LE<C && C<RE) && t<=C-1) || (RE<C && t<=RE-1) ) {
	scratchvec(0)=1;//A->A  
	scratchvec(1)=0;//A->N 
	scratchvec(3)= (1 - freqC * (1-p) * exp( -tau * x[t+1]) )/
	  (1 - freqC * (1-p) * exp( -tau * x[t]) );//N->N
	scratchvec(2)=1-scratchvec(3);//N->A
      }
      else {
	if (RE<C && t==RE) {
	  scratchvec(0)=1;//A->A  
	  scratchvec(1)=0;//A->N
	  scratchvec(3)= 
	    (1 - freqC * (1-p))/
	    (1 - freqC * (1-p) * exp( -tau * x[RE]) );//N->N
	  scratchvec(2)=1-scratchvec(3); //N->A
	}
	else {
	  //RHS
	  if (C<LE && t==C) {
	    scratchvec(0)=1;//A<-A  
	    scratchvec(2)=0;//N<-A
	    scratchvec(3)= 
	      (1 - freqC * (1-p))/
	      (1 - freqC * (1-p) * exp( -tau * x[LE]) );//N<-N
	    scratchvec(1)=1-scratchvec(3);//A<-N
	  }
	  else {
	    if ( ((LE<C && C<RE) && C<=t) || (C<LE && LE<=t) ) {
	      scratchvec(0)=1;//A<-A  
	      scratchvec(2)=0;//N<-A
	      scratchvec(3)= (1 - freqC * (1-p) * exp( -tau * x[t]) )/
		(1 - freqC * (1-p) * exp( -tau * x[t+1]) );//N<-N
	      scratchvec(1)=1-scratchvec(3); //A<-N
	    }
	  }
	}
      }


      for (l=0;l<2;l++) {
	for (k=0;k<2;k++) {
	  for (i=0;i<2;i++) {
	    for (j=0;j<2;j++) {
	      
	      transstore(8*l+4*k+2*i+j,t)=
		scratchvec(2*l+i)*scratchvec(2*k+j);
	      
	    }
	  }
	}
      }
      

      //multiply stuff by 0.25 (accounts for phase inicators in statespace)
      if ( (LE<C && C<RE) && (t==C-1 || t==C) ) {
	if (t==C-1) {
	  for (i=0;i<4;i++) {
	    for (j=1;j<4;j++) {
	      transstore(4*i+j,t)*=0.25; } } }
	else {
	  for (i=1;i<4;i++) {
	    for (j=0;j<4;j++) {
	      transstore(4*i+j,t)*=0.25; } } } }
      else {
	//multiply stuff by 0.5
	if ( ((LE<C && C<RE) && (t<C-1)) || (RE<C && t<=RE) ) {
	  for (i=0;i<4;i++) {
	    for (j=1;j<4;j++) {
	      transstore(4*i+j,t)*=0.5; } } }
	else {
	  if ( ((LE<C && C<RE) && (C<t)) || (C<LE && C<=t) ) {
	    for (i=1;i<4;i++) {
	      for (j=0;j<4;j++) {
		transstore(4*i+j,t)*=0.5; } } } } }
    }
  }
}

void maketransCgen(dgrid &transstoreC, dvector &scratchvec,
		   const dvector &x, 
		   const double &tau, const double &p, 
		   const int &C, const int &LE, const int &RE)
{
  //makes transition probs for jumps away from C, assuming genotypes
  //known for affecteds

  int col,i,j,k,l;  

  for (col=0;col<=1;col++) {
    if (col==0) {
      if (LE<C)  
	{
	  scratchvec(0)= exp( -tau * x[C-1]);//A<-A  
	  scratchvec(1)=0;//A<-N 
	  scratchvec(2)= 1- scratchvec(0);//N<-A
	  scratchvec(3)=1;//N<-N
	}
      else scratchvec(0)=scratchvec(1)=
	     scratchvec(2)=scratchvec(3)= 0;
    }
    else {
      if (C<RE)  
	{
	  scratchvec(0)= exp( -tau * x[C+1]);//A->A  
	  scratchvec(1)=1- scratchvec(0);//A->N 
	  scratchvec(2)=0;//N->A
	  scratchvec(3)=1;//N->N
	}
      else 
	{
	  scratchvec(0)= exp( -tau * x[RE]);//A<-A  
	  scratchvec(1)=0;//A<-N 
	  scratchvec(2)= 1- scratchvec(0);//N<-A
	  scratchvec(3)=1;//N<-N 
      }
    }
   
    for (l=0;l<2;l++) {
      for (k=0;k<2;k++) {
	for (i=0;i<2;i++) {
	  for (j=0;j<2;j++) {
	    
	    transstoreC(8*l+4*k+2*i+j,col)=
	      scratchvec(2*l+i)*scratchvec(2*k+j);

	  }
	}
      }
    }

    //account for phase indicators in statespace
    if (col==0) {
      if (C-1!=LE) {
	for (i=1;i<4;i++) {
	  for (j=0;j<4;j++) {
	    transstoreC(4*i+j,col)*=0.5; } } } }
    else {
      if (C<RE) {
	if (C+1!=RE) {
	  for (i=0;i<4;i++) {
	    for (j=1;j<4;j++) {
	      transstoreC(4*i+j,col)*=0.5; } } } }
      else {
	for (i=1;i<4;i++) {
	  for (j=0;j<4;j++) {
	    transstoreC(4*i+j,col)*=0.5; } } }
    }
  }
}

void makealphagen(dgrid &alpha, const dgrid &initstore, 
		  const dgrid &transstore, const dgrid &transstoreC,
		  const dgrid &obsprobs, const dgrid &obsprobsCnum,
		  const dgrid &obsprobsCden, 
		  const int &C, const int &LE, const int &RE)
{

  //makes alpha object, when genotypes from affecteds are available

  int i,j,k,g,t;
  const int numgen=alpha.dimx()/13;
  int LHSle, LHSre, RHSle, RHSre, nstatesi, nstatesj, nstatesk,
    tsi,tsj,tsk,indl,indr,multi, multk,adj;
  double sum,sum1;

  //these are endpoints for t, rather than t+1
  if (LE<C && C<RE) {
    LHSle=LE; LHSre=C-2;
    RHSle=C;  RHSre=RE-1; }
  else {
    if (C<LE) {
      LHSle=0,LHSre=-1;
      RHSle=LE;  RHSre=RE-1; }
    else {//RE<C
      LHSle=LE; LHSre=RE-1;
      RHSle=0,RHSre=-1; }
  }

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

    if (C<LE) {
      //make column(s) C
      for (i=0;i<7;i++) {
	alpha(13*g+i,C)=obsprobs(7*g+i,C); }
      
      //make column LE
      nstatesi=7;
      nstatesj=7;
      for (j=0;j<nstatesj;j++) {

	sum=0;
	if (j<4) tsj=0; else tsj=3;//needed to get transstore entry
	for (i=0;i<nstatesi;i++) {
	  if (i<4) tsi=0;
	  else {
	    if (i>9) tsi=9; //needed to get transstore entry
	    else {
	      if (i>6) tsi=6; 
	      else tsi=3; } } 

	  if ( (3<i && i<7) || (9<i) ) {
	    if (j==1 || j==4) adj=1;
	    else {
	      if (j==2 || j==5) adj=-1; 
	      else adj=0; } }
	  else adj=0;

	  sum+=
	    alpha(13*g+i,C)*transstore(4*(i-tsi)+(j-tsj+adj),C)*
	    obsprobs(7*g+j,LE); 
	}
	alpha(13*g+j,LE)=sum; 
      }      
    }

    else {//(LE<C)
      //make column(s) LE
      for (i=0;i<4;i++) {
	alpha(13*g+i,LE)=initstore(i,0)*obsprobs(7*g+i,LE); 

      }
    }

    //make column(s) LE+1:C-1 or LE+1:RE
    for (t=LHSle;t<=LHSre;t++) {
      if (t==LE  && LE<C) nstatesi=4; else nstatesi=7;
      for (j=0;j<7;j++) {

	sum=0;
	if (j<4) tsj=0; else tsj=3;
	for (i=0;i<nstatesi;i++) {
	  if (i<4) tsi=0; else tsi=3;//needed to get transstore entry

	  if (j>3) {
	    if (i==1 || i==4) adj=1;
	    else {
	      if (i==2 || i==5) adj=-1; 
	      else adj=0; } }
	  else adj=0;

	  sum+=
	    alpha(13*g+i,t)*transstore(4*(i-tsi+adj)+(j-tsj),t)*
	    obsprobs(7*g+j,t+1); 
	}
	alpha(13*g+j,t+1)=sum; 
      } 
    }
    
    if (LE<C && C<RE) {
      //make column(s) C (see notes for explanation about why this is messy)

      if (C-1==LE) { nstatesi=4; multi=1; }
      else { nstatesi=7; multi=2; }
      if (C+1==RE) { nstatesk=4; multk=1; }
      else { nstatesk=7; multk=2; }


      //make numerators
      //j=0
      sum=0;
      for (i=0;i<nstatesi;i++) {
	if (i<4) tsi=0; else tsi=3;
	sum1=0;
	for (k=0;k<nstatesk;k++) {
	  if (k<4) tsk=0; else tsk=3;
	  sum1+=transstoreC(4*(0    )+(k-tsk),1)*
	    obsprobsCnum(40*g+k-tsk,0); 	  
	}
	sum+=alpha(13*g+i,C-1)*transstore(4*(i-tsi)+(0    ),C-1)*sum1; 
      }
      alpha(13*g  ,C)=sum;

      //j=1:12
      for (indl=0;indl<=1;indl++) {
	for (indr=0;indr<=1;indr++) {
	  
	  //j=1,4,7,10
	  j=indl*6+indr*3+ 1;
	  alpha(13*g+j,C)=
	    ( alpha(13*g+ 1+indl,C-1) + alpha(13*g+ 4+indl,C-1) )*
	    transstore(4*1+1,C-1)*multk*
	    ( transstoreC(4*1+1,1)*
	      obsprobsCnum(40*g+18*indl+9*indr +4  ,0)+
	      transstoreC(4*1+3,1)*
	      obsprobsCnum(40*g+18*indl+9*indr +4+1,0) )+
	    ( alpha(13*g+ 3     ,C-1) + alpha(13*g+ 6     ,C-1) )*
	    transstore(4*3+1,C-1)*multk*
	    ( transstoreC(4*1+1,1)*
	      obsprobsCnum(40*g+18*indl+9*indr +4+2,0)+
	      transstoreC(4*1+3,1)*
	      obsprobsCnum(40*g+18*indl+9*indr +4+3,0) );

	  //j=2,5,8,11
	  j=indl*6+indr*3+ 2;
	  alpha(13*g+j,C)=
	    ( alpha(13*g+ 2-indl,C-1) + alpha(13*g+ 5-indl,C-1) )*
	    transstore(4*2+2,C-1)*multk*
	    ( transstoreC(4*2+2,1)*
	      obsprobsCnum(40*g+18*indl+9*indr +4+4,0)+
	      transstoreC(4*2+3,1)*
	      obsprobsCnum(40*g+18*indl+9*indr +4+5,0) )+
	    ( alpha(13*g+ 3     ,C-1) + alpha(13*g+ 6     ,C-1) )*
	    transstore(4*3+2,C-1)*multk*
	    ( transstoreC(4*2+2,1)*
	      obsprobsCnum(40*g+18*indl+9*indr +4+6,0)+
	      transstoreC(4*2+3,1)*
	      obsprobsCnum(40*g+18*indl+9*indr +4+7,0) );

	  //j=3,6,9,12
	  j=indl*6+indr*3+ 3;
	  alpha(13*g+j,C)=
	    (alpha(13*g+ 3      ,C-1) + alpha(13*g+ 6     ,C-1) )*
	    transstore(4*3+3,C-1)*multk*
	    ( transstoreC(4*3+3,1)*
	      obsprobsCnum(40*g+18*indl+9*indr +4+8,0) );

	}
      }
    
      //make denominators

      //j=0
      sum=0;
      sum+=initstore(0,1)*transstoreC(4*0+0,1)*obsprobsCden(7*g  ,0);
      for (k=1;k<4;k++) {
	sum+=multk*initstore(0,1)*transstoreC(4*0+k,1)*
	  obsprobsCden(7*g+k,0); }
      
      if (sum==0) {
	if (alpha(13*g,C)!=0) {
	  errorfs << "ERROR: alpha(13*" << g << "+" << 0 << ",C) = "
	       << "P(A|B) = P(A&B)/P(B) =!0/0" << endl; exit(1); } }
      else alpha(13*g,C)/=sum; 
      
      for (indl=0;indl<=1;indl++) {
	for (indr=0;indr<=1;indr++) {

	  for (j=1;j<=3;j++) {
	    if (j==1) {
	      //j=1,4,7,10
	      //j=6*indl+3*indr +1;
	      sum=initstore(1,1)*multk*
		( transstoreC(4*1+1,1)*obsprobsCden(7*g+ 1+3*indr,0)+
		  transstoreC(4*1+3,1)*obsprobsCden(7*g+ 3+3*indr,0) );
	    }
	    else {
	      if (j==2) {
		//j=2,5,8,11
		//j=6*indl+3*indr +2;
		sum=initstore(2,1)*multk*
		  ( transstoreC(4*2+2,1)*obsprobsCden(7*g+ 2+3*indr,0)+
		    transstoreC(4*2+3,1)*obsprobsCden(7*g+ 3+3*indr,0) );
	      }
	      else {
		//j=3,6,9,12
		//j=6*indl+3*indr +3;
		sum=initstore(3,1)*multk*
		  ( transstoreC(4*3+3,1)*obsprobsCden(7*g+ 3+3*indr,0) );
	      }
	    }
	  
	    if (sum==0) {
	      if (alpha(13*g+6*indl+3*indr+j,C)!=0) {
		errorfs << "ERROR: alpha(13*" << g << "+" 
		     << 6*indl+3*indr+j << ",C) = "
		     << "P(A|B) = P(A&B)/P(B) =!0/0" << endl; exit(1); } }
	    else alpha(13*g+6*indl+3*indr+j,C)/=sum; 
	
	  }
	}
      }
    }
  
    //make column(s) C+1:RE or LE+1:RE
    for (t=RHSle;t<=RHSre;t++) {
      if (t==C && (LE<C && C<RE) ) nstatesi=13; else nstatesi=7;
      if (t==RE-1 && C<RE) nstatesj=4; else nstatesj=7;
      for (j=0;j<nstatesj;j++) {

	sum=0;
	if (j<4) tsj=0; else tsj=3;//needed to get transstore entry
	for (i=0;i<nstatesi;i++) {
	  if (i<4) tsi=0;
	  else {
	    if (i>9) tsi=9; //needed to get transstore entry
	    else {
	      if (i>6) tsi=6; 
	      else tsi=3; } } 

	  if ( (3<i && i<7) || (9<i) ) {
	    if (j==1 || j==4) adj=1;
	    else {
	      if (j==2 || j==5) adj=-1; 
	      else adj=0; } }
	  else adj=0;

	  sum+=
	    alpha(13*g+i,t)*transstore(4*(i-tsi)+(j-tsj+adj),t)*
	    obsprobs(7*g+j,t+1); 
	}
	alpha(13*g+j,t+1)=sum; 
      }
    }

    if (RE<C) {
      //make column C 
      nstatesi=7;
      nstatesj=7;
      for (j=0;j<nstatesj;j++) {

	sum=0;
	if (j<4) tsj=0; else tsj=3;//needed to get transstore entry
	for (i=0;i<nstatesi;i++) {
	  if (i<4) tsi=0; else tsi=3;
	  
	  if (j>3) {
	    if (i==1 || i==4) adj=1;
	    else {
	      if (i==2 || i==5) adj=-1; 
	      else adj=0; } }
	  else adj=0;

	  sum+=
	    alpha(13*g+i,RE)*transstore(4*(i-tsi+adj)+(j-tsj),RE)*
	    obsprobs(7*g+j,C); 

	}

	if (initstore(j-tsj,1)==0) {
	  if (sum!=0) {
	    errorfs << "ERROR: alpha(13*" << g << "+" 
		 << j << ",C) = "
		 << "P(A|B) = P(A&B)/P(B) =!0/0" << endl; 
	    exit(1);
	  } }
	else sum/=initstore(j-tsj,1);
	alpha(13*g+j,C)=sum; 

      }
    }

  }
    
  //check that all entries are between 0 and 1
  for (i=0;i<alpha.dimx();i++) {
    if (C<LE || RE<C) { 
      if ( !(-0.00001<alpha(i,C) && alpha(i,C)<1.00001) ) {
	errorfs << "ERROR: alpha(" << i << "," << C << ") = " 
	     << alpha(i,C) << endl;
	exit(1);
      }
    }
    for (j=LE;j<=RE;j++) {
      if ( !(-0.00001<alpha(i,j) && alpha(i,j)<1.00001) ) {
	errorfs << "ERROR: alpha(" << i << "," << j << ") = " 
	     << alpha(i,j) << endl;
	exit(1);
      } 
    }
  }
}

void makebetagen(dgrid &beta, const dgrid &initstore, 
		 const dgrid &transstore, const dgrid &transstoreC,
		 const dgrid &obsprobs, const dgrid &obsprobsCnum,
		 const dgrid &obsprobsCden, const dvector &obsprobsCltRE,
		 const int &C, const int &LE, const int &RE)
{

  //makes beta object, assuming genotype data for affecteds

  int i,j,g,t;
  const int numgen=beta.dimx()/13;
  int LHSle, LHSre, RHSle, RHSre, nstatesi, nstatesj, nstatesk,tsi,tsj;
  double sum,sum2,multi,multk;
  int indl,indr,adj;

  //these are endpoints for t, rather than t+1
  if (LE<C && C<RE) {
    LHSle=LE; LHSre=C-2;
    RHSle=C;  RHSre=RE-1; }
  else {
    if (C<LE) {
      LHSle=0,LHSre=-1;
      RHSle=LE;  RHSre=RE-1; }
    else {//RE<C
      LHSle=LE; LHSre=RE-1;
      RHSle=0,RHSre=-1; }
  }

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

    if (RE<C) {
      //make column(s) C
      for (i=0;i<4;i++) {
	beta(13*g+i,C)=initstore(i,1); }
      for (i=4;i<7;i++) {
	beta(13*g+i,C)=beta(13*g+i-3,C); }


      //make column RE 
      //IMPORTANT Note: I'm commenting this out and setting entries to be 1.0
      //at this locus, because we're conditioning on allele at center.
      //When not used for mapping, this MUST be undone! Code may be correct,
      //though there are numerical difficulties for large tau

      for (i=0;i<7;i++) {
	beta(13*g+i,RE)=1.0; }

      /*
	Q: can you use transstoreC when LE<RE<C?

      //i=0
      sum=beta(13*g  ,C)*transstoreC(0,1)*obsprobsCltRE(14*g  );
      beta(13*g  ,RE)=sum;

      //i=1
      sum=0;
      //add j=0
      sum+=beta(13*g  ,C)*transstoreC(4*1  ,1)*obsprobsCltRE(14*g+1);
      //add j=1,5
      for (indl=0;indl<=1;indl++) {//loop over j's
	j=3*indl+ 1+indl; 
	sum+=beta(13*g+j,C)*transstoreC(4*(1+indl)+(1+indl),1)*
	  obsprobsCltRE(14*g+5*indl+ 4+2*indl);
      }
      beta(13*g+1,RE)=sum;

      //      exit(1);

      //i=2
      sum=0;
      //add j=0
      sum+=beta(13*g  ,C)*transstoreC(4*2  ,1)*obsprobsCltRE(14*g+2);
      //add j=2,4
      for (indl=0;indl<=1;indl++) {//loop over j's
	j=3*indl +2-indl;
	sum+=beta(13*g+j,C)*transstoreC(4*(2-indl)+(2-indl),1)*
	  obsprobsCltRE(14*g+5*indl+ 4+2-2*indl);
      }
      beta(13*g+2,RE)=sum;

      //i=3
      sum=0;
      //add j=0
      sum+=beta(13*g  ,C)*transstoreC(4*3  ,1)*obsprobsCltRE(14*g+3);

      //add j=1:6
      for (indl=0;indl<=1;indl++) {//loop over j's
	//add j=1,5
	j=3*indl+ 1+indl; 
	sum+=beta(13*g+j,C)*transstoreC(4*(3     )+(1+indl),1)*
	  obsprobsCltRE(14*g+5*indl+ 4+1+2*indl);
	//add j=2,4
	j=3*indl +2-indl;
	sum+=beta(13*g+j,C)*transstoreC(4*(3     )+(2-indl),1)*
	  obsprobsCltRE(14*g+5*indl+ 4+3-2*indl);
	//add j=3,6
	j=3*indl+3;
	sum+=beta(13*g+j,C)*transstoreC(4*(3     )+(3     ),1)*
	  obsprobsCltRE(14*g+5*indl+ 4+4);
      } 
      beta(13*g+3,RE)=sum;

      //make denominators
      for (i=0;i<4;i++) {
	sum=initstore(i,3)*obsprobsCden(7*g+i,1);
	if (sum==0) {
	  if (beta(13*g+i,RE)!=0) {
	    errorfs << "ERROR: beta(13*" << g << "+" << i << ",RE = "
		 << "P(A|B) = P(A&B)/P(B) =!0/0" << endl; exit(1); } }
	else beta(13*g+i,RE)/=sum; 
      }
 
      //i=4,5,6
      for (i=4;i<7;i++) {
	beta(13*g+i,RE)=beta(13*g+i -3,RE); }
      */
    }
    else {//C<RE
      //make column(s) RE
      for (i=0;i<4;i++) {
	beta(13*g+i,RE)= initstore(i,2); } 
    }

    //make column(s) RE-1:C or RE-1:LE
    for (t=RHSre;t>=RHSle;t--) {
      if (t==C && (LE<C && C<RE) ) nstatesi=13; else nstatesi=7;
      if (t==RE-1 && C<RE) nstatesj=4; else nstatesj=7;
      for (i=0;i<nstatesi;i++) {
	sum=0;
	if (i>9) tsi=9; //needed to get transstore entry
	else {
	  if (i>6) tsi=6; 
	  else {
	    if (i>3) tsi=3;
	    else tsi=0; } }

	for (j=0;j<nstatesj;j++) {
	  if (j<4) tsj=0; else tsj=3;//needed to get transstore entry

	  if ( (3<i && i<7) || (9<i) ) {
	    if (j==1 || j==4)  adj=1;
	    else {
	      if (j==2 || j==5) adj=-1; 
	      else adj=0; } }
	  else adj=0;
	  
	  sum+=
	    beta(13*g+j,t+1)*transstore(4*(i-tsi)+(j-tsj+adj),t)*
	    obsprobs(7*g+j,t+1); 
	  
	}
	beta(13*g+i,t)=sum; } }

      
    if (LE<C && C<RE) {
      //make column(s) C-1 (see notes for explanation about why this is messy)

      if (C-1==LE) { nstatesi=4; multi=1; }
      else { nstatesi=7; multi=2; }
      if (C+1==RE) { nstatesk=4; multk=1; }
      else { nstatesk=7; multk=2; }


      //make numerators
      //i=0;
      beta(13*g  ,C-1)=beta(13*g  ,C)*transstoreC(0,0)*obsprobsCnum(40*g  ,1);

      //i=1
      sum=0;
      //add j=0
      sum+=beta(13*g  ,C)*transstoreC(4*1  ,0)*obsprobsCnum(40*g+1,1);
      //add j=1,4,8,11
      for (indl=0;indl<=1;indl++) {//loop over j's
	for (indr=0;indr<=1;indr++) {
	  j=3*(2*indl+indr) +1+indl;

	  sum2=initstore(1+indl,1)*multk*
	    ( transstoreC(4*(1+indl)+(1+indl),1)*
	      obsprobsCden(7*g+ (1+indl)+3*indr,0)+
	      transstoreC(4*(1+indl)+(3     ),1)*
	      obsprobsCden(7*g+ (3     )+3*indr,0) );

	  if (sum2!=0) 
	    sum+=beta(13*g+j,C)*transstoreC(4*(1+indl)+(1+indl),0)*multk*
	      ( ( initstore(1+indl,1)* //multk*
		  transstoreC(4*(1+indl)+(1+indl),1)*
		  obsprobsCden(7*g+ (1+indl)+3*indr,0)/sum2 )*
		//transstoreC(4*(1+indl)+(1+indl),1)*
		obsprobsCnum(40*g+18*indl+9*indr +4+  4*indl,1) +
		( initstore(1+indl,1)* //multk*
		  transstoreC(4*(1+indl)+(3     ),1)*
		  obsprobsCden(7*g+ (3     )+3*indr,0)/sum2 )*
		//transstoreC(4*(1+indl)+(3     ),1)*
		obsprobsCnum(40*g+18*indl+9*indr +4+1+4*indl,1) ); 
	} } 
      beta(13*g+1,C-1)=sum;

      //i=2
      sum=0;
      //add j=0
      sum+=beta(13*g  ,C)*transstoreC(4*2  ,0)*obsprobsCnum(40*g+2,1);
      //add j=2,5,7,10
      for (indl=0;indl<=1;indl++) {//loop over j's
	for (indr=0;indr<=1;indr++) {
	  j=3*(2*indl+indr) +2-indl;

	  sum2=initstore(2-indl,1)*multk*
	    ( transstoreC(4*(2-indl)+(2-indl),1)*
	      obsprobsCden(7*g+ (2-indl)+3*indr,0)+
	      transstoreC(4*(2-indl)+(3     ),1)*
	      obsprobsCden(7*g+ (3     )+3*indr,0) );

	  if (sum2!=0)
	    sum+=beta(13*g+j,C)*transstoreC(4*(2-indl)+(2-indl),0)*multk*
	      ( (initstore(2-indl,1)* //multk*
		 transstoreC(4*(2-indl)+(2-indl),1)*
		 obsprobsCden(7*g+ (2-indl)+3*indr,0)/sum2 )*
		//transstoreC(4*(2-indl)+(2-indl),1)*
		obsprobsCnum(40*g+18*indl+9*indr +4+4-4*indl,1) +
		(initstore(2-indl,1)* //multk*
		 transstoreC(4*(2-indl)+(3     ),1)*
		 obsprobsCden(7*g+ (3     )+3*indr,0)/sum2 )*
		//transstoreC(4*(2-indl)+(3     ),1)*
		obsprobsCnum(40*g+18*indl+9*indr +4+5-4*indl,1) ); 
	} } 
      beta(13*g+2,C-1)=sum;

      //i=3
      sum=0;
      //add j=0
      sum+=beta(13*g  ,C)*transstoreC(4*3  ,0)*obsprobsCnum(40*g+3,1);
      //add j=1:12
      for (indl=0;indl<=1;indl++) {
	for (indr=0;indr<=1;indr++) {
	  //add j=1,4,8,11
	  j=3*(2*indl+indr) +1+indl;

	  sum2=initstore(1+indl,1)*multk*
	    ( transstoreC(4*(1+indl)+(1+indl),1)*
	      obsprobsCden(7*g+ (1+indl)+3*indr,0)+
	      transstoreC(4*(1+indl)+(3     ),1)*
	      obsprobsCden(7*g+ (3     )+3*indr,0) );

	  if (sum2!=0) 
	    sum+=beta(13*g+j,C)*transstoreC(4*(3     )+(1+indl),0)*multk*
	      ( ( initstore(1+indl,1)* //multk*
		  transstoreC(4*(1+indl)+(1+indl),1)*
		  obsprobsCden(7*g+ (1+indl)+3*indr,0)/sum2 )*
		//transstoreC(4*(1+indl)+(1+indl),1)*
		obsprobsCnum(40*g+18*indl+9*indr +4+2+4*indl,1) +
		( initstore(1+indl,1)* //multk*
		  transstoreC(4*(1+indl)+(3     ),1)*
		  obsprobsCden(7*g+ (3     )+3*indr,0)/sum2 )*
		//transstoreC(4*(1+indl)+(3     ),1)*
		obsprobsCnum(40*g+18*indl+9*indr +4+3+4*indl,1) ); 
	  //add j=2,5,7,10
	  j=3*(2*indl+indr) +2-indl;

	  sum2=initstore(2-indl,1)*multk*
	    ( transstoreC(4*(2-indl)+(2-indl),1)*
	      obsprobsCden(7*g+ (2-indl)+3*indr,0)+
	      transstoreC(4*(2-indl)+(3     ),1)*
	      obsprobsCden(7*g+ (3     )+3*indr,0) );

	  if (sum2!=0)
	    sum+=beta(13*g+j,C)*transstoreC(4*(3     )+(2-indl),0)*multk*
	      ( (initstore(2-indl,1)* //multk*
		 transstoreC(4*(2-indl)+(2-indl),1)*
		 obsprobsCden(7*g+ (2-indl)+3*indr,0)/sum2 )*
		//transstoreC(4*(2-indl)+(2-indl),1)*
		obsprobsCnum(40*g+18*indl+9*indr +4+6-4*indl,1) +
		(initstore(2-indl,1)* //multk*
		 transstoreC(4*(2-indl)+(3     ),1)*
		 obsprobsCden(7*g+ (3     )+3*indr,0)/sum2 )*
		//transstoreC(4*(2-indl)+(3     ),1)*
		obsprobsCnum(40*g+18*indl+9*indr +4+7-4*indl,1) ); 

	  //add j=3,6,9,12
	  j=3*(2*indl+indr) +3;

	  sum2=initstore(3,1)*multk*
	    ( transstoreC(4*3+3,1)*obsprobsCden(7*g+ 3+3*indr,0) );

	  if (sum2!=0)
	    sum+=beta(13*g+j,C)*transstoreC(4*(3     )+(3     ),0)*multk*
	      ( (initstore(3,1)* //multk*
		 transstoreC(4*3+3,1)*obsprobsCden(7*g+ 3+3*indr,0)/sum2 )*
		//transstoreC(4*(3     )+(3     ),1)*
		obsprobsCnum(40*g+18*indl+9*indr+ 4+8       ,1) ); 
	} } 
      beta(13*g+3,C-1)=sum;


    
      //make denominators
      for (i=0;i<4;i++) {
	sum=initstore(i,3)*obsprobsCden(7*g+i,1);

	if (sum==0) {
	  if (beta(13*g+i,C-1)!=0) {
	    errorfs << "ERROR: beta(13*" << g << "+" << i << ",C-1) = "
		 << "P(A|B) = P(A&B)/P(B) =!0/0" << endl; exit(1); } }
	else beta(13*g+i,C-1)/=sum; 
      }
	
      //i=4,5,6
      if (LE!=C-1) { 
	for (i=4;i<7;i++) {
	  beta(13*g+i,C-1)=beta(13*g+i -3,C-1); } }

    }
    

    //make column(s) C-2:LE or RE-1:LE
    for (t=LHSre;t>=LHSle;t--) {
      if (t==LE  && LE<C) nstatesi=4; else nstatesi=7;
      for (i=0;i<nstatesi;i++) {
	sum=0;
	if (i<4) tsi=0; else tsi=3;//needed to get transstore entry
	for (j=0;j<7;j++) {
	  if (j<4) tsj=0; else tsj=3;

	  if (j>3) {
	    if (i==1 || i==4) adj=1;
	    else {
	      if (i==2 || i==5) adj=-1; 
	      else adj=0; } }
	  else adj=0;
	    
	  sum+=
	    beta(13*g+j,t+1)*transstore(4*(i-tsi+adj)+(j-tsj),t)*
	    obsprobs(7*g+j,t+1); 

	}
	beta(13*g+i,t)=sum; 


      }
    }
  
    if (C<LE) {
      //make column(s) C
      nstatesi=7;
      nstatesj=7;
      for (i=0;i<nstatesi;i++) {
	sum=0;
	if (i>9) tsi=9; //needed to get transstore entry
	else {
	  if (i>6) tsi=6; 
	  else {
	    if (i>3) tsi=3;
	    else tsi=0; } }

	for (j=0;j<nstatesj;j++) {
	  if (j<4) tsj=0; else tsj=3;//needed to get transstore entry

	  if ( (3<i && i<7) || (9<i) ) {
	    if (j==1 || j==4)  adj=1;
	    else {
	      if (j==2 || j==5) adj=-1; 
	      else adj=0; } }
	  else adj=0;
	  
	  sum+=
	    beta(13*g+j,LE)*transstore(4*(i-tsi)+(j-tsj+adj),C)*
	    obsprobs(7*g+j,LE); 
	}
	beta(13*g+i,C)=sum; }
      
    }
  }  

  for (i=0;i<beta.dimx();i++) {
    if (C<LE || RE<C) {
      if ( !(-0.0001<beta(i,C) && beta(i,C)<1.0001) ) {
	errorfs << "ERROR: beta(" << i << "," << C << ") = " 
	     << beta(i,C) << endl;
	exit(1);
      }
    }
    for (j=LE;j<=RE;j++)
      if ( !(-0.0001<beta(i,j) && beta(i,j)<1.0001) ) {
	errorfs << "ERROR: beta(" << i << "," << j << ") = " 
	     << beta(i,j) << endl;
	exit(1);
      }
  }
}

void makegammagen(dgrid &gamma, const dgrid &alpha, const dgrid &beta,
		  const int &C, const int &LE, const int &RE)
{
  //makes gamma object, for genotype data for affecteds

  int i,g,t,nstates;
  double denom, sum;

  const int numgen=gamma.dimx()/13;

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

    //calculate likelihood for genotype
    if (C<LE || RE<C) nstates=7; else nstates=13;
    denom=0;
    for (i=0;i<nstates;i++) { denom+=alpha(13*g+i,C)*beta(13*g+i,C); }

    if (C<LE || RE<C) {
      for (i=0;i<7;i++) { 
	gamma(13*g+i,C)=( alpha(13*g+i,C)*beta(13*g+i,C) )/denom; } }

    for (t=LE;t<=RE;t++) {
      if (t==C) nstates=13;
      else {
	if ( (t==LE && LE<C) || (t==RE && C<RE) ) nstates=4;
	else nstates=7; }
      for (i=0;i<nstates;i++) { 
	gamma(13*g+i,t)=alpha(13*g+i,t)*beta(13*g+i,t)/denom; }
  
      //check that gamma's sum to 1 at each locus
      sum=0;
      for (i=0;i<nstates;i++) sum+=gamma(13*g+i,t);

      if ( !(0.999<sum && sum<1.001) ) {
	errorfs << "ERROR: gamma's do not sum to 1.0 (" << sum <<") for gen. "  
	     << g << "+1 and locus " << t << "+1" << endl; 
	exit(1);

      }
    }
  }

  for (g=0;g<gamma.dimx();g++) {
    if (C<LE || RE<C) 
      assert(-0.00001<gamma(g,C) && gamma(g,C)<1.00001);
    for (t=LE;t<=RE;t++)
      assert(-0.00001<gamma(g,t) && gamma(g,t)<1.00001); }

}

void makecstargen(dvector &cstar, const dgrid &gamma, 
		  const int &C, const int &LE, const int &RE)
{
  //makes the cstar, expected counts of haplotypes still sharing from
  //ancestral haplotype at each marker, for genotype data

  int n,i,g,t;
  const int numgen=gamma.dimx()/13;

  if (C<LE || RE<C) {
    cstar(C)=0;
    for (g=0;g<numgen;g++) {
      cstar(C)+=2*gamma(13*g,C);
      for (i=0;i<2;i++) {
	cstar(C)+=gamma(13*g+1 +3*i,C)+gamma(13*g+2 +3*i,C); } } }

  for (t=LE;t<=RE;t++) {
    if (t==C) n=4;
      else {
	if ( (t==LE && LE<C) || (t==RE && C<RE) ) n=1;
	else n=2; }
    cstar(t)=0;
    for (g=0;g<numgen;g++) {
      cstar(t)+=2*gamma(13*g,t);
      for (i=0;i<n;i++) {
	cstar(t)+=gamma(13*g+1 +3*i,t)+gamma(13*g+2 +3*i,t); } } }


  //assertions
  for (t=LE;t<=RE;t++) {
    assert(cstar(t)<=cstar(C)); }

} 

void makecstarmutgen(dvector &cstarmut, const dvector &cstar,
		     const dgrid &gamma,
		     const igrid &data, const ivector &anchap,
		     const dvector &mutmatch,
		     const int &C, const int &LE, const int &RE)
{
  //makes subset of cstars which do not match ancestral hap., only
  //due to mutation, for genotype data

  int n,g,l,t,hap;
  const int numgen=gamma.dimx()/13;
  double sum;

  //assumes no mutation allowed at center
  cstarmut(C)=0;  

  for (t=LE;t<=RE;t++) {
    if (t!=C) {
      cstarmut(t)=0;
      if ( (t==LE && LE<C) || (t==RE && C<RE) ) n=1; else n=2;
      for (g=0;g<numgen;g++) {
	
	for (hap=0;hap<=1;hap++) {
	  sum=0;
	  if (data(2*g+hap,t)!=anchap(t)) {
	    if (data(2*g+hap,t)==0) { 
	      sum=gamma(13*g  ,t);
	      for (l=0;l<n;l++) {
		sum+=gamma(13*g+ 3*l +1+hap  ,t); }
	      cstarmut(t)+=(1-mutmatch[t])*sum; }
	    else {
	      sum=gamma(13*g  ,t);
	      for (l=0;l<n;l++) {
		sum+=gamma(13*g+ 3*l +1+hap  ,t); }
	      cstarmut(t)+=sum; } }
	  
	}
	
      }
    }
  }
}

void makeloglikgen(dvector &loglikvec, const dgrid &alpha, 
		   const dgrid &initstore,
		   const int &C, const int &RE)
{
  //makes vector of loglikelihoods for each genotype from alpha and
  //unconditional distribution of chain at one locus

  int i,g,locus,nstates,tsi;
  double lik;
  const int numgen=alpha.dimx()/13;
  if (C<RE) { locus=RE; nstates=4; } 
  else { locus=C; nstates=7; }
  
  for (g=0;g<numgen;g++) {
    loglikvec[g]=0;
    lik=0;
    for (i=0;i<nstates;i++) {
      if (i<4) tsi=0; else tsi=3;
      lik+=alpha(13*g+i,locus)*initstore(i-tsi,2); 
    }
    assert(0<lik && lik<=1.0001);
    loglikvec[g]=log(lik);
  }
}