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

/* 

   Program "hwmet.c"          (12/06/97)

   Updated for ANSI C; 18/07/00   - to compile type  cc hwmet.c -lm

   Updated to output summary statistics; 09/09/02


   Karen L. Ayres               
   University of Reading 
   e-mail : k.l.ayres@reading.ac.uk 

   Metropolis-Hastings algorithm for obtaining a sample from the
   posterior distribution of the allele frequencies and the inbreeding
   coefficient f, for the Hardy-Weinberg inbreeding model
   (e.g. Malecot, 1969):

   p_ii=p_i(f+(1-f)p_i)
   p_ij=2p_ip_j(1-f)

   The method is detailed in Ayres & Balding, Heredity 1998

   Allows for negative values of f, provided that the condition
   p_ii>=0 is satisfied - for only non-negative values of f, the
   program "hwmet2.c" is more efficient.

   Reads from input file "hwmet.txt", which is a text file
   WHOSE FIRST LINE CONTAINS ONLY AN INTEGER WHICH IS THE NUMBER
   OF ALLELES IN THE RESULTING DATA SET - this is stored in "nall".
   The following lines in the input file should consist of a "nall" by
   "nall" matrix of observed genotype counts - counts for homozygous
   genotypes are on the diagonal, heterozygous genotypes are off the
   diagonal e.g.
   3
   12 34 21
   0  7  11
   0  0   8

   The program gentable.c will create this matrix.

   Outputs "f" parameter values to file "hwmet.out", "p" vectors
   to "hwmetp.out", and information about the implementation to
   "hwmet.stats" (for many alleles, the number of iterates output
   to files (nout) should be reduced (e.g. 2000), since "hwmetp.out"
   will require a large amount of storage space). Note the all files
   *.out or *.stats are simple ASCII text files, but I have used 
   extensions other than txt to indicate what is in the file (out=
   output) - they can all be opened in notepad, wordpad, or any other
   text editor.

   The observed genotype counts are stored in "n".

   The {p_i} have a multivariate uniform prior, and
   f has a Beta(a,b) prior conditional on the vector p (since the
   permissable range of f is dependent on p) on the interval (A,1),
   where A=-pmin/(1-pmin), and p_min is the smallest value in p
   (assuming all p[i] are non-zero), with parameters

   a = 1-A.(b-1)
   b = g_t(p_min)

   and where g_t(*) is a polynomial function of p_min, whose form is
   determined by the value "t", which represents where the right-hand
   tail of the prior density curve is to become approximately 0. This
   program currently caters for three values of "t" - 0.1, 0.2, and 0.4,
   which are likely to be useful for many human populations since they
   give most support to values of f close to zero. The conditional priors
   all have mode at zero, which may be reasonable a priori. (If these
   priors are not reasonable, it may be preferable to apply the program
   "hwmet2.c" which restricts f to non-negative values and implements
   a standard Beta prior.)

   The values of "t" in the program are coded as follows:

   for t=0.1  : set btail = 1
   for t=0.2  : set btail = 2
   for t=0.4  : set btail = 4

   for a uniform (i.e. non-informative) prior  : set btail = 0
   (Note: the program "fprior.c" can be used to inspect the form of the
   marginal prior of f for each of the above values of "t").

   The parameter "btail" is defined at the beginning of the program,
   and should be set to one of the above integers accordingly.
 
   A Metropolis-Hastings algorithm is run for "numit" iterations,
   the first "discard" of which are discarded as the burn-in.

   Subsequently, every "kp"th iteration is output.

   Note that "numit" is determined by the values of "kp", "discard"
   and "nout", where "nout" is the number of values desired to be
   output

   The limiting bound "epsln" (>0) can be changed as necessary, in order 
   to obtain reasonable acceptance rates.

   Density curves can be produced within a suitable statistical computer
   package by applying density estimation to the output from this
   program.

   Note that, for a population subdivided into smaller subpopulations,
   if the input file "hwmet.txt" is that for a subpopulation "s", and
   allele frequencies p[i] are fixed equal to estimates for the total
   population, the algorithm implemented here will approximate the
   posterior distribution, for the locus in question, of the correlation
   of genes within an individual in "s" relative to the total population. 
*/

