Use NetCDF Use Logistic_Regression Real, parameter :: Miss=-99.0 Integer, Parameter :: Yr1=1983, Yr2=2020 Integer, Parameter :: Nmodel=2, Nmon=12, Nreg=7, Nlead=6, Dr_Lev=3, Nlen=3, Nidx=2 Integer, Parameter :: Ncase=3 ! "mod","cli","ann","mlr" Character (Len=10) :: Model_nm(Nmodel)=(/"RegCM","clWRF"/) Character (Len=10) :: IDXname (Nidx) =(/"spei", "spi"/) Integer, Parameter :: Len_SPI (Nlen) =(/1,3,6/) !Character (Len=3), Parameter :: Cases(Ncase)=(/"mod","cli","ann","mlr"/) Character (Len=3), Parameter :: Cases(Ncase)=(/"cli","ann","mlr"/) !Character (Len=3), Parameter :: Cases(Ncase)=(/"ann","mlr"/) !----------------------------------------------------------------------------- Real, Allocatable :: Regs(:,:), Lon(:), Lat(:), xlon(:), xlat(:) Real, Allocatable :: SPEI_M1 (:,:,:,:,:,:) ! Nlon,Nlat,Nlead,Ntime,Ncase,Nmodel) - Mem 1 only Real, Allocatable :: SPEI_est(:,:,:,:) ! Nlon,Nlat,Nlead,Ntime Real, Allocatable :: SPEI_obs(:,:,:) ! Nlon,Nlat,Ntime Real, Allocatable :: Ymax(:,:,:), Ymin(:,:,:) ! Nlon,Nlat,Nlead Real, Allocatable :: Yest(:) ! Ntime,Myvar Integer :: Nlon, Nlat, Ntime, Nt !----------------------------------------------------------------------------- Integer :: ireg, ilead !Real :: a,b,c,d,Thre Real :: PC0,POD0,FAR0,POFD0, FBI0 !----------------------------------------------------------------------------- Real,Allocatable :: XX(:,:) Real (dp), Allocatable :: X(:,:), beta(:), se_beta(:) Integer, Allocatable :: N(:),s1(:),LG(:) Real (dp) :: tmp, chisq, devnce, tail_prob Logical :: Status !----------------------------------------------------------------------------- Character (Len=255) :: Dir_ANN_Save, Inp_Dir, Out_Dir, Mask_File, Fname !Namelist /ann_info/ Nlayer,learn_rate,Nepoch !Namelist /layer_info/ Layer_size !Namelist /ann_hai/ Nhidden,learning_rate,momentum,training_mode,Nepoc,debug !----------------------------------------------------------------------------- Inp_Dir = "/home/tanpv/Data/562_VNU/0.25_VN/Dr_IDXs_Ope/" Mask_File = "/home/tanpv/Data/AWC/Mask_VN_0.25_Corrected.nc" Out_Dir = "/home/tanpv/Progs/562_VNU/2024/Ann_New/Dr_Ann_Out/" Dir_ANN_Save = "/home/tanpv/Progs/562_VNU/2024/Ann_New/ANN_Save/" Ntime = (Yr2-Yr1+1)*Nmon !----------------------------------------------------------------------------- !Open (1, file="Namelist_Ann_Dr.txt",Status="old") !Read (1, ann_info) !Allocate (Layer_size(Nlayer)) !Read (1, layer_info) !Read (1, ann_hai) !Close (1) Mxvar=Nmodel*Ncase; Myvar=1 !----------------------------------------------------------------------------- Call rd_ncdf_2d (Mask_File,Lon,Lat,Nlon,Nlat,Regs,"reg") !----------------------------------------------------------------------------- Do idx=1,Nidx ! "spi", "spei" Do ilen=2,Nlen ! 1, 3, 6 Call Read_all_Dr_IDXs (trim(IDXname(idx)),ilen) if (Allocated (SPEI_est)) Deallocate (SPEI_est) ! Nlon,Nlat,Nlead,Ntime Allocate (SPEI_est(Nlon,Nlat,Nlead,Ntime)); SPEI_est=Miss Do ilead=1,Nlead Do ilon=1,Nlon Do ilat=1,Nlat if (Regs(ilon,ilat)==2) then Call Calc_Logistic (Dir_ANN_Save,idx,ilen,ilead,ilon,ilat) pause endif Enddo Enddo Enddo Enddo Enddo !----------------------------------------------------------------------------- CONTAINS !-------------------------------------------------------------------------------- Subroutine Calc_Logistic (Dir_ANN_Save,idx,ilen,ilead,ilon,ilat) Character (Len=*) :: Dir_ANN_Save Character (Len=255) :: save_file Integer :: idx,ilen,ilead,ilon,ilat,ispl, i,j,k,icount Integer :: n0, n1, Nspl1 Real :: XX(Ntime,Mxvar),YY(Ntime),xtmp(Mxvar), ytmp Real :: X_diff(Mxvar), X_new(Mxvar) Double precision :: x1d (Mxvar), y1d (Myvar) Real, Allocatable :: X1(:,:), X11(:,:) Integer, Allocatable :: nbors(:),LG1(:) Logical :: OK Integer, Allocatable :: Check(:) Real :: Cor, Prob, lambda if (allocated (LG)) Deallocate (LG); if (allocated (X)) Deallocate (X) if (allocated (N)) Deallocate (N); if (allocated (s1)) Deallocate (s1) if (allocated (beta)) Deallocate (beta) if (allocated (se_beta)) Deallocate (se_beta) if (allocated (Yest)) Deallocate (Yest) if (Allocated (Check)) Deallocate (Check) Allocate (Check(Ntime)) Do i=1,Ntime YY(i)=SPEI_obs(ilon,ilat,i) j=0 Do imod=1,Nmodel Do icas=1,Ncase j=j+1; XX(i,j)=SPEI_M1(ilon,ilat,ilead,i,icas,imod) Enddo Enddo Enddo ! icount=0; Check(:)=0 Do i=1,Ntime ytmp=YY(i); xtmp(1:Mxvar)=XX(i,1:Mxvar) OK=Valid_Ar1(xtmp,Mxvar,Miss).AND.(ytmp/=Miss) if (OK) then icount=icount+1; Check(i)=1 endif Enddo Nspl = icount Allocate (LG(Nspl),X1(Nspl,Mxvar),s1(Nspl)) ! N(Nspl), s1(Nspl)) Allocate (beta(0:Mxvar),se_beta(0:Mxvar)) Allocate (Yest(Nspl)) !Allocate (X11(Nspl,Mxvar)) !N(:) = 1 ! To be used by Logistic !-------------------------------------------------------- Thre=-1.0; Prob0=0.6 !-------------------------------------------------------- !Allocate (LG1(Nspl1), X1(Nspl1,Mxvar),s1(Nspl1),Yest(Nspl1)) ispl=0 Do i=1,Ntime if (Check(i)==1) then ispl=ispl+1 if (YY(i)<=Thre) then LG(ispl) = 1 else LG(ispl) = 0 endif X1(ispl,1:Mxvar) = XX(i,1:Mxvar) ! Origin !X (ispl,1:Mxvar) = XX(i,1:Mxvar) ! Real (DP) endif Enddo !Allocate (X1(Nspl,Mxvar)); X1(:,:)=X(:,:) ! For oversampling !Allocate (nbors(Nspl)) !-------------------------------------------------------- n0 = count(LG==0); n1 = count(LG==1); n_over = n0 - n1; n_nbor = 5 Nspl1 = Nspl + n_over ! Sample size Allocate (LG1(Nspl1), X(Nspl1,Mxvar), N(Nspl1)) !, s1(Nspl)) Allocate (nbors(Nspl)) !Allocate (beta(0:Mxvar),se_beta(0:Mxvar)) !Allocate (Yest(Nspl)) N(:) = 1 X (1:Nspl,Mxvar) = X1(1:Nspl,Mxvar) ! Real to Real(dp) LG1(1:Nspl) = LG(1:Nspl) !-------------------------------------------------------- ! Oversampling for LG=1 !print*,n0,n1,n_over,Nspl,Nspl1; pause i1 = 0 Do j=1,n_over i = i1 Do !i = random_int(1, Nspl) i = i + 1 !print*,i;pause if (i > Nspl) exit if (LG(i) == 1) then xtmp(:) = X1(i,:) i1 = i exit endif Enddo !xtmp(:) = X1(i,:) Call Nearest_point (xtmp,X1,LG,Nspl,Mxvar,n_nbor,nbors) k = random_int(1, n_nbor) !print*,"k=",k X_diff(:) = X1(nbors(k),:) - X1(i,:) !Call random_number (lambda) lambda = rand() X_new(:) = X1(i,:) + lambda * X_diff(:) !X_new(:) = X1(nbors(k),:) X(Nspl+j,:) = X_new(:) LG1(Nspl+j) = LG(i) print*,i,LG1(i),LG1(Nspl+j) !nbors(k)) print"(10F7.3)",X1(i,:) print"(10F7.3)",X_new(:) !pause Enddo !-------------------------------------------------------- CALL logistic(Nspl1,X,Mxvar,LG1,N, chisq, devnce, ndf, beta, se_beta, ier) !-------------------------------------------------------- Do i=1,Nspl ytmp=beta(0) Do j=1,Mxvar ytmp=ytmp+beta(j)*X1(i,j) ! F = exp(b0 + b1.X1 + b2.X2 + ... + bk.Xk) Enddo ytmp=EXP(ytmp); ytmp=ytmp/(1.0+ytmp); ! p = F/(1 + F) Yest(i) = ytmp !if (ytmp>0.30) then; s1(i)=1; else; s1(i)=0; endif Enddo Where (Yest>=Prob0) s1 = 1 Elsewhere s1 = 0 Endwhere !-------------------------------------------------------- ispl=0 Do i=1,Ntime if (Check(i)==1) then ispl=ispl+1 print"(A,2F8.4,2I4)","Y_inp,Yest,s1,LG: ",YY(i),Yest(ispl),s1(ispl),LG(ispl) endif Enddo print"(20F8.3)",beta(:) !print*,Layer_size ! Thre=-0.5 ! Call Correlation(Yest(:,1),YY(:),Ntime,Cor,prob,Miss) ! print*,"------------------------------------------------------------------------------" Call Cate_Veri (LG(:),s1(:),Nspl,PC0,POD0,SR0,FAR0,POFD0,FBI0) print"(10A8)","PC","POD","SR","POFD","FAR","FBI" print"(10F8.3)",PC0,POD0,SR0,POFD0,FAR0,FBI0 !print*,RMSE (Yest(:,1),YY(:),Ntime,Miss) End Subroutine Calc_Logistic !------------------------------------------------------------------------------ integer function random_int(low, high) ! Generate a random real number between 0 and 1 real :: random_value Integer :: low, high Call random_number (random_value) ! Scale and offset to fit the desired integer range random_int = int( (high - low + 1) * random_value ) + low end function random_int !------------------------------------------------------------------------------ Subroutine Nearest_point (X0,X,Y,N,M,k,nb) ! Find a nearest position of k of X0 in X with Y is minority_class Integer :: N,M,k,i,j, temp_index Real :: X(N,M), X0(M) Real, Dimension (M) :: dat, diff Real, Dimension (N) :: dist Integer :: Y(N), nb(N) Real :: tmp Do i=1,N nb(i) = i ! Order of index if (Y(i)==1) then dat(:) = X(i,:) diff(:) = X0(:) - dat(:) !dist(i) = sqrt(sum(diff**2)) ! Euclidean distance tmp = 0. Do j=1,M tmp = tmp + diff(j)*diff(j) Enddo dist(i) = sqrt(tmp) else dist(i) = 999999.0 !real(huge) endif !print*,"x0=",X0(:) !print*,"diff=",diff(:) !print*,"dist=",dist(i) Enddo !pause ! Find k nearest neighbors do i = 1, N-1 do j = 1, N-i if (dist(j) > dist(j+1)) then temp_index = nb(j) nb(j) = nb(j+1) nb(j+1) = temp_index tmp = dist(j) dist(j) = dist(j+1) dist(j+1) = tmp endif enddo enddo ! Select the top k nearest neighbors (excluding empty data points) k = min(k, count(Y==1)) ! Ensure k doesn't exceed the number of valid data points do i = 1, k nb(i) = nb(N - i + 1) print*,nb(i),dist(i) enddo !pause End Subroutine Nearest_point !------------------------------------------------------------------------------ Subroutine Read_all_Dr_IDXs (IDX_nm,ilen) Character (Len=*) :: IDX_nm Character (Len=20) :: Var_nm Integer :: it,i,ilen,ilon,ilat,ilead,Nlo,Nla,Nme,Nt,Nt1,Nt2,icas, it1, it2, icpr, ict, icpet Real, Allocatable :: X4d(:,:,:,:),X3d(:,:,:),xlon(:),xlat(:),xlead(:),xtime(:) if (Allocated (SPEI_M1)) Deallocate (SPEI_M1) Allocate (SPEI_M1(Nlon,Nlat,Nlead,Ntime,Ncase,Nmodel)) write (Fname,"(A,I2.2,A)") trim(Inp_Dir)//"Dr_IDX_1983-2020_"//trim(IDX_nm)//"_",Len_SPI(ilen),".nc" Do imod=1,Nmodel Do icas=1,Ncase Var_nm=trim(Model_nm(imod))//"_"//trim(Cases(icas)) Call rd_ncdf_4d (Fname,xlon,xlat,xlead,xtime,Nlo,Nla,Nle,Nt,X4d,trim(Var_nm)) ! X4d(lon,lat,lead,time) Do ilead=1,Nlead Do it=1,Ntime SPEI_M1(:,:,ilead,it,icas,imod) = X4d(:,:,ilead,it) Enddo Enddo Deallocate (xlon,xlat,xlead,xtime,X4d) Enddo Enddo if (Allocated(SPEI_obs)) Deallocate (SPEI_obs) Allocate (SPEI_obs(Nlon,Nlat,Ntime)) Call rd_ncdf_3d (Fname,xlon,xlat,xtime,Nlo,Nla,Nt,X3d,"obs") SPEI_obs(:,:,1:Ntime) = X3d(:,:,1:Ntime) Deallocate (xlon,xlat,xtime,X3d) Where (IsNaN(SPEI_M1)) SPEI_M1=Miss Where (IsNaN(SPEI_obs)) SPEI_obs=Miss Where (SPEI_M1>10.0) SPEI_M1=Miss Where (SPEI_obs>10.0) SPEI_obs=Miss End Subroutine Read_all_Dr_IDXs !------------------------------------------------------------------------------ Subroutine Cate_Veri (O,F,N,PC0,POD0,SR0,FAR0,POFD0,FBI0) Integer :: N,i Integer :: O(N),F(N) Real :: Thre,PC0,POD0,SR0,FAR0,POFD0,FBI0, a,b,c,d, Thre1 !Where (O/=Miss.AND.O>Thre1) O=0 !Where (F/=Miss.AND.F>Thre) F=0 !Where (O/=Miss.AND.O<=Thre1) O=1 !Where (F/=Miss.AND.F<=Thre) F=1 a=0; b=0; c=0; d=0 Do i=1,N if (O(i)==1.AND.F(i)==1) a=a+1 if (O(i)==1.AND.F(i)==0) c=c+1 if (O(i)==0.AND.F(i)==1) b=b+1 if (O(i)==0.AND.F(i)==0) d=d+1 Enddo print*,int(a),Int(b) print*,int(c),Int(d) PC0 =PC (a,b,c,d) POD0 =POD (a,b,c,d) SR0 =SR (a,b,c,d) FAR0 =FAR (a,b,c,d) POFD0=POFD(a,b,c,d) FBI0 =FBI (a,b,c,d) End Subroutine Cate_Veri !------------------------------------------------------------------------------ Include "/home/tanpv/Mylibs/sub_rd_ncdf_Tan.F90" Include "/home/tanpv/Mylibs/Function_Valid_Ar1.f90" Include "/home/tanpv/Mylibs/Correlation.f90" Include "/home/tanpv/Mylibs/Sub_Verification.f90" !Include "/home/tanpv/Mylibs/Sub_neural_2024.f90" !------------------------------------------------------------------------------ End Include "/home/tanpv/Mylibs/Student_Prob.f90"