C ================================================================== C FULLY AUTOMATIC REGRESSION MODELLING (FARM) TOOL C ------------------------------------------------------------------ C (C) ALEKSEI PARNOWSKI, SPACE RESEARCH INSTITUTE, KYIV, UKRAINE C CONTACT E-MAIL: PARNOWSKI@IKD.KIEV.UA C ------------------------------------------------------------------ C NRT FORECASTING PROCESSOR BETA VERSION 0.2 2013-02-02 C ------------------------------------------------------------------ C THIS PROGRAM IS SUPPOSED TO FORECAST THE OUTPUT IN HANDS-OFF MODE C ------------------------------------------------------------------ C INPUT FILES: C CONFIG FILE 'FARM_NRT.CFG' C NRT INPUT DATA, FILENAMES SET IN CONFIG FILE: C ACE/MAG FROM THE BEGINNING OF PREVIOUS MONTH C ACE/SWEPAM FROM THE BEGINNING OF PREVIOUS MONTH C KP PRESENT AND LAST MONTHS C DST PRESENT AND LAST MONTHS C MODEL STRUCTURE, FILENAME SET IN CONFIG FILE C COVARIANCE MATRIX, FILENAME SET IN CONFIG FILE C OUTPUT FILES: C FORECAST, FORMAT AND FILENAME SET IN CONFIG FILE C METADATA IN XML FORMAT, FILENAME SET IN CONFIG FILE C ------------------------------------------------------------------ C NOTE: C THIS VERSION LACKS THE GAP FILLING CAPABILITY AND CAN FORECAST C ONLY ONE PREDICTAND WITH SEVERAL LEAD TIME VALUES AT A TIME C ================================================================== INCLUDE 'readnrtd.for' INCLUDE 'forecast.for' IMPLICIT NONE C K ENUMERATES INPUTS, M - REGRESSORS, N - DATAPOINTS C T IS TIME INTEGER, PARAMETER :: KIN=12 INTEGER, PARAMETER :: KMAX=KIN+4 CHARACTER*255, PARAMETER :: FCFG='farm_nrt.cfg' INTEGER LMAX(KMAX) INTEGER, ALLOCATABLE :: R(:,:,:),RN(:) INTEGER I,J,K,L,M,N,MAXL,TMAX,DAT(8) C U - INPUTS, UF - FILLED VALUES, R - REGRESSORS, C - COEFFICIENTS C DOUBLE PRECISION IS NECESSARY REAL*8, ALLOCATABLE :: U(:,:),C(:),Y(:),DYM(:),COV(:,:) REAL*8, ALLOCATABLE :: YY(:,:),DYY(:,:) INTEGER, ALLOCATABLE :: YR(:),MO(:),DA(:),HR(:),DOY(:) C PARAMETERS READ FROM THE CONFIG FILE: C KY - NUMBER OF OUTPUT AMONG INPUTS, DT - CADENCE, LT - LEAD TIME INTEGER KY,DT(KMAX),LT,UFLAG(KMAX),NY,NU,NP,NUY,NUP,ST REAL*8 UF(KMAX),MU,MY,MP,DU,DY,DP,RSSY,RSSP,SDU,SDY,SDP, &PY,PP,RY,RP,SS CHARACTER*255 FNRES(12),FNCOV(12),FNDAT,FNMET,FMDAT,FMDATH,FNMAG, &FNSWEPAM,FNKP,FNDST,UN(KIN),UU(KIN) C SET MAXIMAL LAGS AND FILL VALUES INDEPENDENTLY FOR EACH INPUT NAMELIST /FARM_FOR/ KY,DT,LT,LMAX,UF,UFLAG,UN,UU, &FNRES,FNCOV,FMDAT,FMDATH,FNMAG,FNSWEPAM,FNKP,FNDST,FNDAT,FNMET C READ CONFIG FROM FILE OPEN(1,FILE=FCFG,STATUS='OLD',ACTION='READ',DELIM='APOSTROPHE') READ(1,FARM_FOR) CLOSE(1) C DETERMINE MAXIMUM LAG AND DO SOME CHECKS IF((KY.LT.1).OR.(KY.GT.KIN))GOTO 111 C ADD CHECKS FOR COINCIDING NRT DATA FILENAMES HERE, LABEL 333 C ADD CHECKS FOR COINCIDING MODEL AND COV FILENAMES HERE, LABEL 444 MAXL=-1 DO 1 K=1,KMAX IF(DT(K).LE.0)GOTO 666 IF(LMAX(K).LT.0)GOTO 777 IF(UFLAG(K).EQ.0)GOTO 1 IF(LMAX(K).GT.MAXL)MAXL=LMAX(K) 1 CONTINUE IF(MAXL.LT.0)GOTO 1111 C DETERMINE THE NUMBER OF DATAPOINTS OPEN(1,FILE=FNMAG,STATUS='OLD',ACTION='READ') N=0 READ(1,'(3(/))') 2 READ(1,'(I4)',END=3)I N=N+1 GOTO 2 3 CLOSE(1) TMAX=CEILING((N+LT)/24.)*24 C ASSIGN MEMORY FOR DYNAMIC ARRAYS ALLOCATE(YR(TMAX),MO(TMAX),DA(TMAX),DOY(TMAX),HR(TMAX), &Y(TMAX),DYM(TMAX),YY(LT/DT(KY),TMAX),DYY(LT/DT(KY),TMAX), &U(KMAX,TMAX)) C ------------------------------------------------------------------ C INPUT NRT DATA CALL READNRTD(N,TMAX,YR,MO,DA,DOY,HR,U,UF,KMAX, &FNMAG,FNSWEPAM,FNKP,FNDST) IF(N.LE.MAXL)GOTO 888 C ------------------------------------------------------------------ C GET CURRENT DATE AND TIME CALL DATE_AND_TIME(VALUES = DAT) IF(DAT(1).LT.0)GOTO 1333 IF(DAT(2).LT.0)GOTO 1333 IF(DAT(3).LT.0)GOTO 1333 IF(DAT(5).LT.0)GOTO 1333 IF(DAT(6).LT.0)GOTO 1333 IF(DAT(7).LT.0)GOTO 1333 OPEN(2,FILE=FNMET) WRITE(2,900)DAT(1),DAT(2),DAT(3),DAT(5),DAT(6),DAT(7), &YR(MAXL+LT+1),MO(MAXL+LT+1),DA(MAXL+LT+1),HR(MAXL+LT+1), &YR(N),MO(N),DA(N),HR(N) 900 FORMAT(''/, &''/, &1X,'Geomagnetic forecast module'/, &1X,'diverse global'/, &1X,'L2'/, &1X,'AFFECTS'/, &1X,'', &I4.4,2('-',I2.2),'T',I2.2,2(':',I2.2), &''/, &1X,'1.0'/, &1X,'', &I4.4,2('-',I2.2),'T',I2.2,':00:00', &''/, &1X,'', &I4.4,2('-',I2.2),'T',I2.2,':59:59', &''/, &1X,'global'/, &1X,'ASCII'/, &1X,'nc'//, &1X,'', &'Lorem ipsum dolor sit amet, consectetuer adipiscing elit, ', &'sed diam nonummy nibh euismod tincidunt ut laoreet dolore ', &'magna aliquam erat volutpat. Ut wisi enim ad minim veniam, ', &'quis nostrud exerci tation ullamcorper suscipit lobortis nisl ', &'ut aliquip ex ea commodo consequat. Duis autem vel eum iriure ', &'dolor in hendrerit in vulputate velit esse molestie consequat, ', &'vel illum dolore eu feugiat nulla facilisis at vero eros et ', &'accumsan et iusto odio dignissim qui blandit praesent luptatum ', &'zzril delenit augue duis dolore te feugait nulla facilisi. Nam ', &'liber tempor cum soluta nobis eleifend option congue nihil ', &'imperdiet doming id quod mazim placerat facer possim assum. ', &'Typi non habent claritatem insitam; est usus legentis in iis ', &'qui facit eorum claritatem. Investigationes demonstraverunt ', &'lectores legere me lius quod ii legunt saepius. Claritas est ', &'etiam processus dynamicus, qui sequitur mutationem ', &'consuetudium lectorum. Mirum est notare quam littera gothica, ', &'quam nunc putamus parum claram, anteposuerit litterarum ', &'formas humanitatis per seacula quarta decima et quinta decima. ', &'Eodem modo typi, qui nunc nobis videntur parum clari, fiant ', &'sollemnes in futurum.', &''//, &1X,'', &'Aleksei Parnowski: parnowski@ikd.kiev.ua', &''/, &1X,'Space Research Institute, ', &'Prospekt Akademika Glushkova 40 building 4/1, ', &'03680 MSP Kyiv-187, Ukraine', &''/, &1X,'', &'0.1', &''/) WRITE(2,950)TRIM(UN(KY)),TRIM(UU(KY)),DT(KY),UF(KY) 950 FORMAT(1X,''/, &2X,'',A6,''/, &2X,'',A6,''/, &2X,'',I2,''/, &2X,'',F7.1,''/, &2X,'') DO 17 L=1,LT/DT(KY) C DETERMINE THE NUMBER OF REGRESSORS OPEN(1,FILE=FNRES(L),STATUS='OLD',ACTION='READ') M=0 6 READ(1,'()',END=7) M=M+1 GOTO 6 7 REWIND(1) ALLOCATE(C(M),COV(M,M),R(M,KMAX,3),RN(M)) C INPUT THE MODELS DO 8 I=1,M 8 READ(1,100,ERR=999)C(I),RN(I),((R(I,K,J), J=1,3), K=1,RN(I)) 100 FORMAT(5X,D12.6,42X,I1,9(1X,I2,1X,I3,1X,I1)) CLOSE(1) C INPUT THE COVARIANCE MATRIX OPEN(1,FILE=FNCOV(L),STATUS='OLD',ACTION='READ') DO 9 I=1,M 9 READ(1,200,ERR=1222)(COV(I,J), J=1,M) 200 FORMAT(999D13.6) CLOSE(1) CALL FORECAST(M,N,TMAX,KY,KMAX,R,RN,U,UF,DT,LT,MAXL,C,Y,COV,DYM) DEALLOCATE(C,COV,R,RN) DO 10 I=1,TMAX YY(L,I)=Y(I) DYY(L,I)=DYM(I) 10 CONTINUE C OUTPUT RESULTS AND CALCULATE METRICS MU=0._8 MY=0._8 MP=0._8 RSSY=0._8 RSSP=0._8 NU=0 NY=0 NP=0 NUY=0 NUP=0 ST=0 IF(LT.EQ.0)ST=1 DO 13 I=MAXL+LT+1,N+L,DT(KY) IF(U(KY,I).EQ.UF(KY))GOTO 11 NU=NU+1 MU=MU+U(KY,I) 11 IF(U(KY,I-L-ST).EQ.UF(KY))GOTO 12 NP=NP+1 MP=MP+U(KY,I-L-ST) IF(U(KY,I).EQ.UF(KY))GOTO 12 NUP=NUP+1 RSSP=RSSP+(U(KY,I-L-ST)-U(KY,I))**2 12 IF(Y(I).EQ.UF(KY))GOTO 13 NY=NY+1 MY=MY+Y(I) IF(U(KY,I).EQ.UF(KY))GOTO 13 NUY=NUY+1 RSSY=RSSY+(Y(I)-U(KY,I))**2 13 CONTINUE MU=MU/NU IF(NY.NE.0)MY=MY/NY IF(NP.NE.0)MP=MP/NP IF(NY.NE.0)RSSY=RSSY/NUY IF(NP.NE.0)RSSP=RSSP/NUP C NOTE THERE IS NO -M IN THE DENOMINATOR C CALCULATE LINEAR CORRELATION RY=0._8 RP=0._8 SDU=0._8 SDY=0._8 SDP=0._8 DO 16 I=MAXL+LT+1,N+L,DT(KY) IF(U(KY,I).EQ.UF(KY))GOTO 14 DU=U(KY,I)-MU SDU=SDU+DU**2 14 IF(NP.EQ.0)GOTO 15 IF(U(KY,I-L-ST).EQ.UF(KY))GOTO 15 DP=U(KY,I-L-ST)-MP SDP=SDP+DP**2 IF(U(KY,I).EQ.UF(KY))GOTO 15 RP=RP+DU*DP 15 IF(NY.EQ.0)GOTO 16 IF(Y(I).EQ.UF(KY))GOTO 16 DY=Y(I)-MY SDY=SDY+DY**2 IF(U(KY,I).EQ.UF(KY))GOTO 16 RY=RY+DU*DY 16 CONTINUE PY=-.999_8 PP=-.999_8 SS=-.999_8 IF(NP.NE.0)RP=RP/DSQRT(SDU)/DSQRT(SDP) IF(NY.NE.0)RY=RY/DSQRT(SDU)/DSQRT(SDY) SDU=SDU/NU IF(NP.NE.0)PP=1._8-RSSP/SDU IF(NY.NE.0)PY=1._8-RSSY/SDU SDU=DSQRT(SDU) IF(NY.GT.0)SDY=SDY/NY IF(NP.GT.0)SDP=SDP/NP IF((NP.NE.0).AND.(NY.NE.0))SS=1._8-RSSY/RSSP IF(NY.NE.0)RSSY=DSQRT(RSSY) IF(NP.NE.0)RSSP=DSQRT(RSSP) IF(NP.NE.0)SDP=DSQRT(SDP) IF(NY.NE.0)SDY=DSQRT(SDY) C OUTPUT METADATA WRITE(2,300)L,RSSY,RSSP,PY,PP,RY,RP,SS 300 FORMAT(3X,''/, &4X,'',I2,''/, &4X,'' &5X,''/, &6X,'','Mean square error',''/, &6X,'',F6.3,''/, &6X,'',F6.3,''/, &5X,''/, &5X,''/, &6X,'','Prediction efficiency',''/, &6X,'',F6.3,''/, &6X,'',F6.3,''/, &5X,''/, &5X,''/, &6X,'','Pearson''s linear correlation coefficient', &''/, &6X,'',F6.3,''/, &6X,'',F6.3,''/, &5X,''/, &5X,''/, &6X,'','Skill score',''/, &6X,'',F6.3,''/, &5X,''/, &4X,''/, &3X,'') 17 CONTINUE WRITE(2,'(2X,""/, &1X,""/, &"")') OPEN(1,FILE=FNDAT) WRITE(1,FMDATH)TRIM(UN(KY))//' ['//TRIM(UU(KY))//']', &(I,'['//TRIM(UU(KY))//']',I=1,LT,DT(KY)) DO 18 I=MAXL+LT+1,TMAX,DT(KY) 18 WRITE(1,FMDAT)(24*(DOY(I)-1)+HR(I)+1),TMAX-I+1, &YR(I),MO(I),DA(I),DOY(I),HR(I),U(KY,I), &(YY(L,I),L=1,LT/DT(KY)) CLOSE(1) CLOSE(2) GOTO 1000 C ------------------------------------------------------------------ C ERROR HANDLING (TO BE REPLACED BY SOMETHING MORE SOPHISTICATED) 111 WRITE(*,'(1X,"KY MUST BE WITHIN 1..",I2," RANGE, YOU PROVIDED ", &I2)')KIN,KY GOTO 1000 333 WRITE(*,*)'INPUT AND OUTPUT FILES CANNOT BE THE SAME' GOTO 1000 444 WRITE(*,*)'DIFFERENT OUTPUT FILES CANNOT BE THE SAME' GOTO 1000 666 WRITE(*,*)'CADENCE MUST BE POSITIVE' GOTO 1000 777 WRITE(*,*)'MAXIMUM LAG MUST BE 0 OR POSITIVE' GOTO 1000 888 WRITE(*,*)'NOT ENOUGH DATA IN VALIDATION SAMPLE' GOTO 1000 999 WRITE(*,'(1X,"ERROR READING MODEL AT LINE #",I4)')I GOTO 1000 1111 WRITE(*,*)'AT LEAST ONE OF UFLAG VALUES MUST BE NON-ZERO' GOTO 1000 1222 WRITE(*,'(1X,"ERROR READING COVARIANCE MATRIX AT",2I4)')I,J GOTO 1000 1333 WRITE(*,*)'ERROR GETTING SYSTEM TIME' GOTO 1000 C RELEASE MEMORY 1000 DEALLOCATE(YR,MO,DA,DOY,HR,U,Y,DYM,YY,DYY) END