#define btail 4  /* determines right-hand tail of conditional beta prior of F */
#define epsln 0.02  /* limit on amount to update p's by (< 0.5 and > 0) */

#define maxall 35  /* maximum number of alleles allowed at a locus */
#define kp (maxall*4)   /* interval between successive output values */
#define nout 5000L     /* number of output values required */
#define numit (nout*kp) /* number of iterations of alg. after burn-n */

#define discard ((long)(0.05*numit))   /* number of values to initially discard */

#define acrt (numit/10)  /* interval between output of acceptance rates */
#define dim1 ((maxall*(maxall+1))/2)

typedef double pvec[maxall];

int     i, j, k, nsam=0, nall, ngen, p1, p2, x[maxall][maxall], n[dim1], nnn;
long    l1, jmp, ix, iy, iz;
double  y, b1, b2, b3, b4, a1, a2, a3, a4, r, pmin, newpmin, f, newf, gsum;
double  epsln2, genf[dim1], z, newz, fst[nout], tmpd, u, mnf;
pvec    p, newp;
FILE    *inp, *out1, *out2;

/* PROTOTYPES */

double WichHill();
void inpdat();
double loggf(double xx);
double lbetad(double fs, double pms);
double rstgam(double alpha1);
double rbetaf();
double llkhd(double fs,pvec ps);

/* FUNCTIONS */

double WichHill()
{
/* Wichmann & Hill random number generator (Algorithm AS183 in Applied
   Statistics 31 (1982))
*/

ix=(171*ix)%30269;
iy=(172*iy)%30307;
iz=(170*iz)%30323;

return fmod((ix/30269.0)+(iy/30307.0)+(iz/30323.0),1.0);
}

