#include <math.h>
#include <stdio.h>

/* These versions of the routines allow for recombination between the
 loci.
 The function get.2locexp.recomb() takes as input the 2-locus
 sperm-typing data counts (length 16) in order
 -,b,a,ab,B,Bb,Ba,Bab,A,Ab,Aa,Aab,AB,ABb,ABa,ABab
 and the current parameter values (length 11) in order
 s (segregation probability for allele 1, assoc. with AB in parent)
 p=(p1,p2) probabilities of 1 and 2 sperm
 deposited with p0 assumed = 1-p1-p2, am vector of 4 amplification
 probabilities for alleles A,B,a, and b; co vector of 4 contamination
 probabilities for alleles A,B,a, and b.  Also takes as input the
 multilocus recombination probabilities:  Pr(no rec., i.e. 1AB or 2ab),
 Pr(rec btw disease and marker haplotype, i.e. 1ab or 2AB), Pr(1Ab or
 2aB), Pr(1aB or 2Ab), where parental haplotypes are 1AB and 2ab.

 The output of get.2locexp.recomb() is the expected values of the
 following
 quantities (length 21), conditional on the data and the current
 parameter values:
 number of segregations of allele 1, associated in parent with AB,
 number of segregations of allele 2, associated in parent with ab,
 number of times 1 sperm is deposited, number of times 2 sperm are
 deposited, number of times 0 sperm are deposited,
 number of times
 allele A is amplified, number of times allele A is present but not
 amplified, number of times allele B is amplified, number of times B is
 present but not amplified, similarly for a and b, number of times
allele
 A  is a contaminant, number of times allele A is not a contaminant,
 number of times allele B is/is not a contaminant, a is/is not a
 contaminant, b is/is not a contaminant.

 The intention is that this function will be run on each two-locus data
 set, the function get.1locexp.recomb() will be run on each one-locus
data
 set, and then the results of all these runs will be used in the next
 step of the E-M algorithm to estimate the parameters (under selected
 constraints) by maximum likelihood.*/

/*x is 16x4, ith row of x2is i-1 base 2. Represents observed phenotype*/
/* y is 8x2, z is 64x4, code is 4x4, prpheno0 is 16, prpheno1 is 16x8, 
prpheno2 is 16x64, mpp is 4, denom is 16, tmp is 16, tmpt is 16, ltmp is
16 */
/* The values of x, y, z, and code are unchanging and are prespecified
outside this subroutine.  Can think of the prpheno's, mpp, denom, tmp,
tmpt, ltmp as temporary storage passed to this routine.*/

void get2locexprecomb(double *data,double *s, double *p, double *am,
double *co, double *mp, double *output, int **x, int **y, int **z, 
int **code, double *prpheno0, double **prpheno1, double **prpheno2, 
double *mpp,double *denom, double *tmp, double *tmpt, double *ltmp)
{

int i,j,k;

/*prpheno0[i] is prob. of (A,B,a,b)^(x[i,]) and no sperm present*/

/*printf("prpheno0\n");*/
for(i=0;i<16;i++){
	prpheno0[i]=(1.-p[0]-p[1]);
	for(j=0;j<4;j++) prpheno0[i]*=(pow(co[j],(double)x[i][j])*
	pow(1.-co[j],1.-(double)x[i][j]));
	/*printf("%e ",prpheno0[i]);*/
}

mpp[0]=mp[1]; mpp[1]=mp[0]; mpp[2]=mp[3]; mpp[3]=mp[2];
/*printf("\nmpp\n");
for(i=0;i<4;i++) printf("%lf ",mpp[i]);*/

/* Now find prob. of each phenotype event and each genotype involving
1 sperm present.  prpheno1 is 16x8 matrix, ith row corresponds to
phenotype x[i,].  kth column corresponds to genotype y[k,], where 1st
digit of y is 1=mutant allele or 2=normal, 2nd digit is 1=AB, 2=ab,
3=Ab, 4=aB (given by code matrix). */

/*printf("\nprpheno1\n");*/
for(k=0;k<4;k++) for(i=0;i<16;i++){
	prpheno1[i][k]=p[0]**s*mp[k];
	prpheno1[i][4+k]=p[0]*(1.-*s)*mpp[k];
}
for(k=0;k<8;k++) for(i=0;i<16;i++) for(j=0;j<4;j++)
prpheno1[i][k]*=(pow(1.-pow(1.-am[j],(double)code[y[k][1]][j])*
(1.-co[j]),(double)x[i][j])*pow(pow(1.-am[j],(double)code[y[k][1]][j])*
(1.-co[j]),1.-(double)x[i][j]));

/*for(i=0;i<16;i++){
	for(k=0;k<8;k++) printf("%e ",prpheno1[i][k]);
	printf("\n");
}*/

/* Now find prob. of each phenotype event and each genotype involving
2 sperm present.  prpheno2 is 16x64 matrix, ith row corresponds to
phenotype x[i,].  kth column corresponds to genotype z[k,], where 1st
digit of z is 1st sperm: 1=mutant allele or 2=normal, 2nd digit is
2nd sperm: same, 3rd digit is 1st sperm: 1=AB, 2=ab,
3=Ab, 4=aB (given by code matrix), 4th digit is 2nd sperm: same */

for(k=0;k<64;k++) for(i=0;i<16;i++){
	prpheno2[i][k]=p[1]*pow(*s, 2.-(double)(z[k][0]+z[k][1]))*
	pow(1.-*s,(double)(z[k][0]+z[k][1]))*pow(mp[z[k][2]],
	1.-(double)z[k][0])*pow(mpp[z[k][2]],(double)z[k][0])*
	pow(mp[z[k][3]],1.-(double)z[k][1])*
	pow(mpp[z[k][3]],(double)z[k][1]);

	for(j=0;j<4;j++) prpheno2[i][k]*=(pow(1.-pow(1.-am[j],
	(double)(code[z[k][2]][j]+code[z[k][3]][j]))*(1.-co[j]),
	(double)x[i][j])*pow(pow(1.-am[j],(double)(code[z[k][2]][j]+
	code[z[k][3]][j]))*(1.-co[j]),1.-(double)x[i][j]));
}

/*printf("\nprpheno2\n");
for(i=0;i<16;i++){
	for(j=0;j<64;j++) printf("%e ",prpheno2[i][j]);
	printf("\n");
}*/

/*printf("denom\n");*/
for(i=0;i<16;i++){
	denom[i]=prpheno0[i];
	for(j=0;j<8;j++) denom[i]+=prpheno1[i][j];
	for(j=0;j<64;j++) denom[i]+=prpheno2[i][j];
	/*printf("%e ",denom[i]);*/
}

/* number of segregations of 1 (mutant), of 2 (normal) */

output[0]=0.;
for(i=0;i<16;i++){
	for(j=0;j<4;j++) if(denom[i]>0.) output[0]+=(data[i]/denom[i]*
	prpheno1[i][j]);
	for(j=0;j<16;j++) if(denom[i]>0.) output[0]+=(data[i]/denom[i]*2.*
	prpheno2[i][j]);
	for(j=16;j<48;j++) if(denom[i]>0.) output[0]+=(data[i]/denom[i]*
	prpheno2[i][j]);
}

output[1]=0.;
for(i=0;i<16;i++){
        for(j=4;j<8;j++) if(denom[i]>0.) output[1]+=(data[i]/denom[i]*
        prpheno1[i][j]);
        for(j=48;j<64;j++) if(denom[i]>0.) output[1]+=(data[i]/denom[i]*2.*
        prpheno2[i][j]);
        for(j=16;j<48;j++) if(denom[i]>0.) output[1]+=(data[i]/denom[i]*
        prpheno2[i][j]);
}

/* number of times 1 sperm is deposited, 2 sperm, 0 sperm */

output[2]=0.;
for(i=0;i<16;i++) for(j=0;j<8;j++) if(denom[i]>0.) 
output[2]+=(data[i]/denom[i]*
prpheno1[i][j]);

output[3]=0.;
for(i=0;i<16;i++) for(j=0;j<64;j++) if(denom[i]>0.) 
output[3]+=(data[i]/denom[i]*
prpheno2[i][j]);

output[4]=-output[2]-output[3];
for(i=0;i<16;i++) output[4]+=data[i];

/* number of times each allele is amplified/ is present but is not
 amplified */

for(j=0;j<4;j++){
	for(i=0;i<16;i++){
		tmp[i]=0.;
		tmpt[i]=0.;
		ltmp[i]=0.;
	}
	for(k=0;k<8;k++){
		for(i=0;i<16;i++){
			tmp[i]+=((double)code[y[k][1]][j]*
			prpheno1[i][k]*am[j]/(am[j]+(1.-am[j])*
			co[j]));
			tmpt[i]+=((double)code[y[k][1]][j]*
                        prpheno1[i][k]);
		}
	}
	for(k=0;k<64;k++){
		for(i=0;i<16;i++){
			if(pow(1.-pow(1.-am[j],(double)(
			code[z[k][2]][j]+code[z[k][3]][j]))*
			(1.-co[j]),(double)x[i][j])*pow(pow(
			1.-am[j],(double)(code[z[k][2]][j]+
			code[z[k][3]][j]))*(1.-co[j]),1.-
			(double)x[i][j])==0.) ltmp[i]=1.;

			tmp[i]+=((1.-ltmp[i])*(double)(
			code[z[k][2]][j]+code[z[k][3]][j])*
			prpheno2[i][k]*am[j]/(
			pow(1.-pow(1.-am[j],(double)(
                        code[z[k][2]][j]+code[z[k][3]][j]))*
                        (1.-co[j]),(double)x[i][j])*pow(pow(
                        1.-am[j],(double)(code[z[k][2]][j]+
                        code[z[k][3]][j]))*(1.-co[j]),1.-
                        (double)x[i][j])+ltmp[i]));

			tmpt[i]+=((1.-ltmp[i])*(double)(
                        code[z[k][2]][j]+code[z[k][3]][j])*
                        prpheno2[i][k]);
		}
	}
	output[3+2*(j+1)]=0.;
	output[4+2*(j+1)]=0.;
	for(i=0;i<16;i++){
		if(x[i][j]!=1) tmp[i]=0.;
		if (denom[i]>0.)
		output[3+2*(j+1)]+=(data[i]/denom[i]*tmp[i]);
		if(denom[i]>0.)
		output[4+2*(j+1)]+=(data[i]/denom[i]*(tmpt[i]-tmp[i]));
	}
}

/* number of times each allele is/is not a contaminant */

for(j=0;j<4;j++){

	for(i=0;i<16;i++){
		tmp[i]=prpheno0[i];
		ltmp[i]=0.;
	}

	for(k=0;k<8;k++){
		for(i=0;i<16;i++) tmp[i]+=(prpheno1[i][k]*pow((1.-
		am[j])*co[j]/(am[j]+(1.-am[j])*co[j]),
		(double)code[y[k][1]][j]));
	}

	for(k=0;k<64;k++){
		for(i=0;i<16;i++){
			if(1.-pow(1.-am[j],(double)(
                        code[z[k][2]][j]+code[z[k][3]][j]))*
                        (1.-co[j])==0.) ltmp[i]=1.;

			tmp[i]+=((1.-ltmp[i])*prpheno2[i][k]*pow(1.-
			am[j],(double)(code[z[k][2]][j]+
			code[z[k][3]][j]))*co[j]/(1.-pow(1.-am[j],
			(double)(code[z[k][2]][j]+code[z[k][3]][j]))*
			(1.-co[j])+ltmp[i]));
		}
	}

	output[11+2*(j+1)]=0.;
	output[12+2*(j+1)]=0.;
	for(i=0;i<16;i++){
		if(x[i][j]!=1) tmp[i]=0.;
		if(denom[i]>0.)
		output[11+2*(j+1)]+=(data[i]/denom[i]*tmp[i]);
                output[12+2*(j+1)]+=(data[i]*
		(double)(1-x[i][j]));
        }
}
}

void makexyzcode(int **x,int **y,int **z,int **code)

/* x is 16x4, y is 8x2, z is 64x4, code is 4x4 */

{
int i,j;

for(i=0;i<16;i++){
	x[i][3]=(int)fmod((double)i,2.);
	x[i][2]=(int)fmod((double)((i-x[i][3])/2),2.);
	x[i][1]=(int)fmod((double)(((i-x[i][3])/2-x[i][2])/2),2.);
	x[i][0]=(int)fmod((double)((((i-x[i][3])/2-x[i][2])/2-x[i][1])/
	2),2.);
}

for(j=0;j<4;j++){
	y[j][0]=0;
	y[j+4][0]=1;
	y[j][1]=j;
	y[j+4][1]=j;
}

code[0][0]=1;
code[0][1]=1;
code[0][2]=0;
code[0][3]=0;
code[1][0]=0;
code[1][1]=0;
code[1][2]=1;
code[1][3]=1;
code[2][0]=1;
code[2][1]=0;
code[2][2]=0;
code[2][3]=1;
code[3][0]=0;
code[3][1]=1;
code[3][2]=1;
code[3][3]=0;
	
for(i=0;i<32;i++){
	z[i][0]=0;
	z[i+32][0]=1;
}
for(i=0;i<16;i++){
	z[i][1]=0;
	z[i+16][1]=1;
	z[i+32][1]=0;
	z[i+48][1]=1;
}
for(i=0;i<4;i++){
	for(j=0;j<4;j++){
		z[i+j*16][2]=0;
		z[i+4+j*16][2]=1;
		z[i+8+j*16][2]=2;
		z[i+12+j*16][2]=3;
	}
}
for(i=0;i<4;i++) for(j=0;j<16;j++) z[j*4+i][3]=i;

}

