149 lines
4.4 KiB
FortranFixed
149 lines
4.4 KiB
FortranFixed
!This subroutine fills gaps in the y variable based on neural network regression
|
|
subroutine gapfilling(nx,nobs,xsamp,ysamp,nmax,maxdist)
|
|
implicit none
|
|
!Gaps must be represented by -9999
|
|
!It is ok to have missing values in xsamp. If any dimension in xsamp is missing, that dimension
|
|
!is not used as the independent variable for the gap in y. For different gaps in y, the dimensions used
|
|
!in x may be different.
|
|
!If a gap is less than maxdist points away from a previous gap for which a fit was conducted and if
|
|
!the dimension is the same, then the previous fit is used (a new fit is not conducted)
|
|
integer nx,nobs,nmax
|
|
double precision xsamp(1:nobs,1:nx),ysamp(1:nobs)
|
|
!Locals
|
|
integer i,j,k,n,idowhat,nh,nxfit,nobsfit,negobsfit,ixuse(nx),
|
|
&ipregap,iposfit(-nmax:nmax),iuseit,ixuse_pre(nx),maxdist
|
|
parameter(nh=15)
|
|
double precision w(1:nx,1:nh),bph(nh),q(nh),bend,
|
|
&xnew(nx),calvalue(nobs),fn9999,tiny,xfit(-nmax:nmax,nx),
|
|
&yfit(-nmax:nmax),rsq,ysamppred(nobs)
|
|
parameter(fn9999=-9999.0d0,tiny=1.0d-6)
|
|
!
|
|
bend=fn9999
|
|
ipregap=-100000
|
|
do i=1,nobs
|
|
ysamppred(i)=ysamp(i)
|
|
if(dabs(ysamp(i)-fn9999).gt.tiny)goto 1000
|
|
!a gap
|
|
nxfit=0
|
|
do j=1,nx
|
|
if(dabs(xsamp(i,j)-fn9999).lt.tiny)then
|
|
!this x dimension is not used
|
|
ixuse(j)=0
|
|
else
|
|
!this x dimension is used
|
|
nxfit=nxfit+1
|
|
xnew(nxfit)=xsamp(i,j)
|
|
ixuse(j)=1
|
|
endif
|
|
enddo
|
|
if(nxfit.eq.0)goto 1000
|
|
if((i-ipregap).lt.maxdist)then
|
|
iuseit=1
|
|
do j=1,nx
|
|
if(ixuse(j).ne.ixuse_pre(j))iuseit=0
|
|
enddo
|
|
if(iuseit.eq.1)goto 30
|
|
endif
|
|
|
|
!Fill this gap by choosing the nmax valid points that are closest to i for the fitting
|
|
!First pick up nmax points from the lower side, index positive
|
|
nobsfit=0
|
|
n=i-1
|
|
1 if(n.lt.1)goto 2
|
|
if(dabs(ysamp(n)-fn9999).gt.tiny)then
|
|
iuseit=1
|
|
!make sure it has the x dimensions needed
|
|
do j=1,nx
|
|
if(ixuse(j).eq.1)then
|
|
if(dabs(xsamp(n,j)-fn9999).lt.tiny)iuseit=0
|
|
endif
|
|
enddo
|
|
if(iuseit.eq.1)then
|
|
nobsfit=nobsfit+1
|
|
yfit(nobsfit)=ysamp(n)
|
|
iposfit(nobsfit)=n
|
|
k=0
|
|
do j=1,nx
|
|
if(ixuse(j).eq.1)then
|
|
k=k+1
|
|
xfit(nobsfit,k)=xsamp(n,j)
|
|
endif
|
|
enddo
|
|
endif
|
|
endif
|
|
if(nobsfit.lt.nmax)then
|
|
n=n-1
|
|
goto 1
|
|
endif
|
|
!
|
|
!now pick up nmax points form the higher side, index negative
|
|
2 negobsfit=1
|
|
n=i+1
|
|
3 if(n.gt.nobs)goto 4
|
|
if(dabs(ysamp(n)-fn9999).gt.tiny)then
|
|
iuseit=1
|
|
!make sure it has the x dimensions needed
|
|
do j=1,nx
|
|
if(ixuse(j).eq.1)then
|
|
if(dabs(xsamp(n,j)-fn9999).lt.tiny)iuseit=0
|
|
endif
|
|
enddo
|
|
if(iuseit.eq.1)then
|
|
negobsfit=negobsfit-1
|
|
yfit(negobsfit)=ysamp(n)
|
|
iposfit(negobsfit)=n
|
|
k=0
|
|
do j=1,nx
|
|
if(ixuse(j).eq.1)then
|
|
k=k+1
|
|
xfit(negobsfit,k)=xsamp(n,j)
|
|
endif
|
|
enddo
|
|
endif
|
|
endif
|
|
if(negobsfit.gt.-nmax)then
|
|
n=n+1
|
|
goto 3
|
|
endif
|
|
|
|
!finally pick up the nmax closest points
|
|
4 if((nobsfit-negobsfit+1).le.nmax)goto 10
|
|
if((i-iposfit(nobsfit)).gt.(iposfit(negobsfit)-i))then
|
|
nobsfit=nobsfit-1
|
|
else
|
|
negobsfit=negobsfit+1
|
|
endif
|
|
goto 4
|
|
10 do n=negobsfit,0
|
|
nobsfit=nobsfit+1
|
|
yfit(nobsfit)=yfit(n)
|
|
do j=1,nxfit
|
|
xfit(nobsfit,j)=xfit(n,j)
|
|
enddo
|
|
enddo
|
|
idowhat=1
|
|
call NeuralNetRegres(idowhat,nxfit,nobsfit,nh,
|
|
&xfit(1:nobsfit,1:nxfit),yfit,calvalue,rsq,
|
|
&w(1:nxfit,1:nh),bph,q,bend,xnew,ysamppred(i:i))
|
|
|
|
do j=1,nobsfit
|
|
write(122,*)j,yfit(j),calvalue(j)
|
|
enddo
|
|
|
|
ipregap=i
|
|
do j=1,nxfit
|
|
ixuse_pre(j)=ixuse(j)
|
|
enddo
|
|
30 idowhat=2
|
|
call NeuralNetRegres(idowhat,nxfit,1,nh,
|
|
&xfit(1:1,1:nxfit),yfit,calvalue,rsq,
|
|
&w(1:nxfit,1:nh),bph,q,bend,xnew,ysamppred(i:i))
|
|
1000 continue
|
|
enddo
|
|
do i=1,nobs
|
|
ysamp(i)=ysamppred(i)
|
|
enddo
|
|
300 format(10f16.8)
|
|
return
|
|
end subroutine gapfilling
|