void inpdat()
{
int i,j,k;
double gsum, mnf, blb, bpa, bpb;

out1=fopen("hwmet.stats","w");
inp=fopen("hwmet.txt","r");
fscanf(inp,"%d",&nall);

/* Check for invalid input file */

if (nall>maxall){
	fprintf(out1,"ERROR - TOO MANY ALLELES - CHANGE MAXALL\n");
	exit(1);
}

ngen=(nall*(nall+1))/2;
for (i=0;i<nall;++i){p[i]=0;}
for (i=0;i<nall;++i){
	for (j=0,k=0;j<nall;++j){
		fscanf(inp,"%d",&x[i][j]);
		nsam+= x[i][j];
            p[i]+= x[i][j];
            p[j]+= x[i][j];
	}
}

for (i=0;i<nall;++i){
	if (p[i]==0){
		fprintf(out1,"WARNING - unobserved allele %d in the data set will force lower bound for f to be approx. 0\n",i+1);
	}
}

/* observed genotype counts, using genotype identifier for ij: i+(j*(j+1))/2
   where i<=j (modified for i=0,...,k-1 from that given for i=1,...,k in
   Zaykin, Zhivotovsky & Weir (1995) - Genetica 96)
*/

for (i=0;i<nall;++i){
	for (j=i;j<nall;++j){
		n[i+(j*(j+1))/2]=x[i][j];
		if (i!=j){
			n[i+(j*(j+1))/2]+= x[j][i];
                  x[i][j]+= x[j][i];
                  x[j][i]=0;
		}
	}
}

/* set value of epsln2, the limiting bound for the updating distribution
   for f --- setting this > [(nall^2).epsln]/[(nall-1).(nall-1-nall.epsln)]
   ensures valid convergence of the Markov chain
*/

if (epsln>=0.5){
	fprintf(out1,"ERROR - LIMITING BOUND FOR P INVALID - CHANGE EPSLN (<0.5)\n");
	exit(1);
}

epsln2=epsln*nall*nall/((nall-1.0)*(nall-1.0-nall*epsln))+10e-10;

/* output statistics to file */

fprintf(out1,"LOG FILE FOR RUN OF HWMET ON THE FOLLOWING DATA\n\n");

fprintf(out1,"****** Adjusted data table (letting heterozygotes AB=BA)******\n");

for (i=0;i<nall;++i){
	for (j=0,k=0;j<nall;++j){
            fprintf(out1,"%5d ",x[i][j]);
	}
	fprintf(out1,"\n");
}

fprintf(out1,"\n\nSample size = %d\n\n",nsam);

fprintf(out1,"****** Settings for MCMC algorithm ******\n\n");
fprintf(out1,"Burn-in length (iterations): %ld\n",discard);
fprintf(out1,"Thinning interval: %ld\n",kp);
fprintf(out1,"Number of iterations (post burn-in): %ld\n",numit);
fprintf(out1,"Outputting %ld values from posterior\n",nout);
fprintf(out1,"Random number seeds (generated via system time function) are:\n");
fprintf(out1,"     ix= %ld, iy= %ld, iz= %ld\n\n",ix,iy,iz);

fprintf(out1,"Limiting bound for proposal distribution of p_i ... %6.3f\n",epsln);
fprintf(out1,"Limiting bound for proposal distribution of f   ... %6.4f\n\n",epsln2);

fprintf(out1,"****** Priors used ******\n\n");

fprintf(out1,"Using multivariate uniform prior for allele frequency vector\n");

if (btail==0){
	fprintf(out1,"Conditional prior for f - uniform\n");
}
else if (btail==1){
	fprintf(out1,"Conditional prior for f - beta with right-hand tail out to 0.1 approx.\n");
}
else if (btail==2){
	fprintf(out1,"Conditional prior for f - beta with right-hand tail out to 0.2 approx.\n");
}
else if (btail==4){
	fprintf(out1,"Conditional prior for f - beta with right-hand tail out to 0.4 approx.\n");
}
else{
	fprintf(out1,"ERROR - INVALID TAIL PARAMETER FOR BETA PRIOR - CHANGE BTAIL (0,1,2,4)\n");
	exit(1);
}

fprintf(out1,"\n** Summary stats, based on %ld simulated values, for prior of f **\n\n",nout);

fprintf(out1,"---------------------------------------------------------------------------------------------------\n");
fprintf(out1,"|         | Mean     |2.5 %c-ile |  5 %c-ile | 25 %c-ile |  median  | 75 %c-ile | 95 %c-ile |97.5 %c-ile|\n",37,37,37,37,37,37);
fprintf(out1,"|-------------------------------------------------------------------------------------------------|\n");

mnf=0.0;
for (l1=0;l1<nout;++l1){
	gsum=0;
	for (i=0;i<nall;++i){
		p[i]= (-log(WichHill()));
		gsum+= p[i];
	}
	for (i=0;i<nall;++i){
		p[i]= p[i]/gsum;
	}  
	pmin=1.0/nall;
	for (i=0;i<nall;++i){
		if (p[i]<pmin){
			pmin=p[i];
		}
	}
	blb= -pmin/(1.0-pmin);
	if (btail==0){
		bpb=1.0;
	}
	else if (btail==1){
		bpb=61.37+1482.08*pmin-1643.18*pow(pmin,2.0)+41521.78*pow(pmin,3.0)-298223.0*pow(pmin,4.0)+723583.9*pow(pmin,5.0)-579512.0*pow(pmin,6.0);
	}
	else if (btail==2){
		bpb=23.61+241.3*pmin-1048.65*pow(pmin,2.0)+8944.7*pow(pmin,3.0)-20820.2*pow(pmin,4.0)+14112.95*pow(pmin,5.0);
	}
	else if (btail==4){
		bpb=11.15+29.75*pmin+773.84*pow(pmin,2.0)-4409.76*pow(pmin,3.0)+9156.72*pow(pmin,4.0)-6515.75*pow(pmin,5.0);
	}
	else {
		fprintf(out1,"ERROR - BTAIL MUST BE = 0, 1, 2 or 4\n");
		exit(1);
	}
	bpa= 1.0-blb*(bpb-1.0);
	f=rbetaf();
	mnf+= f;
	fst[l1]=f;
}
mnf/= nout;

for (l1=0;l1<nout-1;++l1){
   for (jmp=l1+1;jmp<nout;++jmp){
      if (fst[l1]>fst[jmp]){
         tmpd=fst[l1];
         fst[l1]=fst[jmp];
         fst[jmp]=tmpd;
      }
   }
}

/* CALCULATE PERCENTILES BY SETTING i'th ORDER STATISTIC EQUAL TO
   (i-0.5)/n AND USING LINEAR INTERPOLATION
*/

if (mnf>0){
	fprintf(out1,"|Prior    | %9.7f|",mnf);
}
else {
	fprintf(out1,"|Prior    |%9.7f|",mnf);
}

/* 2.5th */

u=0.025*nout+0.5;
k=(int)(u);
if (k==u){
	if (fst[k]>0){
	   fprintf(out1," %9.7f|",fst[k]);
	}
	else {
	   fprintf(out1,"%9.7f|",fst[k]);
	}
}
else {
	if ((fst[k]+fmod(u,1.0)*(fst[k+1]-fst[k]))>0){
	   fprintf(out1," %9.7f|",fst[k]+fmod(u,1.0)*(fst[k+1]-fst[k]));
	}
	else {
	   fprintf(out1,"%9.7f|",fst[k]+fmod(u,1.0)*(fst[k+1]-fst[k]));
	}
}         

/* 5th */

u=0.05*nout+0.5;
k=(int)(u);
if (k==u){
	if (fst[k]>0){
	   fprintf(out1," %9.7f|",fst[k]);
	}
	else {
	   fprintf(out1,"%9.7f|",fst[k]);
	}
}
else {
	if ((fst[k]+fmod(u,1.0)*(fst[k+1]-fst[k]))>0){
	   fprintf(out1," %9.7f|",fst[k]+fmod(u,1.0)*(fst[k+1]-fst[k]));
	}
	else {
	   fprintf(out1,"%9.7f|",fst[k]+fmod(u,1.0)*(fst[k+1]-fst[k]));
	}
}          

/* 25th */

u=0.25*nout+0.5;
k=(int)(u);
if (k==u){
	if (fst[k]>0){
	   fprintf(out1," %9.7f|",fst[k]);
	}
	else {
	   fprintf(out1,"%9.7f|",fst[k]);
	}
}
else {
	if ((fst[k]+fmod(u,1.0)*(fst[k+1]-fst[k]))>0){
	   fprintf(out1," %9.7f|",fst[k]+fmod(u,1.0)*(fst[k+1]-fst[k]));
	}
	else {
	   fprintf(out1,"%9.7f|",fst[k]+fmod(u,1.0)*(fst[k+1]-fst[k]));
	}
}          

/* 50th */

u=0.5*nout+0.5;
k=(int)(u);
if (k==u){
	if (fst[k]>0){
	   fprintf(out1," %9.7f|",fst[k]);
	}
	else {
	   fprintf(out1,"%9.7f|",fst[k]);
	}
}
else {
	if ((fst[k]+fmod(u,1.0)*(fst[k+1]-fst[k]))>0){
	   fprintf(out1," %9.7f|",fst[k]+fmod(u,1.0)*(fst[k+1]-fst[k]));
	}
	else {
	   fprintf(out1,"%9.7f|",fst[k]+fmod(u,1.0)*(fst[k+1]-fst[k]));
	}
}         

/* 75th */

u=0.75*nout+0.5;
k=(int)(u);
if (k==u){
	if (fst[k]>0){
	   fprintf(out1," %9.7f|",fst[k]);
	}
	else {
	   fprintf(out1,"%9.7f|",fst[k]);
	}
}
else {
	if ((fst[k]+fmod(u,1.0)*(fst[k+1]-fst[k]))>0){
	   fprintf(out1," %9.7f|",fst[k]+fmod(u,1.0)*(fst[k+1]-fst[k]));
	}
	else {
	   fprintf(out1,"%9.7f|",fst[k]+fmod(u,1.0)*(fst[k+1]-fst[k]));
	}
}           

/* 95th */

u=0.95*nout+0.5;
k=(int)(u);
if (k==u){
	if (fst[k]>0){
	   fprintf(out1," %9.7f|",fst[k]);
	}
	else {
	   fprintf(out1,"%9.7f|",fst[k]);
	}
}
else {
	if ((fst[k]+fmod(u,1.0)*(fst[k+1]-fst[k]))>0){
	   fprintf(out1," %9.7f|",fst[k]+fmod(u,1.0)*(fst[k+1]-fst[k]));
	}
	else {
	   fprintf(out1,"%9.7f|",fst[k]+fmod(u,1.0)*(fst[k+1]-fst[k]));
	}
}          

/* 97.5th */

u=0.975*nout+0.5;
k=(int)(u);
if (k==u){
	if (fst[k]>0){
	   fprintf(out1," %9.7f|",fst[k]);
	}
	else {
	   fprintf(out1,"%9.7f|",fst[k]);
	}
}
else {
	if ((fst[k]+fmod(u,1.0)*(fst[k+1]-fst[k]))>0){
	   fprintf(out1," %9.7f|",fst[k]+fmod(u,1.0)*(fst[k+1]-fst[k]));
	}
	else {
	   fprintf(out1,"%9.7f|",fst[k]+fmod(u,1.0)*(fst[k+1]-fst[k]));
	}
}          

fprintf(out1,"---------------------------------------------------------------------------------------------------\n\n");

fprintf(out1,"Acceptance rates for this run of the algorithm ...\n");

fclose(inp);
fflush(out1);
}