void makeonexyzcode(int **x1,int **y1,int **z1,int **code1)

/* x1 is 4x2, y1 is 4x2, z1 is 16x4 code1 is 2x2 */

{
int i,j;

for(i=0;i<4;i++){
	x1[i][1]=(int)fmod((double)i,2.);
	x1[i][0]=(int)fmod((double)((i-x1[i][1])/2),2.);
}

for(i=0;i<2;i++){
	y1[i][0]=0;
	y1[i+2][0]=1;
	y1[i][1]=i;
	y1[i+2][1]=i;
}

code1[0][0]=1;
code1[0][1]=0;
code1[1][0]=0;
code1[1][1]=1;

for(i=0;i<8;i++){
	z1[i][0]=0;
	z1[i+8][0]=1;
}
for(i=0;i<4;i++){
	z1[i][1]=0;
	z1[i+4][1]=1;
	z1[i+8][1]=0;
	z1[i+12][1]=1;
}
for(i=0;i<2;i++){
	for(j=0;j<4;j++){
		z1[i+j*4][2]=0;
		z1[i+2+j*4][2]=1;
	}
}
for(i=0;i<2;i++) for(j=0;j<8;j++) z1[j*2+i][3]=i;

}
	

/* The function get.1locexp.recomb() takes as input the 1-locus sperm
typing
 data counts (length 4) in order -,a,A,Aa
 and the current parameter values (length 7) in order s (segregation
 probability for allele 1, assoc. with AB in parent) p=(p1,p2)
 probabilities of 1 and 2 sperm
 deposited with p0 assumed = 1-p1-p2, am vector of 2 amplification
 probabilities for alleles A and a; co vector of 2 contamination
 probabilities for alleles A and a.  Also takes as input the
 recombination fraction between the gene and the marker.

 The output of get.1locexp.recomb() is the expected values of the
following
 quantities (length 13), condional on the data and the current
 parameter values:
 number of  segregations of allele 1, associated in parent with A,
 number of segregations of allele 2, associated in parent with a,
 number of times 1 sperm is deposited, number of times 2 sperm are
 deposited, number of times 0 sperm are deposited,
 number of times
 allele A is amplified, number of times allele A is present but not
 amplified, similarly for a, number of times allele
 A  is a contaminant, number of times allele A is not a contaminant,
 number of times a is/is not a contaminant

 For the purpose, see under get.2locexp.recomb() above.

x is 4x2, ith row is i-1 base 2. y is 4x2, z is 16x4, code is 2x2,
prpheno0 is 4, prpheno1 is 4x4, prpheno2 is 4x16, denom is 4,
tmp is 4, tmpt is 4, ltmp is 4*/

void get1locexprecomb(double *data,double *s,double *p,double *am,double
*co,double *rf,double *output,int **x,int **y,int **z, int **code,
double *prpheno0, double **prpheno1,double **prpheno2,double *denom,
double *tmp,double *tmpt,double *ltmp)
{

int i,j,k;

/* prpheno0[i] is prob of (A,a)^x[i,] and no sperm present */

/*printf("prpheno0\n");*/
for(i=0;i<4;i++){
        prpheno0[i]=(1.-p[0]-p[1]);
        for(j=0;j<2;j++) prpheno0[i]*=(pow(co[j],(double)x[i][j])*
        pow(1.-co[j],1.-(double)x[i][j]));
        /*printf("%e ",prpheno0[i]);
	fflush(stdout);*/
}

/* Now find prob. of each phenotype event and each genotype involving
 1 sperm present.  prpheno1 is 4x4 matrix, ith row corresponds to
 phenotype x[i,].  kth column corresponds to genotype y[k,], where 1st
 digit of y is 1=mutant allele or 2=normal, 2nd digit is 1=A, 2=a
 (given by code matrix).*/

/*printf("\nprpheno1\n");*/
for(i=0;i<4;i++){
	prpheno1[i][0]=p[0]**s*(1.-*rf);
	prpheno1[i][1]=p[0]**s**rf;
	prpheno1[i][2]=p[0]*(1.-*s)**rf;
	prpheno1[i][3]=p[0]*(1.-*s)*(1.-*rf);
}
for(k=0;k<4;k++) for(i=0;i<4;i++) for(j=0;j<2;j++)
prpheno1[i][k]*=(pow(1.-pow(1.-am[j],(double)code[y[k][1]][j])*
(1.-co[j]),(double)x[i][j])*pow(pow(1.-am[j],(double)code[y[k][1]][j])*
(1.-co[j]),1.-(double)x[i][j]));

/*for(i=0;i<4;i++){
        for(k=0;k<4;k++) printf("%e ",prpheno1[i][k]);
        printf("\n");
	fflush(stdout);
}*/

/* Now find prob. of each phenotype event and each genotype involving
 2 sperm present.  prpheno2 is 4x16 matrix, ith row corresponds to
 phenotype x[i,].  kth column corresponds to genotype z[k,], where 1st
 digit of z is 1st sperm: 1=mutant allele or 2=normal, 2nd digit is
 2nd sperm: same, 3rd digit is 1st sperm: 1=A, 2=a
 (given by code matrix), 4th digit is 2nd sperm: same */

for(k=0;k<16;k++) for(i=0;i<4;i++){
        prpheno2[i][k]=p[1]*pow(*s, 2.-(double)(z[k][0]+z[k][1]))*
        pow(1.-*s,(double)(z[k][0]+z[k][1]))*pow(*rf,fmod(fabs((double)
	(z[k][0]-z[k][2])),2.)+fmod(fabs((double)(z[k][1]-z[k][3])),2.))*
	pow(1.-*rf,fmod(fabs((double)(z[k][0]-z[k][2]+1)),2.)+
	fmod(fabs((double)(z[k][1]-z[k][3]+1)),2.));

        for(j=0;j<2;j++) prpheno2[i][k]*=(pow(1.-pow(1.-am[j],
        (double)(code[z[k][2]][j]+code[z[k][3]][j]))*(1.-co[j]),
        (double)x[i][j])*pow(pow(1.-am[j],(double)(code[z[k][2]][j]+
        code[z[k][3]][j]))*(1.-co[j]),1.-(double)x[i][j]));
}

/*printf("\nprpheno2\n");
for(i=0;i<4;i++){
        for(j=0;j<16;j++) printf("%e ",prpheno2[i][j]);
        printf("\n");
	fflush(stdout);
}*/

/*printf("denom\n");*/
for(i=0;i<4;i++){
        denom[i]=prpheno0[i];
        for(j=0;j<4;j++) denom[i]+=prpheno1[i][j];
        for(j=0;j<16;j++) denom[i]+=prpheno2[i][j];
        /*printf("%e ",denom[i]);
	fflush(stdout);*/
}

/* number of segregations of 1 (mutant), of 2 (normal) */

output[0]=0.;
for(i=0;i<4;i++){
        for(j=0;j<2;j++) if(denom[i]>0.) 
	output[0]+=(data[i]/denom[i]*
        prpheno1[i][j]);
        for(j=0;j<4;j++) if(denom[i]>0.) 
	output[0]+=(data[i]/denom[i]*2.*
        prpheno2[i][j]);
        for(j=4;j<12;j++) if(denom[i]>0.)
	output[0]+=(data[i]/denom[i]*
        prpheno2[i][j]);
}

/*printf("output0 %lf\n",output[0]);
fflush(stdout);*/

output[1]=0.;
for(i=0;i<4;i++){
        for(j=2;j<4;j++) if(denom[i]>0.)
	output[1]+=(data[i]/denom[i]*
        prpheno1[i][j]);
        for(j=12;j<16;j++) if(denom[i]>0.)
	output[1]+=(data[i]/denom[i]*2.*
        prpheno2[i][j]);
        for(j=4;j<12;j++) if(denom[i]>0.)
	output[1]+=(data[i]/denom[i]*
        prpheno2[i][j]);
}

/*printf("output1 %lf\n",output[1]);
fflush(stdout);*/

/* number of times 1 sperm is deposited, 2 sperm, 0 sperm */

output[2]=0.;
for(i=0;i<4;i++) for(j=0;j<4;j++) if(denom[i]>0.)
output[2]+=(data[i]/denom[i]*
prpheno1[i][j]);

/*printf("output2 %lf\n",output[2]);
fflush(stdout);*/

output[3]=0.;
for(i=0;i<4;i++) for(j=0;j<16;j++) if(denom[i]>0.)
output[3]+=(data[i]/denom[i]*
prpheno2[i][j]);

/*printf("output3 %lf\n",output[3]);
fflush(stdout);*/

output[4]=-output[2]-output[3];
for(i=0;i<4;i++) output[4]+=data[i];

/*printf("output4 %lf\n",output[4]);
fflush(stdout);*/

/* number of times each allele is amplified/ is present but is not
 amplified */

for(j=0;j<2;j++){
        for(i=0;i<4;i++){
                tmp[i]=0.;
                tmpt[i]=0.;
		ltmp[i]=0.;
        }
        for(k=0;k<4;k++){
		if(am[j]+(1.-am[j])*co[j]==0.) *ltmp=1.; else *ltmp=0.;
                for(i=0;i<4;i++){
                        tmp[i]+=((1.-*ltmp)*(double)code[y[k][1]][j]*
                        prpheno1[i][k]*am[j]/(am[j]+(1.-am[j])*
                        co[j]+*ltmp));
                        tmpt[i]+=((double)code[y[k][1]][j]*
                        prpheno1[i][k]);
                }
        }
	ltmp[0]=0.;
        for(k=0;k<16;k++){
                for(i=0;i<4;i++){
                        if(pow(1.-pow(1.-am[j],(double)(
                        code[z[k][2]][j]+code[z[k][3]][j]))*
                        (1.-co[j]),(double)x[i][j])*pow(pow(
                        1.-am[j],(double)(code[z[k][2]][j]+
                        code[z[k][3]][j]))*(1.-co[j]),1.-
                        (double)x[i][j])==0.) ltmp[i]=1.; 
			else ltmp[i]=0.;

                        tmp[i]+=((1.-ltmp[i])*(double)(
                        code[z[k][2]][j]+code[z[k][3]][j])*
                        prpheno2[i][k]*am[j]/(
                        pow(1.-pow(1.-am[j],(double)(
                        code[z[k][2]][j]+code[z[k][3]][j]))*
                        (1.-co[j]),(double)x[i][j])*pow(pow(
                        1.-am[j],(double)(code[z[k][2]][j]+
                        code[z[k][3]][j]))*(1.-co[j]),1.-
                        (double)x[i][j])+ltmp[i]));

                        tmpt[i]+=((double)(
                        code[z[k][2]][j]+code[z[k][3]][j])*
                        prpheno2[i][k]);
                }
        }
        output[3+2*(j+1)]=0.;
        output[4+2*(j+1)]=0.;
        for(i=0;i<4;i++){
                if(x[i][j]!=1) tmp[i]=0.;
		if(denom[i]>0.)
                output[3+2*(j+1)]+=(data[i]/denom[i]*tmp[i]);
		if(denom[i]>0.)
                output[4+2*(j+1)]+=(data[i]/denom[i]*(tmpt[i]-tmp[i]));
        }
}

/*printf("output5 %lf\n",output[5]);
printf("output6 %lf\n",output[6]);
printf("output7 %lf\n",output[7]);
printf("output8 %lf\n",output[8]);
fflush(stdout);*/

/* number of times each allele is/is not a contaminant */

for(j=0;j<2;j++){

        for(i=0;i<4;i++){
                tmp[i]=prpheno0[i];
                ltmp[i]=0.;
        }

        for(k=0;k<4;k++){
		if(pow(am[j]+(1.-am[j])*co[j],
		(double)code[y[k][1]][j])==0.) *ltmp=1.; 
		else *ltmp=0.;

                for(i=0;i<4;i++) tmp[i]+=((1.-*ltmp)*
		prpheno1[i][k]*pow((1.-
                am[j])*co[j]/(am[j]+(1.-am[j])*co[j]+*ltmp),
                (double)code[y[k][1]][j]));
        }

	*ltmp=0.;
        for(k=0;k<16;k++){
                for(i=0;i<4;i++){
                        if(1.-pow(1.-am[j],(double)(
                        code[z[k][2]][j]+code[z[k][3]][j]))*
                        (1.-co[j])==0.) ltmp[i]=1.; else ltmp[i]=0.;

                        tmp[i]+=((1.-ltmp[i])*prpheno2[i][k]*pow(1.-
                        am[j],(double)(code[z[k][2]][j]+
                        code[z[k][3]][j]))*co[j]/(1.-pow(1.-am[j],
                        (double)(code[z[k][2]][j]+code[z[k][3]][j]))*
                        (1.-co[j])+ltmp[i]));
                }
        }

        output[7+2*(j+1)]=0.;
        output[8+2*(j+1)]=0.;
        for(i=0;i<4;i++){
                if(x[i][j]!=1) tmp[i]=0.;
		if(denom[i]>0.)
                output[7+2*(j+1)]+=(data[i]/denom[i]*tmp[i]);
                output[8+2*(j+1)]+=(data[i]*
                (double)(1-x[i][j]));
        }
}
}

/* The function get.2loglike.recomb() takes as input the 2-locus
sperm-typing
 data counts (length 16) in order given above and the current
 parameter values (length 11) in order given above. Also takes as input
 the
 multilocus recombination probabilities:  Pr(no rec., i.e. 1AB or 2ab),
 Pr(rec btw disease and marker haplotype, i.e. 1ab or 2AB), Pr(1Ab or
 2aB), Pr(1aB or 2Ab), where parental haplotypes are 1AB and 2ab.

 The output of get.2loglike.recomb() is the log-likelihood of the data
given
 the parameter values.*/
