#include "k.h"
#define r 1 	/* length of level vector */
/* 	Bayes test for outlier - declare as outlier if prob>.5
*	Compare type I and II rates with 
*	Bonferroni correction. Average over many reps. */
void main() {
double *x,*xt,*pp,level[r]={.50},gamma,F1,F2,
	*dwork,								/* work space arrays for prog 	*/
	out1[r+1],							/* for output, out column1 		*/
	out2[r+1];							/* for output, out column2 		*/
int i,j,iter,*iwork,					/* loop variables		   		*/
	n,R,S,oo,reps;							
long idum;

/* hardcoded constants													*/
R=100;S=100;gamma=2.5;F1=.001;F2=.0001;oo=10;reps=1;n=400;  			       

/* read in R S gamma F1 F2 oo reps */
/* scanf("%d %d %lf %lf %lf %d %d",&R,&S,&gamma,&F1,&F2,&oo,&reps);		*/	

/* allocate work space (A in prog) */
if ((dwork=(double *)malloc(n*sizeof(double)))==NULL) {
  printf("I can't allocate memory. (perform.c)\n");
  exit(1);
}
/* allocate work space (delta in prog) */
if ((iwork=(int *)malloc(n*sizeof(int)))==NULL) {
  printf("I can't allocate memory. (perform.c)\n");
  exit(1);
}
if ((x=(double *)malloc(n*sizeof(double)))==NULL) {
  printf("I can't allocate memory. (perform.c)\n");
  exit(1);
}
if ((xt=(double *)malloc(n*sizeof(double)))==NULL) {
  printf("I can't allocate memory. (perform.c)\n");
  exit(1);
}
if ((pp=(double *)malloc(n*sizeof(double)))==NULL) {
  printf("I can't allocate memory. (perform.c)\n");
  exit(1);
}			
for (i=0;i<r+1;i++) out1[i]=out2[i]=0.;	/* init out1,out2 		*/
for (iter=0;iter<reps;iter++) {
/*        printf("%d\n",iter);	*/
	datagen(n,R,S,gamma,F1,F2,oo,x,xt,&idum);
	for (j=0;j<n;j++) pp[j]=0.;		/* init pp				*/
	prog(x,n,pp,&idum,dwork,iwork); 
	for (i=0;i<r;i++) { 
		for (j=0;j<(n-oo);j++) { 
			/* indicator = 0 */
			out2[i]+=(pp[j]>level[i]); /* add truth to out[i,2] */
		}
		for (j=(n-oo);j<n;j++) { 
			/* indicator = 1 */
			out1[i]+=(pp[j]>level[i]); /* add truth to out[i,1] */
		}	
	}
 	for (j=0;j<(n-oo);j++) { 
		/* indicator = 0 */
		out2[r]+=(x[j]>3.66); /* add BC to out[r+1,2] */
	}
	for (j=(n-oo);j<n;j++) { 
		/* indicator = 1 */
		out1[r]+=(x[j]>3.66); /* add BC to out[r+1,1] */
	}
}
for (i=0;i<r;i++) printf("%5.2f %16.8f %16.8f\n",level[i],out1[i]/(double)reps,out2[i]/(double)reps);
printf("%5.2f %16.8f %16.8f\n",0.,out1[r]/(double)reps,out2[r]/(double)reps);
/* tmp  _getch(); */
return;
}