double loggf(double xx)
{
/* Adapted from the underlying algorithm in Numerical Recipes in C (1988) 
   [by Press et al., published by Camridge University Press], for 
   calculating the log of the gamma function
*/

double lg=(xx+0.5)*log(xx+5.5)-xx-5.5, const1=1.000000000190015;

const1+= 76.1800917297146/(xx+1.0);
const1+= 24.01409824083091/(xx+3.0);
const1+= 0.001208650973866179/(xx+5.0);
const1-= 86.50532032941677/(xx+2.0);
const1-= 1.2317395572450155/(xx+4.0);
const1-= 0.000005395239384953/(xx+6.0);

lg+= log(2.5066282746310005*const1/xx);

return lg;
}

double lbetad(double fs, double pms)
{
/* evaluate log of the conditional prior density for "fs" */

double blb, bpa, bpb;

blb= -pms/(1.0-pms);
if (btail==0){
	bpb=1.0;
}
else if (btail==1){
	bpb=61.37+1482.08*pms-1643.18*pow(pms,2.0)+41521.78*pow(pms,3.0)-298223.0*pow(pms,4.0)+723583.9*pow(pms,5.0)-579512.0*pow(pms,6.0);
}
else if (btail==2){
	bpb=23.61+241.3*pms-1048.65*pow(pms,2.0)+8944.7*pow(pms,3.0)-20820.2*pow(pms,4.0)+14112.95*pow(pms,5.0);
}
else if (btail==4){
	bpb=11.15+29.75*pms+773.84*pow(pms,2.0)-4409.76*pow(pms,3.0)+9156.72*pow(pms,4.0)-6515.75*pow(pms,5.0);
}
bpa= 1.0-blb*(bpb-1.0);

return (-log(1.0-blb)+loggf(bpa+bpb)-loggf(bpa)-loggf(bpb)+(bpa-1.0)*log((fs-blb)/(1.0-blb))+(bpb-1.0)*log((1.0-fs)/(1.0-blb)));
}