/* x is 16x4, ith row of x is i-1 base 2.*/

void get2loglikerecomb(double *data,double *s, double *p,double *am,double
*co,double *mp,double *loglike,int **x, int **y, int **z,
int **code, double *prpheno0, double **prpheno1, double **prpheno2,
double *mpp,double *denom)
{

int i,j,k;

/*prpheno0[i] is prob. of (A,B,a,b)^(x[i,]) and no sperm present*/

/*printf("prpheno0\n");*/
for(i=0;i<16;i++){
        prpheno0[i]=(1.-p[0]-p[1]);
        for(j=0;j<4;j++) prpheno0[i]*=(pow(co[j],(double)x[i][j])*
        pow(1.-co[j],1.-(double)x[i][j]));
        /*printf("%e ",prpheno0[i]);*/
}

mpp[0]=mp[1]; mpp[1]=mp[0]; mpp[2]=mp[3]; mpp[3]=mp[2];
/*printf("\nmpp\n");
for(i=0;i<4;i++) printf("%lf ",mpp[i]);*/

/* Now find prob. of each phenotype event and each genotype involving
1 sperm present.  prpheno1 is 16x8 matrix, ith row corresponds to
phenotype x[i,].  kth column corresponds to genotype y[k,], where 1st
digit of y is 1=mutant allele or 2=normal, 2nd digit is 1=AB, 2=ab,
3=Ab, 4=aB (given by code matrix). */

/*printf("\nprpheno1\n");*/
for(k=0;k<4;k++) for(i=0;i<16;i++){
        prpheno1[i][k]=p[0]**s*mp[k];
        prpheno1[i][4+k]=p[0]*(1.-*s)*mpp[k];
}
for(k=0;k<8;k++) for(i=0;i<16;i++) for(j=0;j<4;j++)
prpheno1[i][k]*=(pow(1.-pow(1.-am[j],(double)code[y[k][1]][j])*
(1.-co[j]),(double)x[i][j])*pow(pow(1.-am[j],(double)code[y[k][1]][j])*
(1.-co[j]),1.-(double)x[i][j]));

/*for(i=0;i<16;i++){
        for(k=0;k<8;k++) printf("%e ",prpheno1[i][k]);
        printf("\n");
}*/

/* Now find prob. of each phenotype event and each genotype involving
2 sperm present.  prpheno2 is 16x64 matrix, ith row corresponds to
phenotype x[i,].  kth column corresponds to genotype z[k,], where 1st
digit of z is 1st sperm: 1=mutant allele or 2=normal, 2nd digit is
2nd sperm: same, 3rd digit is 1st sperm: 1=AB, 2=ab,
3=Ab, 4=aB (given by code matrix), 4th digit is 2nd sperm: same */

for(k=0;k<64;k++) for(i=0;i<16;i++){
        prpheno2[i][k]=p[1]*pow(*s, 2.-(double)(z[k][0]+z[k][1]))*
        pow(1.-*s,(double)(z[k][0]+z[k][1]))*pow(mp[z[k][2]],
        1.-(double)z[k][0])*pow(mpp[z[k][2]],(double)z[k][0])*
        pow(mp[z[k][3]],1.-(double)z[k][1])*
        pow(mpp[z[k][3]],(double)z[k][1]);

        for(j=0;j<4;j++) prpheno2[i][k]*=(pow(1.-pow(1.-am[j],
        (double)(code[z[k][2]][j]+code[z[k][3]][j]))*(1.-co[j]),
        (double)x[i][j])*pow(pow(1.-am[j],(double)(code[z[k][2]][j]+
        code[z[k][3]][j]))*(1.-co[j]),1.-(double)x[i][j]));
}

/*printf("\nprpheno2\n");
for(i=0;i<16;i++){
        for(j=0;j<64;j++) printf("%e ",prpheno2[i][j]);
        printf("\n");
}*/

/*printf("denom\n");*/
for(i=0;i<16;i++){
        denom[i]=prpheno0[i];
        for(j=0;j<8;j++) denom[i]+=prpheno1[i][j];
        for(j=0;j<64;j++) denom[i]+=prpheno2[i][j];
        /*printf("%e ",denom[i]);*/
}

*loglike=0.;
for(i=0;i<16;i++) if(denom[i]>0.) *loglike+=(data[i]*log(denom[i]));

}

/*The function get.1loglike.recomb() takes as input the 1-locus sperm
typing
 data counts (length 4) in order given above and the current
 parameter values (length 7) in order given above. Also takes as input
 the
 recombination fraction between the gene and the marker.

 The output of get.1loglike.recomb() is the log-likelihood of the data
given
 the parameter values.*/

void get1loglikerecomb(double *data,double *s,double *p,double *am,double
*co,double *rf,double *loglike,int **x,int **y,int **z, int **code,
double *prpheno0, double **prpheno1,double **prpheno2,double *denom)
{

int i,j,k;

/* prpheno0[i] is prob of (A,a)^x[i,] and no sperm present */

/*printf("prpheno0\n");*/
for(i=0;i<4;i++){
        prpheno0[i]=(1.-p[0]-p[1]);
        for(j=0;j<2;j++) prpheno0[i]*=(pow(co[j],(double)x[i][j])*
        pow(1.-co[j],1.-(double)x[i][j]));
        /*printf("%e ",prpheno0[i]);
        fflush(stdout);*/
}

/* Now find prob. of each phenotype event and each genotype involving
 1 sperm present.  prpheno1 is 4x4 matrix, ith row corresponds to
 phenotype x[i,].  kth column corresponds to genotype y[k,], where 1st
 digit of y is 1=mutant allele or 2=normal, 2nd digit is 1=A, 2=a
 (given by code matrix).*/

/*printf("\nprpheno1\n");*/
for(i=0;i<4;i++){
        prpheno1[i][0]=p[0]**s*(1.-*rf);
        prpheno1[i][1]=p[0]**s**rf;
        prpheno1[i][2]=p[0]*(1.-*s)**rf;
        prpheno1[i][3]=p[0]*(1.-*s)*(1.-*rf);
}
for(k=0;k<4;k++) for(i=0;i<4;i++) for(j=0;j<2;j++)
prpheno1[i][k]*=(pow(1.-pow(1.-am[j],(double)code[y[k][1]][j])*
(1.-co[j]),(double)x[i][j])*pow(pow(1.-am[j],(double)code[y[k][1]][j])*
(1.-co[j]),1.-(double)x[i][j]));

/*for(i=0;i<4;i++){
        for(k=0;k<4;k++) printf("%e ",prpheno1[i][k]);
        printf("\n");
        fflush(stdout);
}*/

/* Now find prob. of each phenotype event and each genotype involving
 2 sperm present.  prpheno2 is 4x16 matrix, ith row corresponds to
 phenotype x[i,].  kth column corresponds to genotype z[k,], where 1st
 digit of z is 1st sperm: 1=mutant allele or 2=normal, 2nd digit is
 2nd sperm: same, 3rd digit is 1st sperm: 1=A, 2=a
 (given by code matrix), 4th digit is 2nd sperm: same */

for(k=0;k<16;k++) for(i=0;i<4;i++){
        prpheno2[i][k]=p[1]*pow(*s, 2.-(double)(z[k][0]+z[k][1]))*
        pow(1.-*s,(double)(z[k][0]+z[k][1]))*pow(*rf,fmod(fabs((double)

(z[k][0]-z[k][2])),2.)+fmod(fabs((double)(z[k][1]-z[k][3])),2.))*
        pow(1.-*rf,fmod(fabs((double)(z[k][0]-z[k][2]+1)),2.)+
        fmod(fabs((double)(z[k][1]-z[k][3]+1)),2.));

        for(j=0;j<2;j++) prpheno2[i][k]*=(pow(1.-pow(1.-am[j],
        (double)(code[z[k][2]][j]+code[z[k][3]][j]))*(1.-co[j]),
        (double)x[i][j])*pow(pow(1.-am[j],(double)(code[z[k][2]][j]+
        code[z[k][3]][j]))*(1.-co[j]),1.-(double)x[i][j]));
}

/*printf("\nprpheno2\n");
for(i=0;i<4;i++){
        for(j=0;j<16;j++) printf("%e ",prpheno2[i][j]);
        printf("\n");
        fflush(stdout);
}*/

/*printf("denom\n");*/
for(i=0;i<4;i++){
        denom[i]=prpheno0[i];
        for(j=0;j<4;j++) denom[i]+=prpheno1[i][j];
        for(j=0;j<16;j++) denom[i]+=prpheno2[i][j];
        /*printf("%e ",denom[i]);
        fflush(stdout);*/
}

*loglike=0.;
for(i=0;i<4;i++) if(denom[i]>0.) *loglike+=(data[i]*log(denom[i]));
}

/* The function get2expcellprobs() takes as input the current
 parameter values (length 11) in order given above. Also takes as input
 the
 multilocus recombination probabilities:  Pr(no rec., i.e. 1AB or 2ab),
 Pr(rec btw disease and marker haplotype, i.e. 1ab or 2AB), Pr(1Ab or
 2aB), Pr(1aB or 2Ab), where parental haplotypes are 1AB and 2ab.
 The output of get2expcellprobs() is the cell probabilities for
the data given the parameter values.*/

/* x is 16x4, ith row of x is i-1 base 2.*/

void get2expcellprobs(double *s, double *p,double *am,double
*co,double *mp,int **x, int **y, int **z, int **code, 
double *prpheno0, double **prpheno1, double **prpheno2,
double *mpp,double *denom,double *expcellprobs)
{

int i,j,k;

/*prpheno0[i] is prob. of (A,B,a,b)^(x[i,]) and no sperm present*/

/*printf("prpheno0\n");*/
for(i=0;i<16;i++){
        prpheno0[i]=(1.-p[0]-p[1]);
        for(j=0;j<4;j++) prpheno0[i]*=(pow(co[j],(double)x[i][j])*
        pow(1.-co[j],1.-(double)x[i][j]));
        /*printf("%e ",prpheno0[i]);*/
}

mpp[0]=mp[1]; mpp[1]=mp[0]; mpp[2]=mp[3]; mpp[3]=mp[2];
/*printf("\nmpp\n");
for(i=0;i<4;i++) printf("%lf ",mpp[i]);*/

/* Now find prob. of each phenotype event and each genotype involving
1 sperm present.  prpheno1 is 16x8 matrix, ith row corresponds to
phenotype x[i,].  kth column corresponds to genotype y[k,], where 1st
digit of y is 1=mutant allele or 2=normal, 2nd digit is 1=AB, 2=ab,
3=Ab, 4=aB (given by code matrix). */

/*printf("\nprpheno1\n");*/
for(k=0;k<4;k++) for(i=0;i<16;i++){
        prpheno1[i][k]=p[0]**s*mp[k];
        prpheno1[i][4+k]=p[0]*(1.-*s)*mpp[k];
}
for(k=0;k<8;k++) for(i=0;i<16;i++) for(j=0;j<4;j++)
prpheno1[i][k]*=(pow(1.-pow(1.-am[j],(double)code[y[k][1]][j])*
(1.-co[j]),(double)x[i][j])*pow(pow(1.-am[j],(double)code[y[k][1]][j])*
(1.-co[j]),1.-(double)x[i][j]));

/*for(i=0;i<16;i++){
        for(k=0;k<8;k++) printf("%e ",prpheno1[i][k]);
        printf("\n");
}*/

/* Now find prob. of each phenotype event and each genotype involving
2 sperm present.  prpheno2 is 16x64 matrix, ith row corresponds to
phenotype x[i,].  kth column corresponds to genotype z[k,], where 1st
digit of z is 1st sperm: 1=mutant allele or 2=normal, 2nd digit is
2nd sperm: same, 3rd digit is 1st sperm: 1=AB, 2=ab,
3=Ab, 4=aB (given by code matrix), 4th digit is 2nd sperm: same */

for(k=0;k<64;k++) for(i=0;i<16;i++){
        prpheno2[i][k]=p[1]*pow(*s, 2.-(double)(z[k][0]+z[k][1]))*
        pow(1.-*s,(double)(z[k][0]+z[k][1]))*pow(mp[z[k][2]],
        1.-(double)z[k][0])*pow(mpp[z[k][2]],(double)z[k][0])*
        pow(mp[z[k][3]],1.-(double)z[k][1])*
        pow(mpp[z[k][3]],(double)z[k][1]);

        for(j=0;j<4;j++) prpheno2[i][k]*=(pow(1.-pow(1.-am[j],
        (double)(code[z[k][2]][j]+code[z[k][3]][j]))*(1.-co[j]),
        (double)x[i][j])*pow(pow(1.-am[j],(double)(code[z[k][2]][j]+
        code[z[k][3]][j]))*(1.-co[j]),1.-(double)x[i][j]));
}

/*printf("\nprpheno2\n");
for(i=0;i<16;i++){
        for(j=0;j<64;j++) printf("%e ",prpheno2[i][j]);
        printf("\n");
}*/

/*printf("denom\n");*/
for(i=0;i<16;i++){
        denom[i]=prpheno0[i];
        for(j=0;j<8;j++) denom[i]+=prpheno1[i][j];
        for(j=0;j<64;j++) denom[i]+=prpheno2[i][j];
        /*printf("%e ",denom[i]);*/
}

for(i=0;i<16;i++) expcellprobs[i]=denom[i];
}

