program PCA_driver

!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

! A program intended to drive the general routine for performing a PCA

!analysis on a dataset.

!			version 2/5/02          DRH

!________________________________________________________________________

USE DFIMSLSS

USE IMSLf90 

implicit none



integer, parameter :: rows=59  !max size of data arrays

integer, parameter :: p=7  !look at your data file and check the # of predictors



integer, external :: choose



integer :: i,j,k,g,n,q,o,ioerror,trial(p),place,col

integer :: incl_vars(p)



real ::  X(rows,p),y(rows,1), yhat(rows-1,1) !

real :: min,inter,yavg,Xavg,sumX,covar,Sx(p,p),id(1,p)

real :: L(p,1),U(p,p),Uinv(p,p),D(p,p),sumL,sumvarX,Ut(p,p)

real :: X_XbarT(rows,p),bigY(rows,p),Sy(p,p),Ovec(p),sumY



character (len=12) :: model(p),template(p)

data template /'JanTemp','JulyTemp','Rain','PopDensity',&

'Education','NonWhite','logHC'/

data Ovec /'0.0','0.0','0.0','0.0','0.0','0.0','0.0'/





open (4,file='PCA.txt') !output file

open (20,file='Mortality_dataset.txt',status='old') !datafile





n=0

ioerror=0

do i=1,rows

	n=n+1	

	read (20,*,iostat=ioerror) y(i,1),(X(i,j),j=1,p)

	if(ioerror==-1) then

		n=n-1

		write(6,'(a,i3)') 'end of file reached, n= ',n

		exit

	else if (ioerror/=0) then

		write(6,'(a,i3)') 'read error!, line ',i

	endif

	!write(4,'(<p>(F9.2,1x))') (X(i,j),j=1,p)

enddo

!**********************************************************************

write(4,'(a)') 'covariance matrix of x = '

do i=1,p

	do j=1,p

		call covarience(n,X(:,i),X(:,j),covar)

		Sx(i,j) = covar

	end do

	write(4,'(<p>(f10.2,1x))') Sx(i,:)

end do



!We need to do an eigan decomposition of Sx now.

!this routine is for the eigan decomposition of a symmetric matrix

call EVCSF(p,Sx,p,L,U,p)



do i=1,p

	call average(n,X(:,i),Xavg)

	do j=1,n

		!transpose and !subtract mean in one step

		!X-XbarT(j,i) = X(i,j)-Xavg

		!i don't think we are supposed to transpose

		X_XbarT(j,i) = X(j,i)-Xavg

	end do

end do

bigY = X_XbarT .x. U 





!***********************************************************************

!Checks to see if this has worked

sumL = sum(L)

do i=1,p

	call varience(n,X(:,i),Xavg)

	sumvarX = sumvarX+Xavg

end do



write(4,'(/,a,e14.9)') 'the sum of the eigan values, ',sumL

write(4,'(a,e14.9)') "should be equal to the sum of the variences of the X's, ",sumvarX



Uinv = .i. U

Ut= .t. U

write(4,'(/,a)') "the product UU' should = I, UU' =  "

D = U .x. Ut

do i=1,p

	write(4,'(<p>(f8.5,1x))') (D(i,j),j=1,p)

end do



sumY = sum(bigY)

write(4,'(/,a)') 'covariance matrix of PCA scores, '

do i=1,p

	do j=1,p

		call covarience(n,bigY(:,i),bigY(:,j),covar)

		Sy(i,j) = covar

	end do

	write(4,'(<p>(f10.2,1x))') Sy(i,:)

end do

write(4,'(/,a)') 'should equal the diagonal matrix of eigan values'



D = 0.0

do i=1,p

	D(i,i) = L(i,1)

	write(4,'(<p>(f10.2,1x))') (D(i,j),j=1,p)

end do



	stop 'normal completion'

end program PCA_driver