!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
      double precision function annfunc(nparams,params,nh,q,w,bph,bend)
      implicit none
      integer nparams,nh
!nh is the number of hidden nodes in one hiden layer
!params is the inputs
!w is the weighting coefficients for the inputs
      double precision params(nparams),q(nh),
     &w(nparams,nh),bph(nh),bend
      integer i,v
      double precision term,activatefunc
      annfunc=bend
      do i=1,nh
        term=bph(i)
        do v=1,nparams
          term=term+w(v,i)*params(v)
        enddo
        annfunc=annfunc+q(i)*activatefunc(term)
      enddo
      end

      subroutine derannfunc(nparams,params,nh,q,w,bph,
     &  bend,derq,derw,derbph,derbend)
      implicit none
      integer nparams,nh
!nh is the number of hidden nodes in one hiden layer
!params is the inputs
!w is the weighting coefficients for the inputs
      double precision params(nparams),q(nh),w(nparams,nh),
     &bph(nh),bend,derq(nh),derw(nparams,nh),derbph(nh),derbend
      integer i,v
      double precision term,activatefunc,gradactivatefunc
      derbend=1.0d0
      do i=1,nh
        term=bph(i)
        do v=1,nparams
          term=term+w(v,i)*params(v)
        enddo
        derq(i)=activatefunc(term)
        derbph(i)=q(i)*gradactivatefunc(term)
        do v=1,nparams
          derw(v,i)=derbph(i)*params(v)
        enddo
      enddo
      end

      double precision function activatefunc(x)
      implicit none
      double precision x,crit
      parameter(crit=300)

!      activatefunc=2.0d0*datan(x)/3.14159265d0
!      return
 
      if(x.gt.-crit)then
        activatefunc=1.0d0/(1.0d0+dexp(-x))
      else
        activatefunc=dexp(x)/(1.0d0+dexp(x))
      endif
      return
      end
      
      double precision function gradactivatefunc(x)
      implicit none
      double precision x,crit
      parameter(crit=600)

!      gradactivatefunc=2.0d0/(3.14159265d0*(1.0d0+x*x))
!      return

      if(x.gt.-crit.and.x.lt.crit)then
        gradactivatefunc=
     &     (1.0d0/(dexp(x/2.0d0)+dexp(-x/2.0d0)))**2
      else
        gradactivatefunc=0.0d0
      endif
      return
      end

      subroutine gradannfunc(nparams,params,nh,q,w,bph,
     &der_params)
      implicit none
      integer nparams,nh
      double precision params(nparams),der_params(nparams),
     &q(nh),w(nparams,nh),bph(nh)
      integer i,v
      double precision term,dsigdterm,
     &  gradactivatefunc,activatefunc
      do i=1,nparams
        der_params(i)=0.0d0
      enddo
      do i=1,nh
        term=bph(i)
        do v=1,nparams
          term=term+w(v,i)*params(v)
        enddo
        dsigdterm=gradactivatefunc(term)
        do v=1,nparams
          der_params(v)=der_params(v)+q(i)*dsigdterm*w(v,i)
        enddo
      enddo
      return
      end