/* The function get1expcellprobs() takes as input the current
 parameter values (length 7) in order given above. Also takes as input
 the recombination fraction between the gene and the marker.
 The output of get1expcellprobs() is the cell probabilities for
the data given the parameter values.*/

void get1expcellprobs(double *s, double *p,double *am,double
*co,double *rf,int **x, int **y, int **z, int **code,
double *prpheno0, double **prpheno1, double **prpheno2,
double *denom,double *expcellprobs)
{

int i,j,k;

/* prpheno0[i] is prob of (A,a)^x[i,] and no sperm present */

/*printf("prpheno0\n");*/
for(i=0;i<4;i++){
        prpheno0[i]=(1.-p[0]-p[1]);
        for(j=0;j<2;j++) prpheno0[i]*=(pow(co[j],(double)x[i][j])*
        pow(1.-co[j],1.-(double)x[i][j]));
        /*printf("%e ",prpheno0[i]);
        fflush(stdout);*/
}

/* Now find prob. of each phenotype event and each genotype involving
 1 sperm present.  prpheno1 is 4x4 matrix, ith row corresponds to
 phenotype x[i,].  kth column corresponds to genotype y[k,], where 1st
 digit of y is 1=mutant allele or 2=normal, 2nd digit is 1=A, 2=a
 (given by code matrix).*/

/*printf("\nprpheno1\n");*/
for(i=0;i<4;i++){
        prpheno1[i][0]=p[0]**s*(1.-*rf);
        prpheno1[i][1]=p[0]**s**rf;
        prpheno1[i][2]=p[0]*(1.-*s)**rf;
        prpheno1[i][3]=p[0]*(1.-*s)*(1.-*rf);
}
for(k=0;k<4;k++) for(i=0;i<4;i++) for(j=0;j<2;j++)
prpheno1[i][k]*=(pow(1.-pow(1.-am[j],(double)code[y[k][1]][j])*
(1.-co[j]),(double)x[i][j])*pow(pow(1.-am[j],(double)code[y[k][1]][j])*
(1.-co[j]),1.-(double)x[i][j]));

/*for(i=0;i<4;i++){
        for(k=0;k<4;k++) printf("%e ",prpheno1[i][k]);
        printf("\n");
        fflush(stdout);
}*/

/* Now find prob. of each phenotype event and each genotype involving
 2 sperm present.  prpheno2 is 4x16 matrix, ith row corresponds to
 phenotype x[i,].  kth column corresponds to genotype z[k,], where 1st
 digit of z is 1st sperm: 1=mutant allele or 2=normal, 2nd digit is
 2nd sperm: same, 3rd digit is 1st sperm: 1=A, 2=a
 (given by code matrix), 4th digit is 2nd sperm: same */

for(k=0;k<16;k++) for(i=0;i<4;i++){
        prpheno2[i][k]=p[1]*pow(*s, 2.-(double)(z[k][0]+z[k][1]))*
        pow(1.-*s,(double)(z[k][0]+z[k][1]))*pow(*rf,fmod(fabs((double)

(z[k][0]-z[k][2])),2.)+fmod(fabs((double)(z[k][1]-z[k][3])),2.))*
        pow(1.-*rf,fmod(fabs((double)(z[k][0]-z[k][2]+1)),2.)+
        fmod(fabs((double)(z[k][1]-z[k][3]+1)),2.));

        for(j=0;j<2;j++) prpheno2[i][k]*=(pow(1.-pow(1.-am[j],
        (double)(code[z[k][2]][j]+code[z[k][3]][j]))*(1.-co[j]),
        (double)x[i][j])*pow(pow(1.-am[j],(double)(code[z[k][2]][j]+
        code[z[k][3]][j]))*(1.-co[j]),1.-(double)x[i][j]));
}

/*printf("\nprpheno2\n");
for(i=0;i<4;i++){
        for(j=0;j<16;j++) printf("%e ",prpheno2[i][j]);
        printf("\n");
        fflush(stdout);
}*/

/*printf("denom\n");*/
for(i=0;i<4;i++){
        denom[i]=prpheno0[i];
        for(j=0;j<4;j++) denom[i]+=prpheno1[i][j];
        for(j=0;j<16;j++) denom[i]+=prpheno2[i][j];
        /*printf("%e ",denom[i]);
        fflush(stdout);*/
}

for(i=0;i<4;i++) expcellprobs[i]=denom[i];
}


#include <malloc.h>
#include <stdio.h>

void nrerror(error_text)
char error_text[];
{
	void exit();

	fprintf(stderr,"Numerical Recipes run-time error...\n");
	fprintf(stderr,"%s\n",error_text);
	fprintf(stderr,"...now exiting to system...\n");
	exit(1);
}



double *vector(nl,nh)
int nl,nh;
{
	double *v;

	v=(double *)malloc((unsigned) (nh-nl+1)*sizeof(double));
	if (!v) nrerror("allocation failure in vector()");
	return v-nl;
}

int *ivector(nl,nh)
int nl,nh;
{
	int *v;

	v=(int *)malloc((unsigned) (nh-nl+1)*sizeof(int));
	if (!v) nrerror("allocation failure in ivector()");
	return v-nl;
}

double *dvector(nl,nh)
int nl,nh;
{
	double *v;

	v=(double *)malloc((unsigned) (nh-nl+1)*sizeof(double));
	if (!v) nrerror("allocation failure in dvector()");
	return v-nl;
}



double **matrix(nrl,nrh,ncl,nch)
int nrl,nrh,ncl,nch;
{
	int i;
	double **m;

	m=(double **) malloc((unsigned) (nrh-nrl+1)*sizeof(double*));
	if (!m) nrerror("allocation failure 1 in matrix()");
	m -= nrl;

	for(i=nrl;i<=nrh;i++) {
		m[i]=(double *) malloc((unsigned) (nch-ncl+1)*sizeof(double));
		if (!m[i]) nrerror("allocation failure 2 in matrix()");
		m[i] -= ncl;
	}
	return m;
}

double **dmatrix(nrl,nrh,ncl,nch)
int nrl,nrh,ncl,nch;
{
	int i;
	double **m;

	m=(double **) malloc((unsigned) (nrh-nrl+1)*sizeof(double*));
	if (!m) nrerror("allocation failure 1 in dmatrix()");
	m -= nrl;

	for(i=nrl;i<=nrh;i++) {
		m[i]=(double *) malloc((unsigned) (nch-ncl+1)*sizeof(double));
		if (!m[i]) nrerror("allocation failure 2 in dmatrix()");
		m[i] -= ncl;
	}
	return m;
}

int **imatrix(nrl,nrh,ncl,nch)
int nrl,nrh,ncl,nch;
{
	int i,**m;

	m=(int **)malloc((unsigned) (nrh-nrl+1)*sizeof(int*));
	if (!m) nrerror("allocation failure 1 in imatrix()");
	m -= nrl;

	for(i=nrl;i<=nrh;i++) {
		m[i]=(int *)malloc((unsigned) (nch-ncl+1)*sizeof(int));
		if (!m[i]) nrerror("allocation failure 2 in imatrix()");
		m[i] -= ncl;
	}
	return m;
}



double **submatrix(a,oldrl,oldrh,oldcl,oldch,newrl,newcl)
double **a;
int oldrl,oldrh,oldcl,oldch,newrl,newcl;
{
	int i,j;
	double **m;

	m=(double **) malloc((unsigned) (oldrh-oldrl+1)*sizeof(double*));
	if (!m) nrerror("allocation failure in submatrix()");
	m -= newrl;

	for(i=oldrl,j=newrl;i<=oldrh;i++,j++) m[j]=a[i]+oldcl-newcl;

	return m;
}



void free_vector(v,nl,nh)
double *v;
int nl,nh;
{
	free((char*) (v+nl));
}

void free_ivector(v,nl,nh)
int *v,nl,nh;
{
	free((char*) (v+nl));
}

void free_dvector(v,nl,nh)
double *v;
int nl,nh;
{
	free((char*) (v+nl));
}



void free_matrix(m,nrl,nrh,ncl,nch)
double **m;
int nrl,nrh,ncl,nch;
{
	int i;

	for(i=nrh;i>=nrl;i--) free((char*) (m[i]+ncl));
	free((char*) (m+nrl));
}

void free_dmatrix(m,nrl,nrh,ncl,nch)
double **m;
int nrl,nrh,ncl,nch;
{
	int i;

	for(i=nrh;i>=nrl;i--) free((char*) (m[i]+ncl));
	free((char*) (m+nrl));
}

void free_imatrix(m,nrl,nrh,ncl,nch)
int **m;
int nrl,nrh,ncl,nch;
{
	int i;

	for(i=nrh;i>=nrl;i--) free((char*) (m[i]+ncl));
	free((char*) (m+nrl));
}



void free_submatrix(b,nrl,nrh,ncl,nch)
double **b;
int nrl,nrh,ncl,nch;
{
	free((char*) (b+nrl));
}



double **convert_matrix(a,nrl,nrh,ncl,nch)
double *a;
int nrl,nrh,ncl,nch;
{
	int i,j,nrow,ncol;
	double **m;

	nrow=nrh-nrl+1;
	ncol=nch-ncl+1;
	m = (double **) malloc((unsigned) (nrow)*sizeof(double*));
	if (!m) nrerror("allocation failure in convert_matrix()");
	m -= nrl;
	for(i=0,j=nrl;i<=nrow-1;i++,j++) m[j]=a+ncol*i-ncl;
	return m;
}



void free_convert_matrix(b,nrl,nrh,ncl,nch)
double **b;
int nrl,nrh,ncl,nch;
{
	free((char*) (b+nrl));
}
#include <math.h>
#include <stdio.h>

/* sfix is a vector of length ndon indicating which of the s variables
are
 to be fixed (1 means fix it, 0 means let it vary).  The
corresponding 
 values are given in sval (vector of length ndon). The entry for a
 variable will be ignored if that variable is not fixed.  pfix is a 
 ndonx2 matrix, similarly for pval. amfix, amval,
 cofix, coval are vectors of length nal. sstart, pstar, amstart,
costart are of the same type as s, p, am, and co, and contain the
starting values for these parameters.  Now the probabilities of 1 sperm
deposited and the probabilities of 2 sperm deposited can have equality
relationships that are separately specified.  Also have changed the way
the EM algorithm decides when it's done.*/

