subroutine NeuralNetRegres(idowhat,nx0,nobs0,nh0,xsamp0,
     &ysamp0,yatxsamp0,rsq,w,bph,q,bend,xnew,ypred)
      implicit none
      include 'NeuralNetRegres.h'
!
!=============Inputs regardless of idowhat=========================
!idowhat: =1, fit the data and estimate the coefficients. Provide the
!             initial guess for the coefficients or set bend to -9999
!         =2, coefficients are already available, calculate y at xnew
!nx0: the number of independent (x) variables
!nobs0: the total number of samples
!nh0: the total number of hidden nodes to use. One hidden layer is
!     assumed.
!============When idowhat=1========================================
!       --------Inputs--------
!xsamp0: the values of the independent (x) variables
!ysamp0: the values of the dependent (y) variable. y is one dimension.
!       --------Outputs-------
!w: the slope coefficient to time the normalized x in the activation function
!bph: the intercept coefficient in the activation function
!q: the coefficient to time the value of the activation function
!bend: the residual constant in the neural network regression
!yatxsamp0: the predicted y value at xsamp0
!rsq: R squared
!============When idowhat=2=========================================
!       --------Inputs--------
!w: the slope coefficient to time the normalized x in the activation function
!bph: the intercept coefficient in the activation function
!q: the coefficient to time the value of the activation function
!bend: the residual constant in the neural network regression
!xnew: the new x point who y value is to be estimated (when idowhat=2)
!       --------Outputs-------
!ypred: the predicted y value at xnew
!
      integer idowhat,nx0,nobs0,nh0
      double precision xsamp0(nobs0,nx0),ysamp0(nobs0),
     &  yatxsamp0(nobs0),rsq,w(nx0,nh0),bph(nh0),q(nh0),
     &  bend,xnew(nx0),ypred
!============Locals=========================================
      integer i,j,ndim,ny,INFO,iderivative,iregrestype
!iregrestype=0, ordinary distance regression
      double precision xnormk(nx0),xnormb(nx0),std,fmean,
     &  xmin,xmax,fatbeta,fatbeta0,fatbetacp,ftol,
     &  beta(nx0*nh0+2*nh0+1),betacp(nx0*nh0+2*nh0+1),rms,
     &  agrind,ran2,annfunc,weitx(1:nobs0,1:nx0),
     &  weity(1:nobs0),shortx(1:nobs0,1:nx0),
     &  shorty(1:nobs0),yv(nobs0),fn9999,tiny
      parameter(ftol=1.0d-8,iderivative=1,iregrestype=0,
     &fn9999=-9999.0d0,tiny=1.0d-8)
      external funkmin_neural,f1dim_neural,FCN_neural  
!
      if(idowhat.eq.1)then
!Regression
        nx=nx0
        nobs=nobs0
        nh=nh0
!xnormk: the slope of the linear transformation for xsamp0
!xnormb: the intercept of the linear transformation for xsamp0
!Transform xsamp to become bounded (-1,1) so that different independent variables
!are comparable in magnitude
!xmin ~ -1
!xmax ~ +1
        do i=1,nobs
          ysamp(i)=ysamp0(i)
          weity(i)=1.0d0
          do j=1,nx
            weitx(i,j)=1.0d0
          enddo
        enddo
        do i=1,nx
          call stdmaxmeanmin(nobs,xsamp0(1:nobs,i:i),
     &                     std,fmean,xmin,xmax)
          if(xmax.eq.xmin)then
            xnormk(i)=1.0d0
            xnormb(i)=0.0d0
          else
            xnormk(i)=2.0d0/(xmax-xmin)
            xnormb(i)=-(xmax+xmin)/(xmax-xmin)
          endif
          do j=1,nobs
            xsamp(j,i)=xnormk(i)*xsamp0(j,i)+xnormb(i)
          enddo
        enddo
        ndim=2*nh+nh*nx+1
        do i=1,ndim
          betamin(i)=-1.0d+20
          betamax(i)=1.0d+20
        enddo
        if(dabs(bend-fn9999).lt.tiny)then
!no initial guess. Use the general guess
          do i=1,ndim
            beta(i)=(ran2()-0.5d0)*2.0d0
          enddo
        else
!initial guess provided. transform the guessed bph and w coefficients to correspond
!to the transformed x.
          do i=1,nh0
            do j=1,nx0
              w(j,i)=w(j,i)/xnormk(j)
            enddo
          enddo
          do i=1,nh0
            do j=1,nx0
              bph(i)=bph(i)-w(j,i)*xnormb(j)
            enddo
          enddo
          call coeff_beta(2,nx,nh,beta,w(1:nx,1:nh),bph,q,bend)
        endif
        do i=1,ndim
          betacp(i)=beta(i)
          if(beta(i).lt.betamin(i).or.beta(i).gt.betamax(i))then
!            write(*,*)'Inproper initial guess in NeuralNetRegres.f'
            beta(i)=ran2()
          endif
          betacp(i)=beta(i)
        enddo
        call funkmin_neural(ndim,beta,fatbeta0)
        INFO=iregrestype
        ny=1
        fatbeta=fatbeta0
90      call odr_leastsquare(ndim,FCN_neural,beta,nobs,
     &xsamp(1:nobs,1:nx),nx,ysamp(1:nobs),ny,weitx(1:nobs,1:nx),
     &weity(1:nobs),iderivative,shortx(1:nobs,1:nx),shorty(1:nobs),
     &fatbeta,INFO)
        call funkmin_neural(ndim,beta,fatbeta)
!        if((fatbeta0-fatbeta).gt.ftol)then
!          fatbeta0=fatbeta
!          do i=1,ndim
!            betacp(i)=beta(i)
!          enddo
!          goto 90
!        endif
        if(fatbeta.gt.fatbeta0)then
          j=0
          do i=1,ndim
            beta(i)=betacp(i)
            if(beta(i).lt.betamin(i).or.beta(i).gt.betamax(i))j=1
          enddo
          fatbeta=fatbeta0
          if(j.ne.0)then
            do i=1,ndim
              beta(i)=(ran2()-0.5d0)*2.0d0
            enddo
            call funkmin_neural(ndim,beta,fatbeta)
          endif
        endif
100     fatbetacp=fatbeta
        call nongradopt(ndim,funkmin_neural,f1dim_neural,
     &beta,betamin,betamax,ftol,fatbeta)
!        call RepeatCompassSearch(ndim,beta,fatbeta,
!     &betamin,betamax,funkmin_neural,f1dim_neural,ftol)
!        if(dabs(fatbetacp-fatbeta).gt.ftol)goto 100
        call funkmin_neural(ndim,beta,fatbeta)
        call coeff_beta(idowhat,nx,nh,beta,w(1:nx,1:nh),bph,q,bend)
!transform the estimated bph and w coefficients so that the original x 
!values can be used directly.
        do i=1,nh0
          do j=1,nx0
            bph(i)=bph(i)+w(j,i)*xnormb(j)
          enddo
        enddo
        do i=1,nh0
          do j=1,nx0
            w(j,i)=w(j,i)*xnormk(j)
          enddo
        enddo
        do i=1,nobs0
          yatxsamp0(i)=annfunc(nx0,xsamp0(i:i,1:nx0),nh0,q,
     &      w(1:nx0,1:nh0),bph,bend)
        enddo
        call rsq_rms(ysamp0,yatxsamp0,nobs0,rsq,rms,agrind)
      endif
      if(idowhat.eq.2)then
!Predict y at x with the regression coefficients already estimated
        ypred=annfunc(nx0,xnew,nh0,q,w(1:nx0,1:nh0),bph,bend)
      endif
      return
      end
!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&