!$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
      double precision function 
     &annfunc(nparams,params,nh,q,w,bph,bend,c)
      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,c(nh)
      integer i,v
      double precision term,activatefunc1,activatefunc2
      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)*activatefunc1(term)+
     &c(i)*activatefunc2(term)
      enddo
      end

      subroutine derannfunc(nparams,params,nh,q,w,bph,
     &  bend,c,derq,derw,derbph,derbend,derc)
      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,c(nh),derq(nh),
     &  derw(nparams,nh),derbph(nh),derbend,derc(nh)
      integer i,v
      double precision term,activatefunc1,gradactivatefunc1,
     &activatefunc2,gradactivatefunc2

      derbend=1.0d0
      do i=1,nh
        term=bph(i)
        do v=1,nparams
          term=term+w(v,i)*params(v)
        enddo
        derq(i)=activatefunc1(term)
        derc(i)=activatefunc2(term)
        derbph(i)=q(i)*gradactivatefunc1(term)+
     &            c(i)*gradactivatefunc2(term)
        do v=1,nparams
          derw(v,i)=derbph(i)*params(v)
        enddo
      enddo
      end

      subroutine gradannfunc(nparams,params,nh,q,w,bph,
     &  c,der_params)
      implicit none
      integer nparams,nh
      double precision params(nparams),der_params(nparams),
     &  q(nh),w(nparams,nh),bph(nh),c(nh)
      integer i,v
      double precision term,gradactivatefunc1,gradactivatefunc2
      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
        do v=1,nparams
          der_params(v)=der_params(v)+(q(i)*gradactivatefunc1(term)
     &+c(i)*gradactivatefunc2(term))*w(v,i)
        enddo
      enddo
      return
      end

      double precision function activatefunc1(x)
      implicit none
      double precision x,crit
      parameter(crit=300)
      if(x.gt.-crit)then
        activatefunc1=1.0d0/(1.0d0+dexp(-x))
      else
        activatefunc1=dexp(x)/(1.0d0+dexp(x))
      endif
      return
      end

      double precision function gradactivatefunc1(x)
      implicit none
      double precision x,crit
      parameter(crit=600)
      if(x.gt.-crit.and.x.lt.crit)then
        gradactivatefunc1=
     &     (1.0d0/(dexp(x/2.0d0)+dexp(-x/2.0d0)))**2
      else
        gradactivatefunc1=0.0d0
      endif
      return
      end
 
      double precision function activatefunc2(x)
      implicit none
      double precision x
      
!      activatefunc2=2.0d0*datan(x)/3.14159265d0
!      activatefunc2=1.001d0+dsin(x)
      activatefunc2=x+x*x
      return
      end

      double precision function gradactivatefunc2(x)
      implicit none
      double precision x,crit
      parameter(crit=600)
      
!      gradactivatefunc2=2.0d0/(3.14159265d0*(1.0d0+x*x))
!      gradactivatefunc2=dcos(x)
      gradactivatefunc2=1.0d0+2.0d0*x
      return
      end