void spermanalysisrfssdem(int ndon,int nal,int n2loc,double **twoloc,
int *sind2,int **pind2,int **amal2,int **coal2,double **twolocmp,
int n1loc,double **oneloc,int *sind1,int **pind1,int **amal1,int **coal1,
double *onelocrf,int **x,int **y,int **z,int **code,double *output,
double *prpheno0,double **prpheno1,double **prpheno2,double *mpp,
double *denom,double *tmp,double *tmpt,double *ltmp,double *s,
double **p,double *am,double *co,double *pin,double *amin,double *coin,
double **prpheno11,double **prpheno21,int **x1,int **y1,int **z1,
int **code1,double *loglike,double **exp2,double **exp1,double *sden,
double **pden,double *amden,double *coden,double *maxloglike,int *sfix,
int **pfix,int *amfix,int *cofix,double *sval,double **pval,
double *amval,double *coval,double *sstart,double **pstart,double
*amstart,double *costart,double *dconv,int *niter,int simind,int
**ptype,double *dtmp)
{
int i,j,k,l, iconv;

void makexyzcode(int **x,int **y,int **z,int **code);

void makeonexyzcode(int **x1,int **y1,int **z1,int **code1);

void get2locexprecomb(double *data,double *s, double *p, double *am,
double *co, double *mp, double *output, int **x, int **y, int **z,
int **code, double *prpheno0, double **prpheno1, double **prpheno2,
double *mpp,double *denom, double *tmp, double *tmpt, double *ltmp);

void get1locexprecomb(double *data,double *s,double *p,double *am,double
*co,double *rf,double *output,int **x,int **y,int **z, int **code,
double *prpheno0, double **prpheno1,double **prpheno2,double *denom,
double *tmp,double *tmpt,double *ltmp);

void get2loglikerecomb(double *data,double *s, double *p, double *am,
double *co,double *mp,double *loglike,int **x, int **y, int **z,
int **code, double *prpheno0, double **prpheno1, double **prpheno2,
double *mpp,double *denom);

void get1loglikerecomb(double *data,double *s,double *p,double *am,
double *co,double *rf,double *loglike,int **x,int **y,int **z, 
int **code, double *prpheno0, double **prpheno1,double **prpheno2,
double *denom);

makexyzcode(x,y,z,code);
makeonexyzcode(x1,y1,z1,code1);

for(i=0;i<ndon;i++){
	if(sfix[i]==1) s[i]=sval[i]; else s[i]=sstart[i];
	if(pfix[i][0]==1) p[i][0]=pval[i][0]; 
	else p[i][0]=pstart[i][0];
	if(pfix[i][1]==1) p[i][1]=pval[i][1]; 
	else p[i][1]=pstart[i][1];
}
/* Amplification an contamination params are in order A,B,a,b */
for(i=0;i<nal;i++){
        if(amfix[i]==1) am[i]=amval[i]; else am[i]=amstart[i];
        if(cofix[i]==1) co[i]=coval[i]; else co[i]=costart[i];
}

iconv=0;
*niter=0;
while(iconv==0){
*niter+=300;
for(k=0;k<300;k++){
	/*printf("d2 1 %f\n",p[1][1]);
	fflush(stdout);*/
	for(i=0;i<n2loc;i++){
		
		for(j=0;j<4;j++) amin[j]=am[amal2[i][j]];
		for(j=0;j<4;j++) coin[j]=co[coal2[i][j]];
		for(j=0;j<2;j++) pin[j]=p[pind2[i][j]][j];

/*if(*niter==2400){ 
printf("iteration %d\n",k);
printf("twoloc\n");
for(j=0;j<16;j++) printf("%f ",twoloc[i][j]);
printf("\ns p1 p2 a1 a2 a3 a4 c1 c2 c3 c4\n");*/
/*printf("%f %f %f %f %f %f %f %f %f %f %f\n",s[sind2[i]],
pin[0],pin[1],amin[0],amin[1],amin[2],amin[3],coin[0],coin[1],
coin[2],coin[3]);*/
/*printf("twolocmp %f %f %f %f\n",twolocmp[i][0],twolocmp[i][1],
twolocmp[i][2],twolocmp[i][3]);
}*/

		get2locexprecomb(twoloc[i],s+sind2[i],pin,
		amin,coin,twolocmp[i],output,x,y,z,code,prpheno0,
		prpheno1,prpheno2,mpp,denom,tmp,tmpt,ltmp);

/*if(*niter==2400){
printf("Output of get2locexprecomb\n");
for(j=0;j<21;j++) printf("%f ",output[j]);
printf("\n");
fflush(stdout);
}*/

		for(j=0;j<21;j++) exp2[i][j]=output[j];
	}

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

		for(j=0;j<2;j++) amin[j]=am[amal1[i][j]];
		for(j=0;j<2;j++) coin[j]=co[coal1[i][j]];
		for(j=0;j<2;j++) pin[j]=p[pind1[i][j]][j];

		get1locexprecomb(oneloc[i],s+sind1[i],pin,
		amin,coin,onelocrf+i,output,x1,y1,z1,code1,prpheno0,
		prpheno11,prpheno21,denom,tmp,tmpt,ltmp);

		for(j=0;j<13;j++) exp1[i][j]=output[j];
	}

	for(i=0;i<ndon;i++){
        	s[i]=0.;
		sden[i]=0.;
        	p[i][0]=0.;
        	pden[i][0]=0.;
        	p[i][1]=0.;
        	pden[i][1]=0.;
	}
	for(i=0;i<nal;i++){
        	am[i]=0.;
        	amden[i]=0.;
        	co[i]=0.;
        	coden[i]=0.;
	}

	for(j=0;j<n2loc;j++){
		s[sind2[j]]+=exp2[j][0];
		sden[sind2[j]]+=(exp2[j][0]+exp2[j][1]);
		p[pind2[j][0]][0]+=exp2[j][2];
		p[pind2[j][1]][1]+=exp2[j][3];
		for(i=0;i<2;i++) for(l=2;l<5;l++) 
		pden[pind2[j][i]][i]+=exp2[j][l];
		am[amal2[j][0]]+=exp2[j][5];
		am[amal2[j][1]]+=exp2[j][7];
		am[amal2[j][2]]+=exp2[j][9];
		am[amal2[j][3]]+=exp2[j][11];
		amden[amal2[j][0]]+=(exp2[j][5]+exp2[j][6]);
		amden[amal2[j][1]]+=(exp2[j][7]+exp2[j][8]);
		amden[amal2[j][2]]+=(exp2[j][9]+exp2[j][10]);
		amden[amal2[j][3]]+=(exp2[j][11]+exp2[j][12]);
		co[coal2[j][0]]+=exp2[j][13];
		co[coal2[j][1]]+=exp2[j][15];
		co[coal2[j][2]]+=exp2[j][17];
		co[coal2[j][3]]+=exp2[j][19];
		coden[coal2[j][0]]+=(exp2[j][13]+exp2[j][14]);
		coden[coal2[j][1]]+=(exp2[j][15]+exp2[j][16]);
		coden[coal2[j][2]]+=(exp2[j][17]+exp2[j][18]);
		coden[coal2[j][3]]+=(exp2[j][19]+exp2[j][20]);
	}

	for(j=0;j<n1loc;j++){
		s[sind1[j]]+=exp1[j][0];
                sden[sind1[j]]+=(exp1[j][0]+exp1[j][1]);
                p[pind1[j][0]][0]+=exp1[j][2];
                p[pind1[j][1]][1]+=exp1[j][3];
                for(i=0;i<2;i++) for(l=2;l<5;l++)
                pden[pind1[j][i]][i]+=exp1[j][l];
                am[amal1[j][0]]+=exp1[j][5];
                am[amal1[j][1]]+=exp1[j][7];
                amden[amal1[j][0]]+=(exp1[j][5]+exp1[j][6]);
                amden[amal1[j][1]]+=(exp1[j][7]+exp1[j][8]);
		/*printf("%lf %lf\n",co[coal1[j][0]],co[coal1[j][1]]);*/
                co[coal1[j][0]]+=exp1[j][9];
                co[coal1[j][1]]+=exp1[j][11];
		/*printf("%lf %lf\n",co[coal1[j][0]],co[coal1[j][1]]);*/
		/*printf("%lf %lf\n",coden[coal1[j][0]],coden[coal1[j][1]]);*/
                coden[coal1[j][0]]+=(exp1[j][9]+exp1[j][10]);
                coden[coal1[j][1]]+=(exp1[j][11]+exp1[j][12]);
		/*printf("%lf %lf\n",coden[coal1[j][0]],coden[coal1[j][1]]);*/
        }

	for(j=0;j<ndon;j++){
		if(sden[j]>0.) s[j]=s[j]/sden[j];

		if(sfix[j]==1) s[j]=sval[j];
	}

	for(i=0;i<ndon;i++) for(j=0;j<2;j++){
		if(ptype[i][4*j+2]==1){
			/*printf("i=%d, ptype[i][4*j+2]==1\n",i);*/
			p[i][j]=pval[i][j];
		}
		else if(ptype[i][4*j+1]==1 && ptype[i][4*j+3]==0){
/*printf("i=%d, ptype[i][4*j+1]==1 && ptype[i][4*j+3]==0\n",i);*/
			if(pden[i][j]>0.) p[i][j]=p[i][j]/pden[i][j];
		}
		else if(ptype[i][4*j+1]==1 && ptype[i][4*j+3]==1){

			dtmp[0]=pden[i][j];
			dtmp[2]=p[i][j];
			dtmp[3]=pden[i][j];
			dtmp[4]=p[i][j];
			for(l=0;l<n2loc && !(pind2[l][j]==i &&
			pfix[pind2[l][1-j]][1-j]==1);l++); 
			if(l<n2loc){
				dtmp[5]=pval[pind2[l][1-j]][1-j];
				dtmp[2]*=(1.-dtmp[5]);
			}
			else{ 
				for(l=0;l<n1loc && !(pind1[l][j]==i &&
				pfix[pind1[l][1-j]][1-j]==1);l++);
				if(l<n1loc){
                                	dtmp[5]=pval[pind1[l][1-j]][1-j];
                                	dtmp[2]*=(1.-dtmp[5]);
                        	}
			}
			for(l=0;l<n2loc;l++) if(pind2[l][j]==i &&
			pfix[pind2[l][1-j]][1-j]==1){
				dtmp[0]-=(exp2[l][3-j]);
				dtmp[3]-=(exp2[l][3-j]+exp2[l][4]);
				dtmp[4]+=(exp2[l][4]);
			}
			for(l=0;l<n1loc;l++) if(pind1[l][j]==i &&
			pfix[pind1[l][1-j]][1-j]==1){
                                	dtmp[0]-=(exp1[l][3-j]);
                                	dtmp[3]-=(exp1[l][3-j]+exp1[l][4]);
                                	dtmp[4]+=(exp1[l][4]);
			}
			dtmp[1]=-(1.-dtmp[5])*dtmp[3]-dtmp[4];
			if(dtmp[0]<=0.) p[i][j]=0.;
			else if((dtmp[3]=pow(dtmp[1],2.)-4.*dtmp[0]*
			dtmp[2])<0.) p[i][j]=1.-dtmp[5];
			else p[i][j]=(-dtmp[1]-pow(dtmp[3],.5))/
			(2.*dtmp[0]);
	/*printf("i=%d, ptype[i][4*j+1]==1 && ptype[i][4*j+3]==1\n",i);*/
		}
	}

	for(i=0;i<ndon;i++) for(j=0;j<2;j++){
		if(ptype[i][4*j+1]==0 && ptype[i][4*j]==1 &&
		ptype[i][4*j+2]==0){
/*printf("i=%d, ptype[i][4*j+1]==0 && ptype[i][4*j]==1 && ptype[i][4*j+2]==0\n",
i);*/
			p[i][j]=0.;
			for(l=0;l<n2loc;l++) if(pind2[l][j]==i){
				p[i][j]+=(exp2[l][j+2]*(1.-
			p[pind2[l][1-j]][1-j]));
				pden[i][j]-=(exp2[l][3-j]);
                /*printf("p[i][j] %f\n",p[i][j]);
                printf("pden[i][j] %f\n",pden[i][j]);*/

			}
			for(l=0;l<n1loc;l++) if(pind1[l][j]==i){
				p[i][j]+=(exp1[l][j+2]*(1.-
				p[pind1[l][1-j]][1-j]));
				pden[i][j]-=(exp1[l][3-j]);
                /*printf("p[i][j] %f\n",p[i][j]);
                printf("pden[i][j] %f\n",pden[i][j]);*/
			}
			if(pden[i][j]>0.) p[i][j]=p[i][j]/pden[i][j];
		}
	}
			
	for(j=0;j<nal;j++){
		if(amden[j]>0.) am[j]=am[j]/amden[j];
		if(coden[j]>0.) co[j]=co[j]/coden[j];

		if(amfix[j]==1) am[j]=amval[j];
		if(cofix[j]==1) co[j]=coval[j];
	}

/*for(i=0;i<ndon;i++) printf("%7e ",s[i]);
printf("\n");
for(i=0;i<ndon;i++){
        for(j=0;j<2;j++) printf("%7e ",p[i][j]);
        printf("\n");
}
for(i=0;i<nal;i++) printf("%7e ",am[i]);
printf("\n");
for(i=0;i<nal;i++) printf("%7e ",co[i]);
printf("\n\n");*/

/* compare the parameters after 250th and 300th iterations and look
for a change. */
if(k==249){
	j=0;
	for(i=0;i<ndon;i++){
		dconv[j]=s[i];
		j++;
	}
	for(i=0;i<ndon;i++) for(l=0;l<2;l++){
		dconv[j]=p[i][l];
		j++;
	}
	for(i=0;i<nal;i++){
		dconv[j]=am[i];
		j++;
	}
	for(i=0;i<nal;i++){
		dconv[j]=co[i];
		j++;
	}
*maxloglike=0.;
for(l=0;l<n2loc;l++){
 
        for(i=0;i<4;i++) amin[i]=am[amal2[l][i]];
        for(i=0;i<4;i++) coin[i]=co[coal2[l][i]];
        for(i=0;i<2;i++) pin[i]=p[pind2[l][i]][i];
 
        get2loglikerecomb(twoloc[l],s+sind2[l],pin,amin,
        coin,twolocmp[l],loglike,x,y,z,code,prpheno0,prpheno1,
        prpheno2,mpp,denom);
 
        *maxloglike+=*loglike;
        maxloglike[l+1]=*loglike;
}
 
for(i=0;i<n1loc;i++){
 
        for(l=0;l<2;l++) amin[l]=am[amal1[i][l]];
        for(l=0;l<2;l++) coin[l]=co[coal1[i][l]];
        for(l=0;l<2;l++) pin[l]=p[pind1[i][l]][l];
 
        get1loglikerecomb(oneloc[i],s+sind1[i],pin,amin,
        coin,onelocrf+i,loglike,x1,y1,z1,code1,prpheno0,prpheno11,
        prpheno21,denom);
 
        *maxloglike+=*loglike;
        maxloglike[n2loc+i+1]=*loglike;
}

	dconv[3*ndon+2*nal]=*maxloglike;
}

} /* End of loop over number of iterations */


iconv=1;
j=0;
for(i=0;i<ndon;i++){
	if(!(s[i]>pow(10.,-10.) || s[i]>dconv[j]) 
|| !(fabs(dconv[j]-s[i])/s[i]>pow(10.,-6.)) ){
/*|| 
!(fabs(dconv[j]-s[i])/s[i]>pow(10.,-3.) || fabs(dconv[j]-s[i])>
pow(10.,-6.))){*/
		iconv*=1;
		/*printf("%f %f\n",dconv[j],s[i]);*/
	}
	else{
		iconv*=0;
		/*printf("%d %f %f\n",*niter,dconv[j],s[i]);*/
	}
	j++;
}
for(i=0;i<ndon;i++) for(l=0;l<2;l++){
	if(!(p[i][l]>pow(10.,-10.) || p[i][l]>dconv[j])
|| !(fabs(dconv[j]-p[i][l])/p[i][l]>pow(10.,-6.))){
/* ||
!(fabs(dconv[j]-p[i][l])/p[i][l]>pow(10.,-3.) || fabs(dconv[j]-p[i][l])>
pow(10.,-6.))){*/
		iconv*=1;
		/*printf("%f %f\n",dconv[j],p[i][l]);*/
	}
	else{
		iconv*=0;
		/*printf("%d %e %e\n",*niter,dconv[j],p[i][l]);*/
	}
        j++;
}
for(i=0;i<nal;i++){
	if(!(am[i]>pow(10.,-10.) || am[i]>dconv[j])
|| !(fabs(dconv[j]-am[i])/am[i]>pow(10.,-6.))){
/* ||
!(fabs(dconv[j]-am[i])/am[i]>pow(10.,-3.) || fabs(dconv[j]-am[i])>
pow(10.,-6.))){*/
		iconv*=1;
		/*printf("%f %f\n",dconv[j],am[i]);*/
	}
	else{ 
		iconv*=0;
		/*printf("%d %f %f\n",*niter,dconv[j],am[i]);*/
	}
        j++;
}
for(i=0;i<nal;i++){
        if(!(co[i]>pow(10.,-10.) || co[i]>dconv[j])
|| !(fabs(dconv[j]-co[i])/co[i]>pow(10.,-6.))){
/* ||
!(fabs(dconv[j]-co[i])/co[i]>pow(10.,-3.) || fabs(dconv[j]-co[i])>
pow(10.,-6.))){*/
		iconv*=1;
		/*printf("%f %f\n",dconv[j],co[i]);*/
	}
	else{ 
		iconv*=0;
		/*printf("%d %e %e\n",*niter,dconv[j],co[i]);*/
	}
        j++;
}
fflush(stdout);

if(iconv==0 && simind==1){
*maxloglike=0.;
for(i=0;i<n2loc;i++){

        for(j=0;j<4;j++) amin[j]=am[amal2[i][j]];
        for(j=0;j<4;j++) coin[j]=co[coal2[i][j]];
        for(j=0;j<2;j++) pin[j]=p[pind2[i][j]][j];

        get2loglikerecomb(twoloc[i],s+sind2[i],pin,amin,
        coin,twolocmp[i],loglike,x,y,z,code,prpheno0,prpheno1,
        prpheno2,mpp,denom);
 
        *maxloglike+=*loglike;
        maxloglike[i+1]=*loglike;
}
 
for(i=0;i<n1loc;i++){
 
        for(j=0;j<2;j++) amin[j]=am[amal1[i][j]];
        for(j=0;j<2;j++) coin[j]=co[coal1[i][j]];
        for(j=0;j<2;j++) pin[j]=p[pind1[i][j]][j];
 
        get1loglikerecomb(oneloc[i],s+sind1[i],pin,amin,
        coin,onelocrf+i,loglike,x1,y1,z1,code1,prpheno0,prpheno11,
        prpheno21,denom);
 
        *maxloglike+=*loglike;
        maxloglike[n2loc+i+1]=*loglike;
}

if(fabs(dconv[ndon*3+2*nal]-*maxloglike)<pow(10.,-4.)) iconv=1;
}


if(*niter>90000){
printf("Warning: exiting EM algorithm --- too many iterations\n");
for(i=0;i<n2loc;i++){
	for(j=0;j<16;j++) printf("%f ",twoloc[i][j]);
	printf("\n");
}
for(i=0;i<n1loc;i++){
	for(j=0;j<4;j++) printf("%f ",oneloc[i][j]);
	printf("\n");
}
for(i=0;i<ndon;i++) printf("%f ",sstart[i]);
printf("\n");
for(i=0;i<ndon;i++) printf("%f ",pstart[i][0]);
printf("\n");
for(i=0;i<ndon;i++) printf("%f ",pstart[i][1]);
printf("\n");
for(i=0;i<nal;i++) printf("%f ",amstart[i]);
printf("\n");
for(i=0;i<nal;i++) printf("%f ",costart[i]);
printf("\n");
	break;
}

}/* end of while statement on iconv */