double rstgam(double alpha1)
{
/* Cheng & Feast algorithm (algorithm 3.20 in Ripley, 1988) for generating
   a standard Gamma random variable with shape parameter alpha>1
*/

int     acpt;
double  c1, c2, c3, c4, c5, u1, u2, w;

c1=alpha1-1.0;
c2=(alpha1-(1.0/(6.0*alpha1)))/c1;
c3=2.0/c1;
c4=c3+2.0;
c5=1.0/sqrt(alpha1);
do{
	do{
		u1=WichHill();
		u2=WichHill();
		if (alpha1>2.5){
			u1=u2+c5*(1.0-(1.86*u1));
		}
	}
	while (u1>=1.0 || u1<=0);

	w=(c2*u2)/u1;
	if ((c3*u1+w+(1.0/w))<=c4){
		acpt=1;
	}
	else if ((c3*log(u1)-log(w)+w)>=1.0){ /* repeat outer "do" loop */
		acpt=0;
	}
	else{
		acpt=1;
	}
	}
while (acpt!=1);

return c1*w;
}

double rbetaf()
{
/* returns a Beta random variable (from the prior distribution) to
   initialize f, with parameters determined by the minimum of the allele
   frequencies "pmin", on the interval [-pmin/(1-pmin),1]
*/

double rgm1, rgm2, bpa, bpb, blb;

blb= -pmin/(1.0-pmin);

if (btail==0){
	bpb=1.0;
}
else if (btail==1){
	bpb=61.37+1482.08*pmin-1643.18*pow(pmin,2.0)+41521.78*pow(pmin,3.0)-298223.0*pow(pmin,4.0)+723583.9*pow(pmin,5.0)-579512.0*pow(pmin,6.0);
}
else if (btail==2){
	bpb=23.61+241.3*pmin-1048.65*pow(pmin,2.0)+8944.7*pow(pmin,3.0)-20820.2*pow(pmin,4.0)+14112.95*pow(pmin,5.0);
}
else if (btail==4){
	bpb=11.15+29.75*pmin+773.84*pow(pmin,2.0)-4409.76*pow(pmin,3.0)+9156.72*pow(pmin,4.0)-6515.75*pow(pmin,5.0);
}
bpa= 1.0-blb*(bpb-1.0);

if (bpa>1.0){
	rgm1=rstgam(bpa);
}
else if (bpa<1.0){    /* Ex. 3.22 in Ripley(1988) for alpha<1 */
	rgm1=rstgam(bpa+1.0)*pow(WichHill(),1.0/bpa);
}
else {    /* Exponential */
	rgm1=(-log(WichHill()));
}
if (bpb>1.0){
	rgm2=rstgam(bpb);
}
else if (bpb<1.0){    /* Ex. 3.22 in Ripley(1988) for alpha<1 */
	rgm2=rstgam(bpb+1.0)*pow(WichHill(),1.0/bpb);
}
else {    /* Exponential */
	rgm2=(-log(WichHill()));
}

return (rgm1/(rgm1+rgm2))*(1.0-blb)+blb;
}

