Real, Parameter :: Miss=-99. Real, Allocatable :: Xbar(:), Std(:), Rx(:), D(:), V(:) Real, Allocatable :: X(:,:), Y(:,:), YT(:,:) Real, Allocatable :: Fx(:,:), Ax(:,:) ! X = Fx*A'x Real, Allocatable :: EOF(:,:), PCA(:,:) !For testing: Real, Allocatable :: FF(:,:), AA(:,:) Real, Allocatable :: X1(:,:), EVtor2(:,:) Real, Allocatable :: Modul_XE(:), YdX(:,:) Real, Allocatable :: EVtor(:,:), Lamda(:,:) Real, Allocatable :: EValues(:) Integer :: Nt, Nx Integer :: Mcount, Ncount, switch, Standard, Nth0, Nthre Real :: Threshold Character (Len=250) :: Fname Real :: Factor Namelist /Model_Info/ Fname, switch, standard, Threshold, Nth0 Integer :: i,j,k Open (1,file="Namelist_PCA_New.txt",Status="old") Read (1, Model_Info) Nt=10; Nx=15 Factor = Sqrt(Real(Nt)) Allocate (X(Nt,Nx), Y(Nt,Nx), YT(Nx,Nt)) ! SL ban dau va SL quy tam hoac chuan hoa Open (1,file=fname) read (1,*) Do i=1, Nt read (1,*) (X(i,j),j=1,Nx) Enddo if (switch==1) then Mcount=Nt; Ncount=Nx else if (switch==0) then Mcount=Nx; Ncount=Nt endif !For testing: Allocate (FF(Mcount,Mcount), AA(Mcount,Mcount)) Allocate (X1(Nt,Nx), EVtor2(Mcount,Mcount)) Allocate (YdX(Nt,Nx)) !Allocate (X(Nt,Nx), Y(Nt,Nx), YT(Nx,Nt)) ! SL ban dau va SL quy tam hoac chuan hoa Allocate (Xbar(Mcount), Std(Mcount), Rx(Mcount*Mcount)) ! TB, DLC, ma tran Covar hoac Corre Allocate (EValues(Mcount), V(Mcount*Mcount)) ! Tri rieng va Vecto rieng (1D) Allocate (Fx(Ncount,Mcount), Ax(Mcount,Mcount)) ! EOF va PCA Allocate (D(Mcount), EVtor(Mcount,Mcount)) ! Working arr va matran Vecto rieng Allocate (Lamda(Mcount,Mcount)) Call Covar(Nt,Nx, Miss) Call EIGEN(Mcount,MV) Call Calc_Thresholds (Nth0, Threshold, CON) ! !tmp = 0. !do I=1,Mcount ! L=I+(I*I-I)/2 ! EValues(I) = Rx(L) ! tmp = tmp+EValues(I) ! Tong cac tri rieng !enddo !print*,"Tong tri rieng=",tmp !tmp1 = 0; Nthre = 0 ! De so voi Tong cac tri rieng !Do i=1,Mcount ! tmp1 = tmp1+EValues(I); Nthre=Nthre+1 ! if (tmp1/tmp > Threshold .OR. Nthre >= Nth0) exit ! Vuot nguong nhung khong qua Nth0 nhan to !Enddo !CON = EValues(Nthre) ! Xac dinh lai nguong gia tri cua Tri rieng duoc giu lai !print*,"Evalues=",EValues(1:Nth0) !print*,"CON=",CON !WRITE(*,*)' TRI RIENG 1,2,3,4,5=',Rx(1),Rx(3),Rx(6),Rx(10),Rx(15) !write(*,'(A,F10.3)') "Gia tri nguong CON: ", CON !write(*,'(A,F10.3)') "So phan tram luong thong tin duoc giu lai: ",tmp1/tmp*100. Call TRACE(Mcount,CON,K) Call EigenValues_and_EigenVectors Allocate (EOF(Nt,K), PCA(Nx,K)) !print*,"K from Trace=",K !DO I=1,K ! L=I+(I*I-I)/2 ! EValues(I)=Rx(L) !ENDDO !L=0 !DO J=1,K ! DO I=1,Mcount ! L=L+1 ! D(I)=V(L) ! EVtor(i,j)=V(L) ! Ma tran Vector rieng, moi cot 1 vector ! ENDDO !print "(A,I3,20F8.3)","EVT=",j, (EVtor(i,j),i=1,Mcount) !ENDDO CALL LOAD(Mcount,K) Call Calc_AX !DO I=1,Mcount ! DO J=1,K ! L=Mcount*(J-1)+I ! D(J)=V(L) ! Ax(I,J)=V(L) ! Ma tran PCs, moi cot 1 PC !print*,i,j,L ! ENDDO !Ax(I,1:K)=D(1:K) !ENDDO !If (K > 1) Then ! NV=NC+1 Call Calc_EOF_PCA_Standard !EndIF ! For reference: ! Calculating the EOFs (F): X(N,M) x A(M,K) = F(N,K) ! Xpre(j) = (Xorg(iyr,j)-Mean(imon,j))/Std(imon,j)/tmp; ! Fpre=Matmul(Xpre,A) ! Do j=1, NEigenMax; Fpre(j)=Fpre(j)/Eig(imon,j); Enddo ! EOFs(iyr,1:NEigenMax,imon)=Fpre(1:NEigenMax) !Y =MatMul(Fx(1:Ncount,1:K),Transpose(Ax(1:Mcount,1:K))) FF=MatMul(Transpose(EOF),EOF) AA=MatMul(Transpose(PCA),PCA) if (switch==0) then X1=MatMul(EOF,Transpose(PCA))*Factor else if (switch==1) then X1=MatMul(EOF,Transpose(PCA))*Factor Do j=1,Nx Do i=1,Nt YdX(i,j)=Y(i,j)/X1(i,j) Enddo Enddo endif if (standard==1) then Do j=1, Mcount ! X1(:,j)=X1(:,j)*Factor Enddo endif EVtor2=MatMul(Transpose(EVtor),EVtor) Call Write_test Deallocate (D) Contains !--------------------------------------------------- Subroutine Calc_AX Integer :: i,j DO I=1,Mcount DO J=1,K L=Mcount*(J-1)+I D(J)=V(L) Ax(I,J)=V(L) ! Ma tran PCs, moi cot 1 PC !Ax(i,j) = EVtor(i,j) ! Thu voi PCs chuan hoa ENDDO ENDDO End Subroutine Calc_AX !--------------------------------------------------- Subroutine EigenValues_and_EigenVectors Integer :: i, L DO I=1,K L=I+(I*I-I)/2 EValues(I)=Rx(L) ENDDO L=0 DO J=1,K DO I=1,Mcount L=L+1 D(I)=V(L) EVtor(i,j)=V(L) ! Ma tran Vector rieng, moi cot 1 vector ENDDO ENDDO End Subroutine EigenValues_and_EigenVectors !--------------------------------------------------- Subroutine Calc_Thresholds (Nth0, Threshold, CON) ! Input: Number of Eigenvalues (Nth0) and Threshold value (Threshold) of Eigenvalues ! Output: CON satisfying N-eigenvalue>=Nth0 OR accumulated Eigenvalues>Threshold Real :: tmp, tmp1, Threshold, CON Integer :: i, L, Nthre, Nth0 tmp = 0. do I=1,Mcount L=I+(I*I-I)/2 EValues(I) = Rx(L) tmp = tmp+EValues(I) ! Tong cac tri rieng enddo print*,"Tong tri rieng=",tmp tmp1 = 0; Nthre = 0 ! De so voi Tong cac tri rieng Do i=1,Mcount tmp1 = tmp1+EValues(I); Nthre=Nthre+1 if (tmp1/tmp > Threshold .OR. Nthre >= Nth0) exit ! Vuot nguong nhung khong qua Nth0 nhan to Enddo CON = EValues(Nthre) ! Xac dinh lai nguong gia tri cua Tri rieng duoc giu lai print*,"Evalues=",EValues(1:Nth0) print*,"CON=",CON WRITE(*,*)' TRI RIENG 1,2,3,4,5=',Rx(1),Rx(3),Rx(6),Rx(10),Rx(15) write(*,'(A,F10.3)') "Gia tri nguong CON: ", CON write(*,'(A,F10.3)') "So phan tram luong thong tin duoc giu lai: ",tmp1/tmp*100. End Subroutine Calc_Thresholds !--------------------------------------------------- Subroutine Calc_EOF_PCA_Standard Integer :: i,j Real :: Modul_Fx if (switch==0) then Ax(1:Nx,1:K)=EVtor(1:Nx,1:K) Fx(1:Nt,1:K) = MatMul(Y(1:Nt,1:Nx),Ax(1:Nx,1:K)) Do j=1,K !Fx(:,j)=Fx(:,j)/(EValues(j))/Factor ! Fx Chuan hoa Fx(:,j)=Fx(:,j)/Factor ! Voi PCs chuan hoa (Ax=EVtor) Enddo EOF(1:Nt,1:K) = Fx(1:Nt,1:K) PCA(1:Nx,1:K) = Ax(1:Nx,1:K) else if (switch==1) then Do j=1, K !Ax(:,j)=Ax(:,j)/sqrt(EValues(j)) Lamda(j,j)=Sqrt(EValues(j)) Enddo YT = Transpose(Y) !Fx(1:Nx,1:K) = MatMul(YT(1:Nx,1:Nt),Ax(1:Nt,1:K)) Fx(1:Nx,1:K) = MatMul(YT(1:Nx,1:Nt),EVtor(1:Nt,1:K)) Do j=1,K Modul_Fx=0 Do i=1, Ncount Modul_Fx = Modul_Fx + Fx(i,j)*Fx(i,j) Enddo Modul_Fx = Sqrt(Modul_Fx) Fx(:,j)=Fx(:,j)/Modul_Fx Enddo !Fx=MatMul(Fx(1:Ncount,1:K),Lamda(1:K,1:K)) EOF(1:Nt,1:K) = Ax(1:Nt,1:K) PCA(1:Nx,1:K) = Fx(1:Ncount,1:K) endif End Subroutine Calc_EOF_PCA_Standard !--------------------------------------------------- Subroutine Calc_EOF_PCA Integer :: i,j Real :: Modul_Fx if (switch==0) then Fx(1:Nt,1:K) = MatMul(Y(1:Nt,1:Nx),Ax(1:Nx,1:K)) Do j=1,K Fx(:,j)=Fx(:,j)/(EValues(j))/Factor ! Fx Chuan hoa !Fx(:,j)=Fx(:,j)/Factor ! Voi PCs chuan hoa (Ax=EVtor) Enddo EOF(1:Nt,1:K) = Fx(1:Nt,1:K) PCA(1:Nx,1:K) = Ax(1:Nx,1:K) else if (switch==1) then Do j=1, K Ax(:,j)=Ax(:,j)/sqrt(EValues(j)) Lamda(j,j)=Sqrt(EValues(j)) Enddo YT = Transpose(Y) Fx(1:Nx,1:K) = MatMul(YT(1:Nx,1:Nt),Ax(1:Nt,1:K)) Do j=1,K Modul_Fx=0 Do i=1, Ncount Modul_Fx = Modul_Fx + Fx(i,j)*Fx(i,j) Enddo Modul_Fx = Sqrt(Modul_Fx) Fx(:,j)=Fx(:,j)/Modul_Fx Enddo Fx=MatMul(Fx(1:Ncount,1:K),Lamda(1:K,1:K)) EOF(1:Nt,1:K) = Ax(1:Nt,1:K) PCA(1:Nx,1:K) = Fx(1:Ncount,1:K) endif End Subroutine Calc_EOF_PCA !--------------------------------------------------- Subroutine Write_test !Open (3,file="axfx1.txt",status="unknown") print*,"PCA" Do i=1,Nx print "(20F10.2)", (PCA(i,j),j=1,K) ! Mcount) Enddo print*,"EOF" Do i=1, Nt print "(20F10.2)", (EOF(i,j),j=1,K) ! Mcount) Enddo print*,"AA" Do i=1, K !Mcount print "(20F10.2)", (AA(i,j),j=1,K) ! Mcount) Enddo print*,"FF" Do i=1, K print "(20F10.2)", (FF(i,j),j=1,K)! Mcount) Enddo print*,"Y" Do i=1, Ncount print "(20F10.2)", (Y(i,j),j=1, Mcount) Enddo print*,"X1" Do i=1, Ncount print "(20F10.2)", (X1(i,j),j=1, Mcount) Enddo print*,"Y/X1" Do i=1, Ncount print "(20F10.2)", (YdX(i,j),j=1, Mcount) Enddo print*,"EVtor2" Do i=1, Mcount print "(20F10.2)", (EVtor2(i,j),j=1, Mcount) Enddo End Subroutine Write_test !--------------------------------------------------- SUBROUTINE Covar(Nt,Nx, Miss) Real, Allocatable :: R(:,:) Integer :: Nt, Nx, N1, L, J, K, IO, I Real :: Miss, Sx, Xtb !--------------------------------------------------- Allocate (R(Mcount,Mcount)) Do j=1, Nx Sx=0; Xtb=0; N1=0 Do i=1, Nt if (X(i,j)/=Miss) then n1=n1+1; Xtb=Xtb+X(i,j); Sx=Sx+X(i,j)*X(i,j) Endif Enddo if (n1>nt/2) then XBAR(j)=Xtb/real(n1); STD(j)=Sqrt(Sx/real(n1)-XBAR(j)*XBAR(j)) endif Do i=1, Nt if (X(i,j)==Miss) X(i,j)=XBAR(j) ! Fill-up Miss by mean values Y(i,j) = (X(i,j)-XBAR(j)) if (standard==1) Y(i,j) = Y(i,j)/STD(j) Enddo Enddo if (switch==1) then R = MatMul(Y,Transpose(Y)) else if (switch==0) then R = MatMul(Transpose(Y),Y) endif R(:,:)=R(:,:)/real(Nt) ! Store R in 1D-Array of Rx to be returned Do I=1,Mcount DO J=1,Mcount if ((I-J) < 0) then L=I+(J*J-J)/2 else L=J+(I*I-I)/2 endif RX(L)=R(i,j) ENDDO ENDDO Deallocate (R) End SUBROUTINE Covar !--------------------------------------------------- SUBROUTINE EIGEN(N,MV) ! Input: Covariance matrix of Rx(:) with real size is M*M elements ! switch: =1 neu tinh tuong quan theo cac lat cat thoi gian (dao chieu) ! Output: Rx with Eigenvalues are stored on diagonal ! V(:) matrix with real size is M*M elements stored Eigenvectors in column by column Real, Allocatable :: A(:), R(:) REAL ANORM,ANRMX,THR,X,Y,SINX,SINX2,COSX,COSX2,SINCS,XRANGE Integer :: N, I, J, M, icount WRITE(*,*)' BAT DAU EIGEN' icount=0 Allocate (R(Mcount*Mcount), A(Mcount*Mcount)) A(:) = Rx(:) 5 XRANGE=1.E-6 IF(MV-1) 10,25,10 10 IQ=-N DO 20 J=1,N IQ=IQ+N DO 20 I=1,N IJ=IQ+I R(IJ)=0. IF(I-J) 20,15,20 15 R(IJ)=1. 20 CONTINUE 25 ANORM=0. DO 35 I=1,N DO 35 J=1,N IF(I-J) 30,35,30 30 IA=I+(J*J-J)/2 ANORM=ANORM+A(IA)*A(IA) 35 CONTINUE IF(ANORM)165,165,40 40 ANORM=1.414*SQRT(ANORM) ANRMX=ANORM*XRANGE/FLOAT(N) IND=0 THR=ANORM 45 THR=THR/FLOAT(N) !icount=icount+1; print*,"Lan lap: ",icount 50 L=1 55 M=L+1 60 MQ=(M*M-M)/2 LQ=(L*L-L)/2 LM=L+MQ 62 IF(ABS(A(LM))-THR)130,65,65 65 IND=1 LL=L+LQ MM=M+MQ X=.5*(A(LL)-A(MM)) 68 Y=-A(LM)/SQRT(A(LM)*A(LM)+X*X) IF(X) 70,75,75 70 Y=-Y 75 SINX=Y/SQRT(2.*(1.+(SQRT(1.-Y*Y)))) SINX2=SINX*SINX 78 COSX=SQRT(1.-SINX2) COSX2=COSX*COSX SINCS=SINX*COSX ILQ=N*(L-1) IMQ=N*(M-1) DO 125 I=1,N IQ=(I*I-I)/2 IF(I-L) 80,115,80 80 IF(I-M) 85,115,90 85 IM=I+MQ GOTO 95 90 IM=M+IQ 95 IF(I-L) 100,105,105 100 IL=I+LQ GOTO 110 105 IL=L+IQ 110 X=A(IL)*COSX-A(IM)*SINX A(IM)=A(IL)*SINX+A(IM)*COSX A(IL)=X 115 IF(MV-1) 120,125,120 120 ILR=ILQ+I IMR=IMQ+I X=R(ILR)*COSX-R(IMR)*SINX R(IMR)=R(ILR)*SINX+R(IMR)*COSX R(ILR)=X 125 CONTINUE X=2.*A(LM)*SINCS Y=A(LL)*COSX2+A(MM)*SINX2-X X=A(LL)*SINX2+A(MM)*COSX2+X A(LM)=(A(LL)-A(MM))*SINCS+A(LM)*(COSX2-SINX2) A(LL)=Y A(MM)=X 130 IF(M-N) 135,140,135 135 M=M+1 GOTO 60 140 IF(L-(N-1)) 145,150,145 145 L=L+1 GOTO 55 150 IF(IND-1) 160,155,160 155 IND=0 GOTO 50 160 IF(THR-ANRMX) 165,165,45 165 IQ=-N DO 185 I=1,N IQ=IQ+N LL=I+(I*I-I)/2 JQ=N*(I-2) DO 185 J=I,N JQ=JQ+N MM=J+(J*J-J)/2 IF(A(LL)-A(MM)) 170,185,185 170 X=A(LL) A(LL)=A(MM) A(MM)=X IF(MV-1) 175,185,175 175 DO 180 K=1,N ILR=IQ+K IMR=JQ+K X=R(ILR) R(ILR)=R(IMR) 180 R(IMR)=X 185 CONTINUE V(:)=R(:) Rx(:) = A(:) Deallocate (A,R) WRITE(*,*)' KET THUC EIGEN' END SUBROUTINE EIGEN !--------------------------------------------------- ! CHUONG TRINH TIM K TRI RIENG SUBROUTINE TRACE(M,CON,K) Real, Allocatable :: R(:) Integer :: M, K, I, L Real :: CON WRITE(*,*)' BAT DAU C.T TRACE' Allocate (R(Mcount*Mcount)) R(:) = Rx(:) FM=M L=0 DO I=1,M L=L+I D(I)=R(L) ENDDO K=0 DO I=1,M IF (D(I) < CON) Exit K=K+1 D(I)=D(I)/FM ENDDO DO I=2,K D(I)=D(I)+D(I-1) ENDDO Deallocate (R) END SUBROUTINE TRACE !----------------------------------------------------------------------- !CHUONG TRINH TINH TAI TRONG NHAN TO SUBROUTINE LOAD(M,K) Real, Allocatable :: R(:) !, V(:) Real :: SQ Integer :: M, K, I, J, L Allocate (R(Mcount*Mcount)) R(:)=Rx(:) L=0 JJ=0 DO J=1,K JJ=JJ+J SQ=SQRT(R(JJ)) DO I=1,M L=L+1 V(L)=SQ*V(L) ENDDO ENDDO Deallocate (R) END SUBROUTINE LOAD !----------------------------------------------------------------------- END