/* Now calculate the likelihood */

*maxloglike=0.;
for(i=0;i<n2loc;i++){

	for(j=0;j<4;j++) amin[j]=am[amal2[i][j]];
        for(j=0;j<4;j++) coin[j]=co[coal2[i][j]];
	for(j=0;j<2;j++) pin[j]=p[pind2[i][j]][j];

	get2loglikerecomb(twoloc[i],s+sind2[i],pin,amin,
	coin,twolocmp[i],loglike,x,y,z,code,prpheno0,prpheno1,
	prpheno2,mpp,denom);

	*maxloglike+=*loglike;
	maxloglike[i+1]=*loglike;
}

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

        for(j=0;j<2;j++) amin[j]=am[amal1[i][j]];
        for(j=0;j<2;j++) coin[j]=co[coal1[i][j]];
	for(j=0;j<2;j++) pin[j]=p[pind1[i][j]][j];

	get1loglikerecomb(oneloc[i],s+sind1[i],pin,amin,
	coin,onelocrf+i,loglike,x1,y1,z1,code1,prpheno0,prpheno11,
	prpheno21,denom);

	*maxloglike+=*loglike;
	maxloglike[n2loc+i+1]=*loglike;
}

if(!(*maxloglike<0.)){
for(i=0;i<n2loc;i++){
        for(j=0;j<16;j++) printf("%f ",twoloc[i][j]);
        printf("\n");
}
for(i=0;i<n1loc;i++){
        for(j=0;j<4;j++) printf("%f ",oneloc[i][j]);
        printf("\n");
}
for(i=0;i<ndon;i++) printf("%f ",sstart[i]);
printf("\n");
for(i=0;i<ndon;i++) printf("%f ",pstart[i][0]);
printf("\n");
for(i=0;i<ndon;i++) printf("%f ",pstart[i][1]);
printf("\n");
for(i=0;i<nal;i++) printf("%f ",amstart[i]);
printf("\n");
for(i=0;i<nal;i++) printf("%f ",costart[i]);
printf("\n");
}

}
#include <math.h>
#include <stdio.h>

/* maxvals[3*nd(n+2*nal+1] is the vector of parameter values that 
maximize the
likelihood, with last entry the maximized likelihood, quan is
qchisq(.95,1) for a 95% confidence interval, qchisq(.68,1) to get
standard errors */

