#include "k.h"
/*  C implementation of Gibbs sampler  
	input - kdata.txt
	      - tables of A1 alleles (1 row of 6 values per table)
	        cases followed by controls
	        r0 r1 r2 s0 s1 s2 
    output - rslt.txt
           - a vector of prob obs is an outlier (length n) */
void main() {
double	jnk[6];				/* Work for fetching data	*/
double	xsumsq=0.,lambda,p,num,r,fun1,fun2,cc,dd,rsig,
		*x,*A,*out,*deltot,zeta,kappa,tau2,epsilon;
FILE *pinfile;
FILE *rslttxt; 
int i,j,R,N,
	lcnt=0,				/* input file line counter	*/
	n=0,				/* data counter				*/
	nu=0,*delta,
	df=0,				/* rchisq degrees of freedom*/
	itmp=0,	
	o[3],e[3],          /* work for allele table    */			
	ngib,burn;	 
long idum=2348;			/* random seed 				*/

/* hardcoded constants									*/
zeta=1000.;kappa=4.;tau2=1.;epsilon=.01;ngib=500;burn=50;	 

rslttxt=fopen("rslt.txt","w+"); 
 
if ((pinfile=fopen("kdata.txt","rt"))==NULL) {
    printf("I can't open file.");
    exit(1);
} 
/* count no. of tables in file */

n=getnum(pinfile,jnk);
while (!feof(pinfile)) n+=getnum(pinfile,jnk);
rewind(pinfile);
		
/* allocate memory 	*/

if ((x=(double *)malloc(n*sizeof(double)))==NULL) {
  printf("I can't allocate memory.");
  exit(1);
}
if ((A=(double *)malloc(n*sizeof(double)))==NULL) {
  printf("I can't allocate memory.\n");
  exit(1);
}
if ((deltot=(double *)malloc(n*sizeof(double)))==NULL) {
  printf("I can't allocate memory.\n");
  exit(1);
}
if ((delta=(int *)malloc(n*sizeof(int)))==NULL) {
  printf("I can't allocate memory.\n");
  exit(1);
}
if ((out=(double *)malloc((burn+ngib)*sizeof(double)))==NULL) {
  printf("I can't allocate memory.\n");
  exit(1);
}
for (i=0;i<n;i++) {
	x[i]=0.; A[i]=0.; deltot[i]=0.; delta[i]=0;;
}
/* read data  */

rewind(pinfile);
n=0;
while (!feof(pinfile)) {
    itmp=getnum(pinfile,jnk);
	if (itmp!=6) printf("Bad data (n=%d) at line %d\n",itmp,n+1);
	else { 
	    R=0;N=0;
		for (i=0;i<3;i++) {
			o[i]=(int)jnk[i]; 
			e[i]=(int)jnk[i+3];
			R+=o[i];
			N+=o[i]+e[i];
		}
		x[n++]=chi(o,e,R,N);
	}
}
fclose(pinfile);
printf("\n%d lines read.\n\n",n);
rsig=median(x,&n)/.675;
kappa=kappa*rsig;
tau2=tau2*rsig*rsig;
for (i=0;i<n;i++) xsumsq+=x[i]*x[i];
lambda=xsumsq/n;
dd=1./(1./tau2+1./lambda);
for (i=0;i<n;i++) { 
	cc=dd*(x[i]/lambda+kappa/tau2);
	A[i]=dd*rnorm(&idum)+cc;
	fun1=(x[i]-A[i])/sqrt(lambda);
	fun2=x[i]/sqrt(lambda);
	p=dnorm(fun1)*epsilon/(dnorm(fun1)*epsilon+dnorm(fun2)*(1.-epsilon));
	delta[i]=(p>runif(&idum));
}
/* gibbs update */
for (i=0;i<burn+ngib;i++) {
	/* calc num <- sum((x-delta*A)^2+nu*zeta  */
	/* and df <- sum(1-delta)+nu for rchisq		*/				 
	num=(double)nu*zeta;
	df=nu;
	for (j=0;j<n;j++) {
		num+=pow(x[j]-delta[j]*A[j],2); 
	    df+=(1-delta[j]);
	} 
	lambda=num/rchisq((long)df,&idum);
	dd=1/(1/tau2+1/lambda); /* dd is used in loop below */
	r=0;
	for (j=0;j<n;j++) {
		fun1=(x[j]-A[j])/sqrt(lambda);
		fun2=x[j]/sqrt(lambda);
		p=dnorm(fun1)*epsilon/(dnorm(fun1)*epsilon+dnorm(fun2)*(1.-epsilon));
		delta[j]=(p>runif(&idum));
		if (i>=burn) deltot[j]+=p;
		if(i>=burn && j==0)	printf("%14.9f \n",p);
		/* dd=1/(1/tau2+1/lambda) doesn't change so it's before the loop */
		cc=dd*(x[j]/lambda+kappa/tau2);
		A[j]=delta[j]*(sqrt(dd)*rnorm(&idum)+cc)+
			(1-delta[j])*(sqrt(tau2)*rnorm(&idum)+kappa);
		r+=delta[j];
	} 
}
for (j=0;j<n;j++) {
 if(j==0)	printf("%d %14.9f \n",ngib,deltot[j]);
	deltot[j]=deltot[j]/ngib;
	fprintf(rslttxt,"%14.9f\n",deltot[j]);
}
fclose(rslttxt);  	
return;
}		   
/* Put numbers fm inputfile line into double array; return number read.*/
int getnum(FILE *f, double *x) {
    static char buffer[MAX_INPUT+1];
    if (fgets(buffer,MAX_INPUT,f)==NULL) return(0);
    return(sscanf(buffer,"%lf %lf %lf %lf %lf %lf",x,x+1,x+2,x+3,x+4,x+5));
}
double dnorm(double x) {
  return(exp(-.5*x*x)/sqrt(2.*PI));
}