c
c	FORTRAN PROGRAM TO ESTIMATE THE AREA UNDER THE ROC CURVE		
c		FOR CLUSTERED DATA
c
c
c
	implicit real (a-z)
	integer n_clust,cl_size,i,j,k,l,truth(500,10),id(500)
	integer count,sim,abn,nor,kount,compar
	dimension test(500,10),vx(500,10),vy(500,10)
	logical d_nor,d_abn
c
c
	open(unit=1, file='data.dat', status='old')
c	open(unit=2, file='cluster.out', status='new')
c
	do 1000, sim=1,1
c
	print *,'Enter the total number of patients (i.e. # clusters)'
	read *,n_clust
	print *,'Enter the maximum cluster size'
	read *,cl_size
	abn=0
	nor=0
c
 20	format(i3,1x,i2)
 21	format(i3,1x,i1,1x,i3)
	print *,'n_clust',n_clust,'cl_size',cl_size
c
	do 50, i=1, n_clust
	do 45, k=1,cl_size
	read(1,21)id(i),truth(i,k),test(i,k)
c	print *,id(i),truth(i,k),test(i,k)
 45	continue
 50	continue
c
c
c	Compute the X-structural components
c
	do 60, i=1,n_clust
	do 59, k=1,cl_size
	vx(i,k)=99
 59	continue
 60	continue
c
c
c
	do 100, i=1,n_clust
	do 95, k=1,cl_size
	sum=0.0
	count=0
c
	if(truth(i,k) .eq. 1)goto 95
	abn=abn+1
c 
	do 90, j=1,n_clust
	do 85, l=1,cl_size
c
	if(truth(j,l) .eq. 2)goto 85
c
	count=count+1
	if(test(i,k) .gt. test(j,l))sum=sum+1
	if(test(i,k) .eq. test(j,l))sum=sum+0.5
c
c	
 85	continue
 90	continue
c
c
	vx(i,k)=sum/real(count)
c
c
 95	continue
 100	continue
c
c
c	Compute the Y-structural components
c
c
	do 120, i=1,n_clust
	do 119, k=1,cl_size
	vy(i,k)=99
 119	continue
 120	continue
c
c
	do 200, i=1,n_clust
	do 195, k=1,cl_size
	sum=0.0
	count=0
c
	if(truth(i,k) .eq. 2)goto 195
	nor=nor+1
c 
	do 190, j=1,n_clust
	do 185, l=1,cl_size
c
	if(truth(j,l) .eq. 1)goto 185
c
	count=count+1
	if(test(i,k) .lt. test(j,l))sum=sum+1
	if(test(i,k) .eq. test(j,l))sum=sum+0.5
c
c	
 185	continue
 190	continue
c
c
	vy(i,k)=sum/real(count)
c
c
 195	continue
 200	continue
c
c
c	Compute the area
c
c
	count=0
	sum=0.0
c
	do 250, i=1,n_clust
	do 249, k=1,cl_size
	if(vx(i,k) .gt. 1.0)goto 249
	sum=sum+vx(i,k)
	count=count+1
 249	continue
 250	continue
c
	area=sum/real(count)
	print *,'Estimated Area is ',area
c
c
c
c	Compute the variance term S_10
c
c
	sqsum=0.0
	kount=0
c
c
	do 275, i=1,n_clust
	count=0
	sum=0.0
	do 265, k=1,cl_size
	if(vx(i,k) .le. 1.0)then
		count=count+1
		sum=sum+vx(i,k)
	endif
 265	continue
	if(count .eq. 0)goto 275
	kount=kount+1
	sqsum=sqsum+(sum-real(count)*area)**2
 275	continue
c
	s_10=sqsum*real(kount)/real((abn)*(kount-1))
	print *,'s_10',s_10
	print *,'I10',kount
c
c
c
c	Compute the variance term S_01
c
c
	sqsum=0.0
	kount=0
c
c
	do 295, i=1,n_clust
	count=0
	sum=0.0
c
	do 285, k=1,cl_size
	if(vy(i,k) .le. 1.0)then
		count=count+1
		sum=sum+vy(i,k)
	endif
 285	continue
	if(count .eq. 0)goto 295
	kount=kount+1
	sqsum=sqsum+(sum-real(count)*area)**2
 295	continue
c
	s_01=sqsum*real(kount)/real((nor)*(kount-1))
	print *,'s_01',s_01
	print *,'I01',kount
c
c
c
c	Compute the covariance term, s_11
c
c
	sqsum=0.0
	kount=0
	compar=0
c
c
	do 320, i=1,n_clust
c
	count1=0
	sum1=0.0
	count2=0
	sum2=0.0
	d_abn=.false.
	d_nor=.false.
c
	do 310, k=1,cl_size
	if(vx(i,k) .le. 1.0)then
		d_abn=.true.
		count1=count1+1
		sum1=sum1+vx(i,k)
	endif
	if(vy(i,k) .le. 1.0)then
		d_nor=.true.
		count2=count2+1
		sum2=sum2+vy(i,k)
	endif
 310	continue
c
	if(d_abn .eq. .false. .or. d_nor .eq. .false.)goto 320
	kount=kount+1
	compar=compar+(count1*count2)
	sqsum=sqsum+(sum1-real(count1)*area)*(sum2-
     + real(count2)*area)
 320	continue
c
	if(kount .eq. 0)then
		s_11=0.00
		compar=99
		goto 325
	endif
	s_11=sqsum*real(n_clust)/real(n_clust-1)
	print *,'s_11',s_11
c	print *,'sqsum',sqsum
c
c
 325	continue
c
c
c	Compute the Variance of the area
c
c
	var=(1.0/real(abn))*s_10+(1.0/real(nor))*s_01+
     + (2.0/real(abn*nor))*s_11
	print *,'Estimated variance of area is',var
	print *,'N',abn,'M',nor
c
c
c
c	write(2,500)sim,area,vara,nor,abn,compar,kount,
c    + s_01,s_10,s_11a,varb,varc
c
 500	format(i4,1x,f6.4,1x,f8.6,1x,4(i3,1x),5(f8.6,1x))
c
c
c
 1000	continue
c
	stop
	end