void spermanalysisfssdioemcis(int ndon,int nal,int n2loc,double **twoloc,
int *sind2,int **pind2,int **amal2,int **coal2,double **twolocmp,
int n1loc,double **oneloc,int *sind1,int **pind1,int **amal1,int **coal1,
double *onelocrf,int **x,int **y,int **z,int **code,double *output,
double *prpheno0,double **prpheno1,double **prpheno2,double *mpp,
double *denom,double *tmp,double *tmpt,double *ltmp,double *s,
double **p,double *am,double *co,double *pin,double *amin,double *coin,
double **prpheno11,double **prpheno21,int **x1,int **y1,int **z1,
int **code1,double *loglike,double **exp2,double **exp1,double *sden,
double **pden,double *amden,double *coden,double *maxloglike,int *sfix,
int **pfix,int *amfix,int *cofix,double *sval,double **pval,
double *amval,double *coval,double *maxvals,double *quan,
int *sfix2,int **pfix2,int *amfix2,int *cofix2,double *sval2,
double **pval2,double *amval2,double *coval2,double *sstart,double
**pstart,double *amstart,double *costart,double *dconv,int *niter, int
nstep1,int nstep2,int nstep3,double *sstart2,double
**pstart2,double *amstart2,double *costart2,int **ptype,int **ptype2,
double *dtmp)
{
int i,j,k,l,ll,m,nstep,nt;
double eps,min,max,a,b,upp,low;

void spermanalysisrfssdem(int ndon,int nal,int n2loc,double **twoloc,
int *sind2,int **pind2,int **amal2,int **coal2,double **twolocmp,
int n1loc,double **oneloc,int *sind1,int **pind1,int **amal1,int **coal1,
double *onelocrf,int **x,int **y,int **z,int **code,double *output,
double *prpheno0,double **prpheno1,double **prpheno2,double *mpp,
double *denom,double *tmp,double *tmpt,double *ltmp,double *s,
double **p,double *am,double *co,double *pin,double *amin,double *coin,
double **prpheno11,double **prpheno21,int **x1,int **y1,int **z1,
int **code1,double *loglike,double **exp2,double **exp1,double *sden,
double **pden,double *amden,double *coden,double *maxloglike,int *sfix,
int **pfix,int *amfix,int *cofix,double *sval,double **pval,
double *amval,double *coval,double *sstart,double **pstart,double
*amstart,double *costart,double *dconv,int *niter,int simind,int
**ptype,double *dtmp);

eps=pow(10.,-10.);
min=eps;
max=1-eps;

for(i=0;i<ndon;i++){
sfix2[i]=sfix[i];
sval2[i]=sval[i];
pfix2[i][0]=pfix[i][0];
pval2[i][0]=pval[i][0];
pfix2[i][1]=pfix[i][1];
pval2[i][1]=pval[i][1];
}

for(i=0;i<nal;i++){
amfix2[i]=amfix[i];
amval2[i]=amval[i];
cofix2[i]=cofix[i];
coval2[i]=coval[i];
}

for(i=0;i<ndon;i++) for(j=0;j<8;j++) ptype2[i][j]=ptype[i][j];

for(i=0;i<ndon;i++){
sstart2[i]=maxvals[i];
for(j=0;j<2;j++) pstart2[i][j]=maxvals[ndon+2*i+j];
}
for(i=0;i<nal;i++){
amstart2[i]=maxvals[3*ndon+i];
costart2[i]=maxvals[3*ndon+nal+i];
}

printf("Confidence intervals for estimated parameters\n\n");
printf("Segregation parameters\n\n");
printf("Parameter     Estimate     Interval\n");

for(i=0;i<0;i++){
/*for(i=0;i<ndon;i++){*/

for(j=0;j<ndon;j++){ 
	sfix2[j]=sfix[j];
	sval2[j]=sval[j];
}

nt=1;
for(j=0;j<n2loc;j++) if(sind2[j]==i && sfix[i]==0) nt=0;
for(j=0;j<n1loc;j++) if(sind1[j]==i && sfix[i]==0) nt=0;
if(nt==0){
	sfix2[i]=1;
	a=maxvals[i];
	b=max;
	if(a>b) upp=a; else{
		for(j=0;j<nstep1;j++){
			sval2[i]=(a+b)/2.;
			spermanalysisrfssdem(ndon,nal,n2loc,
twoloc,sind2,pind2,amal2,coal2,twolocmp,n1loc,oneloc,sind1,pind1,amal1,
coal1,onelocrf,x,y,z,code,output,prpheno0,prpheno1,prpheno2,mpp,denom,
tmp,tmpt,ltmp,s,p,am,co,pin,amin,
coin,prpheno11,prpheno21,x1,y1,z1,code1,loglike,exp2,exp1,sden,
pden,amden,coden,maxloglike,sfix2,pfix2,amfix2,cofix2,sval2,pval2,
amval2,coval2,sstart2,pstart2,amstart2,costart2,dconv,niter,1,ptype,
dtmp);

			/*for(k=0;k<ndon;k++) printf("%7e ",s[k]);
			printf("\n");
			for(k=0;k<ndon;k++){
				for(l=0;l<2;l++) 
				printf("%7e ",p[k][l]);
				printf("\n");
			}
			for(k=0;k<nal;k++) printf("%7e ",am[k]);
			printf("\n");
			for(k=0;k<nal;k++) printf("%7e ",co[k]);
			printf("\n");
			printf("%7e\n",*maxloglike);*/

			if(*maxloglike>maxvals[3*ndon+2*nal]-
			*quan/2.) a=(a+b)/2.;
			else b=(a+b)/2.;
			/*printf("%7e %7e\n",a,b);
			fflush(stdout);*/
		}
		upp=b;
	}
	/*printf("%7e\n",upp);*/

	a=min;
	b=maxvals[i];
	if(a>b) low=b; else{
		for(j=0;j<nstep1;j++){
			sval2[i]=(a+b)/2.;
                        spermanalysisrfssdem(ndon,nal,n2loc,
twoloc,sind2,pind2,amal2,coal2,twolocmp,n1loc,oneloc,sind1,pind1,amal1,
coal1,onelocrf,x,y,z,code,output,prpheno0,prpheno1,prpheno2,mpp,denom,
tmp,tmpt,ltmp,s,p,am,co,pin,amin,
coin,prpheno11,prpheno21,x1,y1,z1,code1,loglike,exp2,exp1,sden,
pden,amden,coden,maxloglike,sfix2,pfix2,amfix2,cofix2,sval2,pval2,
amval2,coval2,sstart2,pstart2,amstart2,costart2,dconv,niter,1,ptype,
dtmp);

			/*for(k=0;k<ndon;k++) printf("%7e ",s[i]);
			printf("\n");
			for(k=0;k<ndon;k++){
				for(l=0;l<2;l++)
				printf("%7e ",p[k][l]);
				printf("\n");
			}
			for(k=0;k<nal;k++) printf("%7e ",am[k]);
			printf("\n");
			for(k=0;k<nal;k++) printf("%7e ",co[k]);
			printf("\n");
			printf("%7e\n",*maxloglike);*/

			if(*maxloglike>maxvals[3*ndon+2*nal]-
			*quan/2.) b=(a+b)/2.;
			else a=(a+b)/2.;
			/*printf("%7e %7e\n",a,b); 
			fflush(stdout);*/
		} 
		low=a;
	}
	/*printf("%7e\n",low);*/

	printf("s %4d %15f %9f %7f\n",i,maxvals[i],low,upp);
	fflush(stdout);
}
}

for(j=0;j<ndon;j++){ 
	sfix2[j]=sfix[j];
	sval2[j]=sval[j];
}

printf("\nDeposit parameters\n\n");
printf("Parameter     Estimate     Interval\n");

for(i=0;i<ndon;i++){
	/*for(j=0;j<2;j++) if(pfix[i][j]==0){*/
	for(j=0;j<1;j++) if(pfix[i][j]==0){
		nt=1;
		b=max;
		for(k=0;k<n2loc;k++) if(pind2[k][j]==i){
			nt=0;
			if(pfix[pind2[k][1-j]][1-j]==1)
if(b>1.-pval[pind2[k][1-j]][1-j]) b=1.-pval[pind2[k][1-j]][1-j];
		}
		for(k=0;k<n1loc;k++) if(pind1[k][j]==i){
			nt=0;
			if(pfix[pind1[k][1-j]][1-j]==1)
if(b>1.-pval[pind1[k][1-j]][1-j]) b=1.-pval[pind1[k][1-j]][1-j];
		}
		if(nt==0){
		if(j==0) nstep=nstep1;
		if(j==1) nstep=nstep2;
		for(k=0;k<ndon;k++) for(l=0;l<2;l++){ 
			pfix2[k][l]=pfix[k][l];
			pval2[k][l]=pval[k][l];
			ptype2[k][l*4+2]=ptype[k][l*4+2];
			ptype2[k][l*4+3]=ptype[k][l*4+3];
		}
		pfix2[i][j]=1;
		ptype2[i][j*4+2]=1;
for(k=0;k<n2loc;k++) if(pind2[k][j]==i){
if(ptype[pind2[k][1-j]][4*(1-j)+3]==1){
printf("Sorry, confidence interval calculation for this parameter\n");
printf("not yet implemented (deposit parameter d%d %d paired with\n",
j+1,i);
printf("d%d %d which is paired with another d%d parameter that is\n",
2-j,pind2[k][1-j],j+1);
printf("set to a fixed value.\n");
nt=1;
}
ptype2[pind2[k][1-j]][4*(1-j)+3]=1;
}
for(k=0;k<n1loc;k++) if(pind1[k][j]==i){
if(ptype[pind1[k][1-j]][4*(1-j)+3]==1){
printf("Sorry, confidence interval calculation for this parameter\n");
printf("not yet implemented (deposit parameter d%d %d paired with\n",
j+1,i);
printf("d%d %d which is paired with another d%d parameter that is\n",
2-j,pind1[k][1-j],j+1);
printf("set to a fixed value.\n");
nt=1;
}
ptype2[pind1[k][1-j]][4*(1-j)+3]=1;
}
if(nt==0){
		a=maxvals[ndon+2*i+j];
		if(a>b) upp=a; else{
			for(k=0;k<nstep;k++){
				pval2[i][j]=(a+b)/2.;
				for(l=0;l<ndon;l++) 
				for(ll=0;ll<2;ll++) 
				pstart2[l][ll]=maxvals[ndon+2*l+ll];
for(l=0;l<n2loc;l++) if(pind2[l][j]==i && 
pstart2[pind2[l][1-j]][1-j]+pval2[i][j]>1.){
	if(j==0) pstart2[pind2[l][1-j]][1-j]=.2*(1.-pval2[i][j]);
	if(j==1) pstart2[pind2[l][1-j]][1-j]=.95*(1.-pval2[i][j]);
}
for(l=0;l<n1loc;l++) if(pind1[l][j]==i && 
pstart2[pind1[l][1-j]][1-j]+pval2[i][j]>1.){
        if(j==0) pstart2[pind1[l][1-j]][1-j]=.2*(1.-pval2[i][j]); 
        if(j==1) pstart2[pind1[l][1-j]][1-j]=.95*(1.-pval2[i][j]); 
}        
                        spermanalysisrfssdem(ndon,nal,n2loc,
twoloc,sind2,pind2,amal2,coal2,twolocmp,n1loc,oneloc,sind1,pind1,amal1,
coal1,onelocrf,x,y,z,code,output,prpheno0,prpheno1,prpheno2,mpp,denom,
tmp,tmpt,ltmp,s,p,am,co,pin,amin,
coin,prpheno11,prpheno21,x1,y1,z1,code1,loglike,exp2,exp1,sden,
pden,amden,coden,maxloglike,sfix2,pfix2,amfix2,cofix2,sval2,pval2,
amval2,coval2,sstart2,pstart2,amstart2,costart2,dconv,niter,1,ptype2,
dtmp);

				/*for(l=0;l<ndon;l++) 
				printf("%7e ",s[l]);
				printf("\n");
				for(m=0;m<ndon;m++){
					for(l=0;l<2;l++)
					printf("%7e ",p[m][l]);
					printf("\n");
				}
				for(l=0;l<nal;l++) 
				printf("%7e ",am[l]);
				printf("\n");
				for(l=0;l<nal;l++) 
				printf("%7e ",co[l]);
				printf("\n");
				printf("%7e\n",*maxloglike);*/

				if(*maxloglike>
				maxvals[3*ndon+2*nal]- *quan/2.) 
				a=(a+b)/2.;
				else b=(a+b)/2.;
				/*printf("%7e %7e\n",a,b);*/
			fflush(stdout);
			}
			upp=b;
		}
		/*printf("%7e\n",upp);*/

		a=min;
		b=maxvals[ndon+2*i+j];
		if(a>b) low=b; else{
			if(j==1) nstep=nstep3;
			for(k=0;k<nstep;k++){
				pval2[i][j]=(a+b)/2.;
                        spermanalysisrfssdem(ndon,nal,n2loc,
twoloc,sind2,pind2,amal2,coal2,twolocmp,n1loc,oneloc,sind1,pind1,amal1,
coal1,onelocrf,x,y,z,code,output,prpheno0,prpheno1,prpheno2,mpp,denom,
tmp,tmpt,ltmp,s,p,am,co,pin,amin,
coin,prpheno11,prpheno21,x1,y1,z1,code1,loglike,exp2,exp1,sden,
pden,amden,coden,maxloglike,sfix2,pfix2,amfix2,cofix2,sval2,pval2,
amval2,coval2,sstart2,pstart2,amstart2,costart2,dconv,niter,1,ptype2,
dtmp);

				/*for(l=0;l<ndon;l++)
				printf("%7e ",s[l]);
				printf("\n");
				for(m=0;m<ndon;m++){
					for(l=0;l<2;l++)
					printf("%7e ",p[m][l]);
					printf("\n");   
				}
				for(l=0;l<nal;l++) 
				printf("%7e ",am[l]);
				printf("\n");
				for(l=0;l<nal;l++) 
				printf("%7e ",co[l]);
				printf("\n");
				printf("%7e\n",*maxloglike);*/

				if(*maxloglike>
				maxvals[3*ndon+2*nal]- *quan/2.) 
				b=(a+b)/2.;
				else a=(a+b)/2.;
				/*printf("%7e %7e\n",a,b);
			fflush(stdout);*/
			}
			low=a;
		}
		/*printf("%7e\n",low);*/

			printf("d%d %3d %15f %9f %7f\n",j+1,i,
			maxvals[ndon+2*i+j],low,upp);
			fflush(stdout);
}
		}
	}
}

for(k=0;k<ndon;k++) for(l=0;l<2;l++){
	pfix2[k][l]=pfix[k][l];
	pval2[k][l]=pval[k][l];
}

printf("\nAmplification parameters\n\n");
printf("Parameter     Estimate     Interval\n");

/*for(i=0;i<nal;i++){*/
for(i=0;i<0;i++){
	for(j=0;j<nal;j++){ 
		amfix2[j]=amfix[j];
		amval2[j]=amval[j];
	}
	nt=1;
	for(j=0;j<n2loc;j++) for(l=0;l<4;l++) if(amal2[j][l]==i &&
	amfix[i]==0) nt=0;
	for(j=0;j<n1loc;j++) for(l=0;l<2;l++) if(amal1[j][l]==i &&
	amfix[i]==0) nt=0;
	if(nt==0){
		amfix2[i]=1;
		a=maxvals[3*ndon+i];	
                b=max;
                if(a>b) upp=a; else{
                        for(j=0;j<nstep1;j++){
				amval2[i]=(a+b)/2.;
                        spermanalysisrfssdem(ndon,nal,n2loc,
twoloc,sind2,pind2,amal2,coal2,twolocmp,n1loc,oneloc,sind1,pind1,amal1,
coal1,onelocrf,x,y,z,code,output,prpheno0,prpheno1,prpheno2,mpp,denom,
tmp,tmpt,ltmp,s,p,am,co,pin,amin,
coin,prpheno11,prpheno21,x1,y1,z1,code1,loglike,exp2,exp1,sden,
pden,amden,coden,maxloglike,sfix2,pfix2,amfix2,cofix2,sval2,pval2,
amval2,coval2,sstart2,pstart2,amstart2,costart2,dconv,niter,1,ptype,
dtmp);

                                /*for(k=0;k<ndon;k++) printf("%7e ",s[i]);
                                printf("\n");
                                for(k=0;k<ndon;k++){
                                        for(l=0;l<2;l++)
                                        printf("%7e ",p[k][l]);
                                        printf("\n");
                                }
                                for(k=0;k<nal;k++) printf("%7e ",am[k]);
                                printf("\n");
                                for(k=0;k<nal;k++) printf("%7e ",co[k]);
                                printf("\n");
                                printf("%7e\n",*maxloglike);*/
 
                                if(*maxloglike>maxvals[3*ndon+2*nal]-
                                *quan/2.) a=(a+b)/2.;
                                else b=(a+b)/2.;
                                /*printf("%7e %7e\n",a,b);
			fflush(stdout);*/
                        }
                        upp=b;
                }
                /*printf("%7e\n",upp);*/
 
                a=min;
                b=maxvals[3*ndon+i];
                if(a>b) low=b; else{
                        for(j=0;j<nstep1;j++){
				amval2[i]=(a+b)/2.;
                        spermanalysisrfssdem(ndon,nal,n2loc,
twoloc,sind2,pind2,amal2,coal2,twolocmp,n1loc,oneloc,sind1,pind1,amal1,
coal1,onelocrf,x,y,z,code,output,prpheno0,prpheno1,prpheno2,mpp,denom,
tmp,tmpt,ltmp,s,p,am,co,pin,amin,
coin,prpheno11,prpheno21,x1,y1,z1,code1,loglike,exp2,exp1,sden,
pden,amden,coden,maxloglike,sfix2,pfix2,amfix2,cofix2,sval2,pval2,
amval2,coval2,sstart2,pstart2,amstart2,costart2,dconv,niter,1,ptype,
dtmp);

                                /*for(k=0;k<ndon;k++) printf("%7e ",s[i]);
                                printf("\n"); 
                                for(k=0;k<ndon;k++){ 
                                        for(l=0;l<2;l++) 
                                        printf("%7e ",p[k][l]); 
                                        printf("\n"); 
                                } 
                                for(k=0;k<nal;k++) printf("%7e ",am[k]);

                                printf("\n"); 
                                for(k=0;k<nal;k++) printf("%7e ",co[k]);

                                printf("\n"); 
                                printf("%7e\n",*maxloglike); */
  
                                if(*maxloglike>maxvals[3*ndon+2*nal]- 
                                *quan/2.) b=(a+b)/2.;
                                else a=(a+b)/2.;
                                /*printf("%7e %7e\n",a,b);
				fflush(stdout);*/
                        }
                        low=a;
                }
                /*printf("%7e\n",low);*/
 
printf("a %4d %15f %9f %7f\n",i,maxvals[3*ndon+i],low,upp);
		fflush(stdout);
        }
}

for(j=0;j<nal;j++){ 
	amfix2[j]=amfix[j];
	amval2[j]=amval[j];
}

printf("\nContamination parameters\n\n");
printf("Parameter     Estimate     Interval\n");

/*for(i=0;i<nal;i++){*/
for(i=0;i<0;i++){
        for(j=0;j<nal;j++){ 
		cofix2[j]=cofix[j];
		coval2[j]=coval[j];
	}
        nt=1;
        for(j=0;j<n2loc;j++) for(l=0;l<4;l++) if(coal2[j][l]==i &&
	cofix[i]==0) nt=0;
        for(j=0;j<n1loc;j++) for(l=0;l<2;l++) if(coal1[j][l]==i &&
	cofix[i]==0) nt=0;
        if(nt==0){
                cofix2[i]=1;
                a=maxvals[3*ndon+nal+i];
                b=max;
                if(a>b) upp=a; else{
                        for(j=0;j<nstep2;j++){
                                coval2[i]=(a+b)/2.;
                        spermanalysisrfssdem(ndon,nal,n2loc,
twoloc,sind2,pind2,amal2,coal2,twolocmp,n1loc,oneloc,sind1,pind1,amal1,
coal1,onelocrf,x,y,z,code,output,prpheno0,prpheno1,prpheno2,mpp,denom,
tmp,tmpt,ltmp,s,p,am,co,pin,amin,
coin,prpheno11,prpheno21,x1,y1,z1,code1,loglike,exp2,exp1,sden,
pden,amden,coden,maxloglike,sfix2,pfix2,amfix2,cofix2,sval2,pval2,
amval2,coval2,sstart2,pstart2,amstart2,costart2,dconv,niter,1,ptype,
dtmp);

                                /*for(k=0;k<ndon;k++) printf("%7e ",s[i]);
                                printf("\n");
                                for(k=0;k<ndon;k++){
                                        for(l=0;l<2;l++)
                                        printf("%7e ",p[k][l]);
                                        printf("\n");
                                }
                                for(k=0;k<nal;k++) printf("%7e ",am[k]);
                                printf("\n");
                                for(k=0;k<nal;k++) printf("%7e ",co[k]);
                                printf("\n");
                                printf("%7e\n",*maxloglike);*/
 
                                if(*maxloglike>maxvals[3*ndon+2*nal]-
                                *quan/2.) a=(a+b)/2.;
                                else b=(a+b)/2.;
                                /*printf("%7e %7e\n",a,b);
				fflush(stdout);*/
                        }
                        upp=b;
                }
                /*printf("%7e\n",upp);*/
 
                a=min;
                b=maxvals[3*ndon+nal+i];
                if(a>b) low=b; else{
                        for(j=0;j<nstep3;j++){
                                coval2[i]=(a+b)/2.;
                        spermanalysisrfssdem(ndon,nal,n2loc,
twoloc,sind2,pind2,amal2,coal2,twolocmp,n1loc,oneloc,sind1,pind1,amal1,
coal1,onelocrf,x,y,z,code,output,prpheno0,prpheno1,prpheno2,mpp,denom,
tmp,tmpt,ltmp,s,p,am,co,pin,amin,
coin,prpheno11,prpheno21,x1,y1,z1,code1,loglike,exp2,exp1,sden,
pden,amden,coden,maxloglike,sfix2,pfix2,amfix2,cofix2,sval2,pval2,
amval2,coval2,sstart2,pstart2,amstart2,costart2,dconv,niter,1,ptype,
dtmp);

                                /*for(k=0;k<ndon;k++) printf("%7e ",s[i]);
                                printf("\n");
                                for(k=0;k<ndon;k++){
                                        for(l=0;l<2;l++)
                                        printf("%7e ",p[k][l]);
                                        printf("\n");
                                }
                                for(k=0;k<nal;k++) printf("%7e ",am[k]);
 
                                printf("\n");
                                for(k=0;k<nal;k++) printf("%7e ",co[k]);
 
                                printf("\n");
                                printf("%7e\n",*maxloglike);*/
 
                                if(*maxloglike>maxvals[3*ndon+2*nal]-
                                *quan/2.) b=(a+b)/2.;
                                else a=(a+b)/2.;
                                /*printf("%7e %7e\n",a,b);
				fflush(stdout);*/
                        }
                        low=a;
                }
                /*printf("%7e\n",low);*/
 
printf("c %4d %15f %9f %7f\n",i,maxvals[3*ndon+nal+i],low,upp);
		fflush(stdout);
        }
}

}
#include <math.h>
#include <stdio.h>