double llkhd(double fs,pvec ps)
{
/* evaluates the log of the joint likelihood (excluding multinomial
   constant)
*/

double  lld=0.0;

for (i=0;i<nall;++i){
	for (j=i+1;j<nall;++j){
		genf[i+(j*(j+1))/2]=2.0*(1-fs)*ps[i]*ps[j];
	}
}
for (i=0;i<nall;++i){
	genf[(i*(i+3))/2]=ps[i]*(fs+(1-fs)*ps[i]);
}

/* log likelihood */

for (i=0;i<ngen;++i){
	if (n[i]!=0){
		lld+= n[i]*log(genf[i]);
	}
}

return lld;
}

main()
{
srand((unsigned) time(NULL));
out2=fopen("hwmet.out","w");

/* initialize seeds for random number generator */

ix=(long)(rand()%30000)+1;
iy=(long)(rand()%30000)+1;
iz=(long)(rand()%30000)+1;

inpdat();

/* initialize "p" and "f" values */

nnn=0;
for (i=0,gsum=0;i<nall;++i){
	p[i]= (-log(WichHill()));
	gsum+= p[i];
}
for (i=0;i<nall;++i){
	p[i]=p[i]/gsum;
}

pmin=1.0/nall;
for (i=0;i<nall;++i){
	if (p[i]<pmin){
		pmin=p[i];
	}
}
f=rbetaf();

/* initialize density value "z" */

z=lbetad(f,pmin)+llkhd(f,p);

/* begin Metropolis-Hastings algorithm */

for (l1=1;l1<=numit+discard;++l1){

/* choose new candidate values "newp" and "newf" */

	p1=(int)(WichHill()*nall);
	do {
		p2=(int)(WichHill()*nall);
	}
	while (p2==p1);

	b1=(p[p1]+epsln<p[p1]+p[p2])?p[p1]+epsln:p[p1]+p[p2];
	if (b1>1.0){  /* can only occur through rounding error */
		b1=1.0;
	}
	a1=(p[p1]-epsln>0)?p[p1]-epsln:0;
	for (i=0;i<nall;++i){
		newp[i]=p[i];
	}
	newp[p1]= (b1-a1)*WichHill()+a1;
	newp[p2]-= newp[p1]-p[p1];

	newpmin=1.0/nall;
	for (i=0;i<nall;++i){
		if (newp[i]<newpmin){
			newpmin=newp[i];
		}
	}

	b3=(1.0<f+epsln2)?1.0:f+epsln2;
	a3=(-newpmin/(1.0-newpmin)>f-epsln2)?-newpmin/(1.0-newpmin):f-epsln2;
	newf=a3+WichHill()*(b3-a3);

/* evaluate density "newz" at candidate values */

	newz=lbetad(newf,newpmin)+llkhd(newf,newp);

/* decide whether or not to accept candidate values */

	b2=(newp[p1]+epsln<newp[p1]+newp[p2])?newp[p1]+epsln:newp[p1]+newp[p2];
	if (b1>1.0){  /* can only occur through rounding error */
		b1=1.0;
	}
	a2=(newp[p1]-epsln>0)?newp[p1]-epsln:0;

	b4=(1.0<newf+epsln2)?1.0:newf+epsln2;
	a4=(-pmin/(1.0-pmin)>newf-epsln2)?-pmin/(1.0-pmin):newf-epsln2;

	r=newz-z+log(b1-a1)-log(b2-a2)+log(b3-a3)-log(b4-a4);

	if (log(WichHill())<r){ /* accept candidate values */
		pmin=newpmin;
		z=newz;
		for (i=0;i<nall;++i){
			p[i]=newp[i];
		}
		f=newf;
		jmp++;
	}

/* output */

    if (((l1/kp)*kp==l1) && (l1>discard)){
		fprintf(out2,"%10.6f\n",f);
		fflush(out2);
        fst[nnn]=f;
        nnn++;
    }            

	if ((l1/acrt)*acrt==l1){     /* output acceptance rate */
        printf("%ld iterations of %ld completed ...\n",l1,numit);
		fprintf(out1,"%12.6f\n",((double)(jmp))/l1);
		fflush(out1);
	}
}

printf("PROGRAM COMPLETE - f VALUES STORED IN hwmet.out\n");
printf("                 - LOG AND SUMMARY STATISTICS STORED IN hwmet.stats\n");

mnf=0.0;
for (l1=0;l1<nout;++l1){
	mnf+= fst[l1];
}
mnf/= nout;

for (l1=0;l1<nout-1;++l1){
   for (jmp=l1+1;jmp<nout;++jmp){
      if (fst[l1]>fst[jmp]){
         tmpd=fst[l1];
         fst[l1]=fst[jmp];
         fst[jmp]=tmpd;
      }
   }
}

/* CALCULATE PERCENTILES BY SETTING i'th ORDER STATISTIC EQUAL TO
   (i-0.5)/n AND USING LINEAR INTERPOLATION
*/

fprintf(out1,"\n\n****** POSTERIOR DISTRIBUTION ******\n");

fprintf(out1,"\n** Summary stats, based on %ld simulated values, for posterior of f **\n\n",nout);

fprintf(out1,"---------------------------------------------------------------------------------------------------\n");
fprintf(out1,"|         | Mean     |2.5 %c-ile |  5 %c-ile | 25 %c-ile |  median  | 75 %c-ile | 95 %c-ile |97.5 %c-ile|\n",37,37,37,37,37,37);
fprintf(out1,"|-------------------------------------------------------------------------------------------------|\n");

if (mnf>0){
	fprintf(out1,"|Posterior| %9.7f|",mnf);
}
else {
	fprintf(out1,"|Posterior|%9.7f|",mnf);
}

/* 2.5th */

u=0.025*nout+0.5;
k=(int)(u);
if (k==u){
	if (fst[k]>0){
	   fprintf(out1," %9.7f|",fst[k]);
	}
	else {
	   fprintf(out1,"%9.7f|",fst[k]);
	}
}
else {
	if ((fst[k]+fmod(u,1.0)*(fst[k+1]-fst[k]))>0){
	   fprintf(out1," %9.7f|",fst[k]+fmod(u,1.0)*(fst[k+1]-fst[k]));
	}
	else {
	   fprintf(out1,"%9.7f|",fst[k]+fmod(u,1.0)*(fst[k+1]-fst[k]));
	}
}         

/* 5th */

u=0.05*nout+0.5;
k=(int)(u);
if (k==u){
	if (fst[k]>0){
	   fprintf(out1," %9.7f|",fst[k]);
	}
	else {
	   fprintf(out1,"%9.7f|",fst[k]);
	}
}
else {
	if ((fst[k]+fmod(u,1.0)*(fst[k+1]-fst[k]))>0){
	   fprintf(out1," %9.7f|",fst[k]+fmod(u,1.0)*(fst[k+1]-fst[k]));
	}
	else {
	   fprintf(out1,"%9.7f|",fst[k]+fmod(u,1.0)*(fst[k+1]-fst[k]));
	}
}          

/* 25th */

u=0.25*nout+0.5;
k=(int)(u);
if (k==u){
	if (fst[k]>0){
	   fprintf(out1," %9.7f|",fst[k]);
	}
	else {
	   fprintf(out1,"%9.7f|",fst[k]);
	}
}
else {
	if ((fst[k]+fmod(u,1.0)*(fst[k+1]-fst[k]))>0){
	   fprintf(out1," %9.7f|",fst[k]+fmod(u,1.0)*(fst[k+1]-fst[k]));
	}
	else {
	   fprintf(out1,"%9.7f|",fst[k]+fmod(u,1.0)*(fst[k+1]-fst[k]));
	}
}         

/* 50th */

u=0.5*nout+0.5;
k=(int)(u);
if (k==u){
	if (fst[k]>0){
	   fprintf(out1," %9.7f|",fst[k]);
	}
	else {
	   fprintf(out1,"%9.7f|",fst[k]);
	}
}
else {
	if ((fst[k]+fmod(u,1.0)*(fst[k+1]-fst[k]))>0){
	   fprintf(out1," %9.7f|",fst[k]+fmod(u,1.0)*(fst[k+1]-fst[k]));
	}
	else {
	   fprintf(out1,"%9.7f|",fst[k]+fmod(u,1.0)*(fst[k+1]-fst[k]));
	}
}           

/* 75th */

u=0.75*nout+0.5;
k=(int)(u);
if (k==u){
	if (fst[k]>0){
	   fprintf(out1," %9.7f|",fst[k]);
	}
	else {
	   fprintf(out1,"%9.7f|",fst[k]);
	}
}
else {
	if ((fst[k]+fmod(u,1.0)*(fst[k+1]-fst[k]))>0){
	   fprintf(out1," %9.7f|",fst[k]+fmod(u,1.0)*(fst[k+1]-fst[k]));
	}
	else {
	   fprintf(out1,"%9.7f|",fst[k]+fmod(u,1.0)*(fst[k+1]-fst[k]));
	}
}        

/* 95th */

u=0.95*nout+0.5;
k=(int)(u);
if (k==u){
	if (fst[k]>0){
	   fprintf(out1," %9.7f|",fst[k]);
	}
	else {
	   fprintf(out1,"%9.7f|",fst[k]);
	}
}
else {
	if ((fst[k]+fmod(u,1.0)*(fst[k+1]-fst[k]))>0){
	   fprintf(out1," %9.7f|",fst[k]+fmod(u,1.0)*(fst[k+1]-fst[k]));
	}
	else {
	   fprintf(out1,"%9.7f|",fst[k]+fmod(u,1.0)*(fst[k+1]-fst[k]));
	}
}          

/* 97.5th */

u=0.975*nout+0.5;
k=(int)(u);
if (k==u){
	if (fst[k]>0){
	   fprintf(out1," %9.7f|",fst[k]);
	}
	else {
	   fprintf(out1,"%9.7f|",fst[k]);
	}
}
else {
	if ((fst[k]+fmod(u,1.0)*(fst[k+1]-fst[k]))>0){
	   fprintf(out1," %9.7f|",fst[k]+fmod(u,1.0)*(fst[k+1]-fst[k]));
	}
	else {
	   fprintf(out1,"%9.7f|",fst[k]+fmod(u,1.0)*(fst[k+1]-fst[k]));
	}
}          

fprintf(out1,"---------------------------------------------------------------------------------------------------\n\n");

fclose(out1);
fclose(out2);

return 0;
}