!This subroutine fills gaps in the y variable based on neural network regression
      subroutine gapfilling(nx,nobs,xsamp,ysamp,nmax)
      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.
      integer nx,nobs,nmax
      double precision xsamp(1:nobs,1:nx),ysamp(1:nobs)
!Locals
      integer i,j,k,n,idowhat,nh,nxfit,nobsfit,ixuse(nx),
     &iposdif,itakethis,iposfit(nobs),iuseit
      parameter(nh=5)
      double precision w(1:nx,1:nh),bph(nh),q(nh),bend,
     &xnew(nx),calvalue(nobs),fn9999,tiny,xfit(nobs,nx),
     &yfit(nobs),rsq,x1pre(nmax),ysamppred(nobs)
      parameter(fn9999=-9999.0d0,tiny=1.0d-6)
!
      do i=1,nmax
        x1pre(i)=fn9999
      enddo
      bend=fn9999
      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
!Fill this gap by choosing the nmax valid points that are closest to i for the fitting
        nobsfit=0
10      iposdif=10000000
        do n=1,nobs
          if(n.ne.i.and.dabs(ysamp(n)-fn9999).gt.tiny)then
            iuseit=1
!make sure it is not one that has been already selected
            do k=1,nobsfit
              if(n.eq.iposfit(k))iuseit=0
            enddo
            if(iuseit.eq.1)then
!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
            endif
            if(iuseit.eq.1)then
!make sure the distance is smaller than the current miminum
              if(iabs(n-i).lt.iposdif)then
                iposdif=iabs(n-i)
                itakethis=n
              endif
            endif
          endif
        enddo
        nobsfit=nobsfit+1
        iposfit(nobsfit)=itakethis
        yfit(nobsfit)=ysamp(itakethis)
        n=0
        do j=1,nx
          if(ixuse(j).eq.1)then
            n=n+1
            xfit(nobsfit,n)=xsamp(itakethis,j)
          endif
        enddo
        if(nobsfit.lt.nmax)goto 10
!We test to see if the same set has been used in the
!fitting before.
        do n=1,nobsfit
          if(dabs(xfit(n,1)-x1pre(n)).gt.tiny)goto 20
        enddo
!this set has been fit in the previous step
        goto 30
20      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 n=1,nobsfit
          x1pre(n)=xfit(n,1)
        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