/* maxvals[3*nd(n+2*nal+1] is the vector of parameter values that 
maximize the
likelihood, with last entry the maximized likelihood, quan is
qchisq(.95,1) for a 95% confidence interval, qchisq(.68,1) to get
standard errors */

void spermanalysisfssdioemcisonly(int ndon,int nal,int n2loc,double **twoloc,
int *sind2,int **pind2,int **amal2,int **coal2,double **twolocmp,
int n1loc,double **oneloc,int *sind1,int **pind1,int **amal1,int **coal1,
double *onelocrf,int **x,int **y,int **z,int **code,double *output,
double *prpheno0,double **prpheno1,double **prpheno2,double *mpp,
double *denom,double *tmp,double *tmpt,double *ltmp,double *s,
double **p,double *am,double *co,double *pin,double *amin,double *coin,
double **prpheno11,double **prpheno21,int **x1,int **y1,int **z1,
int **code1,double *loglike,double **exp2,double **exp1,double *sden,
double **pden,double *amden,double *coden,double *maxloglike,int *sfix,
int **pfix,int *amfix,int *cofix,double *sval,double **pval,
double *amval,double *coval,double *maxvals,double *quan,
int *sfix2,int **pfix2,int *amfix2,int *cofix2,double *sval2,
double **pval2,double *amval2,double *coval2,double *sstart,double
**pstart,double *amstart,double *costart,double *dconv,int *niter,int
nstep1,double *sstart2,double **pstart2,double *amstart2,
double *costart2,int **ptype,double *dtmp)
{
int i,j,nt;
double eps,min,max,a,b,upp,low;

void spermanalysisrfssdem(int ndon,int nal,int n2loc,double **twoloc,
int *sind2,int **pind2,int **amal2,int **coal2,double **twolocmp,
int n1loc,double **oneloc,int *sind1,int **pind1,int **amal1,int **coal1,
double *onelocrf,int **x,int **y,int **z,int **code,double *output,
double *prpheno0,double **prpheno1,double **prpheno2,double *mpp,
double *denom,double *tmp,double *tmpt,double *ltmp,double *s,
double **p,double *am,double *co,double *pin,double *amin,double *coin,
double **prpheno11,double **prpheno21,int **x1,int **y1,int **z1,
int **code1,double *loglike,double **exp2,double **exp1,double *sden,
double **pden,double *amden,double *coden,double *maxloglike,int *sfix,
int **pfix,int *amfix,int *cofix,double *sval,double **pval,
double *amval,double *coval,double *sstart,double **pstart,double
*amstart,double *costart,double *dconv,int *niter,int simind,int
**ptype,double *dtmp);

eps=pow(10.,-10.);
min=eps;
max=1-eps;

for(i=0;i<ndon;i++){
sfix2[i]=sfix[i];
sval2[i]=sval[i];
pfix2[i][0]=pfix[i][0];
pval2[i][0]=pval[i][0];
pfix2[i][1]=pfix[i][1];
pval2[i][1]=pval[i][1];
}

for(i=0;i<nal;i++){
amfix2[i]=amfix[i];
amval2[i]=amval[i];
cofix2[i]=cofix[i];
coval2[i]=coval[i];
}

for(i=0;i<ndon;i++){
	sstart2[i]=maxvals[i];
	for(j=0;j<2;j++) pstart2[i][j]=maxvals[ndon+2*i+j];
}

for(i=0;i<nal;i++){
	amstart2[i]=maxvals[3*ndon+i];
	costart2[i]=maxvals[3*ndon+nal+i];
}

printf("Confidence intervals for segregation parameters\n\n");
printf("Parameter     Estimate     Interval\n");

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

for(j=0;j<ndon;j++){ 
	sfix2[j]=sfix[j];
	sval2[j]=sval[j];
}
nt=1;
for(j=0;j<n2loc;j++) if(sind2[j]==i && sfix[i]==0) nt=0;
for(j=0;j<n1loc;j++) if(sind1[j]==i && sfix[i]==0) nt=0;
if(nt==0){
	sfix2[i]=1;
	a=maxvals[i];
	b=max;
	if(a>b) upp=a; else{
		for(j=0;j<nstep1;j++){
			sval2[i]=(a+b)/2.;
			spermanalysisrfssdem(ndon,nal,n2loc,
twoloc,sind2,pind2,amal2,coal2,twolocmp,n1loc,oneloc,sind1,pind1,amal1,
coal1,onelocrf,x,y,z,code,output,prpheno0,prpheno1,prpheno2,mpp,denom,
tmp,tmpt,ltmp,s,p,am,co,pin,amin,
coin,prpheno11,prpheno21,x1,y1,z1,code1,loglike,exp2,exp1,sden,
pden,amden,coden,maxloglike,sfix2,pfix2,amfix2,cofix2,sval2,pval2,
amval2,coval2,sstart2,pstart2,amstart2,costart2,dconv,niter,1,ptype,
dtmp);

			/*for(k=0;k<ndon;k++) printf("%7e ",s[k]);
			printf("\n");
			for(k=0;k<ndon;k++){
				for(l=0;l<2;l++) 
				printf("%7e ",p[k][l]);
				printf("\n");
			}
			for(k=0;k<nal;k++) printf("%7e ",am[k]);
			printf("\n");
			for(k=0;k<nal;k++) printf("%7e ",co[k]);
			printf("\n");
			printf("%7e\n",*maxloglike);*/

			if(*maxloglike>maxvals[3*ndon+2*nal]-
			*quan/2.) a=(a+b)/2.;
			else b=(a+b)/2.;
			/*printf("%7e %7e\n",a,b);
			fflush(stdout);*/
		}
		upp=b;
	}
	/*printf("%7e\n",upp);*/

	a=min;
	b=maxvals[i];
	if(a>b) low=b; else{
		for(j=0;j<nstep1;j++){
			sval2[i]=(a+b)/2.;
                        spermanalysisrfssdem(ndon,nal,n2loc,
twoloc,sind2,pind2,amal2,coal2,twolocmp,n1loc,oneloc,sind1,pind1,amal1,
coal1,onelocrf,x,y,z,code,output,prpheno0,prpheno1,prpheno2,mpp,denom,
tmp,tmpt,ltmp,s,p,am,co,pin,amin,
coin,prpheno11,prpheno21,x1,y1,z1,code1,loglike,exp2,exp1,sden,
pden,amden,coden,maxloglike,sfix2,pfix2,amfix2,cofix2,sval2,pval2,
amval2,coval2,sstart2,pstart2,amstart2,costart2,dconv,niter,1,ptype,
dtmp);

			/*for(k=0;k<ndon;k++) printf("%7e ",s[i]);
			printf("\n");
			for(k=0;k<ndon;k++){
				for(l=0;l<2;l++)
				printf("%7e ",p[k][l]);
				printf("\n");
			}
			for(k=0;k<nal;k++) printf("%7e ",am[k]);
			printf("\n");
			for(k=0;k<nal;k++) printf("%7e ",co[k]);
			printf("\n");
			printf("%7e\n",*maxloglike);*/

			if(*maxloglike>maxvals[3*ndon+2*nal]-
			*quan/2.) b=(a+b)/2.;
			else a=(a+b)/2.;
			/*printf("%7e %7e\n",a,b); 
			fflush(stdout);*/
		} 
		low=a;
	}
	/*printf("%7e\n",low);*/

printf("s %4d %15f %9f %7f\n",i,maxvals[i],low,upp);
	fflush(stdout);
}
}

}
/* This routine calculates the chi-square quantile (1 df) associated
with the probability cilevel, i.e. this is the inverse cdf of a
chi-square on 1 df. */

#include <math.h>
#include <stdio.h>

void chisqquan(double cilevel,double *quan,double *tmp){

int i,j;
double gammacdf(double df, double z,double *tmp);

for(i=1;gammacdf(.5,(double)i/2.,tmp+2)<cilevel;i++);
*tmp=(double)(i-1);
tmp[1]=(double)i;
for(j=0;j<24;j++){
	*quan=(*tmp+tmp[1])/2.;
	if(gammacdf(.5,*quan/2.,tmp+2)<cilevel) *tmp=*quan;
	else tmp[1]=*quan;
}
*quan=tmp[1];
}

/* This routine calculates the log of the gamma function */

double lngamma(double z,double *tmp){

*tmp=(z+.5)*log(z+5.5)-z-5.5;
*tmp+=log(2.5066282746310005*(1.000000000190015+76.18009172947146/
(z+1)-86.50532032941677/(z+2)+24.01409824083091/(z+3)-
1.231739572450155/(z+4)+0.1208650973866179e-2/(z+5)-
0.5395239384953e-5/(z+6))/z);
return *tmp;
}

/* This routine calculates the series representation of the incomplete
gamma function */

void igseries(double *ans,double df,double z,double *tmp){

double lngamma(double z,double *tmp);
int i;

if(z<=0.) *ans=0.;
else{
	tmp[2]=1./df;
	tmp[1]=df;
	*tmp=1./df;
	for(i=0;i<200;i++){
		tmp[1]++;
		*tmp*=(z/tmp[1]);
		tmp[2]+=*tmp;
	}
	*ans=exp(-z+df*log(z)-lngamma(df,tmp))*tmp[2];
}
}

/* This routine calculates the continued fraction representation of the
incomplete gamma function */

void igcontfrac(double *ans,double df,double z,double *tmp){

double lngamma(double z,double *tmp);
int i;

tmp[6]=pow(10.,-30.);
tmp[1]=z+1.-df;
tmp[2]=1./tmp[6];
tmp[3]=1./tmp[1];
tmp[5]=tmp[3];
for(i=1;i<=200;i++){
	*tmp=-(double)i*((double)i-df);
	tmp[1]+=2.;
	tmp[3]=*tmp*tmp[3]+tmp[1];
	if(fabs(tmp[3])<tmp[6]) tmp[3]=tmp[6];
	tmp[2]=tmp[1]+*tmp/tmp[2];
	if(fabs(tmp[2])< tmp[6]) tmp[2]=tmp[6];
	tmp[3]=1./tmp[3];
	tmp[4]=tmp[3]*tmp[2];
	tmp[5]*=tmp[4];
	if(fabs(tmp[4]-1.)<pow(10.,-8.)) break;
}
*ans=exp(-z+df*log(z)-lngamma(df,tmp))*tmp[5];
}

/* This routines calculates the cdf of the gamma density */

double gammacdf(double df,double z,double *tmp){

void igcontfrac(double *ans,double df,double z,double *tmp);
void igseries(double *ans,double df,double z,double *tmp);

if(z<(df+1.)){
	igseries(tmp+7,df,z,tmp);
	return tmp[7];
} else{
	igcontfrac(tmp+7,df,z,tmp);
	return 1.-tmp[7];
}
}
/*Must use with seed.c*/

#include <math.h>
#include <stdio.h>

double randomnum(long *seed,long *ab,long *cd)
{
    int j;
    long k;
    double temp;

    if (*seed <= 0 || !*cd){
	if (-(*seed) < 1) *seed=1;
	else *seed = -(*seed);
	for (j=32+7; j >=0; j--){
	    k = (*seed)/127773;
	    *seed = 16807*(*seed-k*127773)-2836*k;
	    if (*seed < 0) *seed += 2147483647;
	    if (j<32) ab[j] = *seed;
	}
	*cd = ab[0];
    }
    k = (*seed)/127773;
    *seed = 16807*(*seed-k*127773)-2836*k;
    if (*seed <0) *seed += 2147483647;
    j = *cd/(1+1073741823/16);
    *cd = ab[j];
    ab[j] = *seed;
    if ((temp = (1./2147483647)**cd) > 1.-1.2e-7) temp = (1.-1.2e-7);
    return temp;
}