      function funmin(p,n)
      implicit double precision (a-h,o-z)
      save
      include 'NET_PARAMS'
      include 'POINT_PARAMS'
      dimension p(n),gan(nt)
      common/tograd/xt(nop),yt(nop)
      data itric/0/
      ONCE: if (itric.eq.0) then
              itric = 1
              if (n.ne.nt) then
                 write(*,*) 'N must equal to: ',nt
                 stop'Rerun'
              endif
              x1 = -4.d0      ! Prepare the training set
              xn =  4.d0
              do i = 1,nop
                xi = x1 + float(i-1)/float(nop-1)*(xn-x1)
                yi = xi*sin(xi**2)
                xt(i) = xi
                yt(i) = yi
              enddo
            endif ONCE
* Regularization (Weight decay, quadratic)
            wpen = 0.d0
            do i=1,nt
              wpen=wpen + p(i)**2
            enddo
            call getfla(2,str)    ! Regularization Penalty
            wpen = wpen*str/nt
            funmin = 0.
            do i = 1,nop
              call twohlnn(xt(i),p,an,gan)
              funmin = funmin + (an-yt(i))**2
            enddo
            funmin = funmin/nop + wpen
* PRINTING
            call getfla(1,pri)
            if (pri.eq.1.d0) then
                call setfla(1,0.d0)
                write(*,100)
                write(*,101)
                write(*,102) M1,M2,LD,NT
                write(*,100)
100   format(40('='))
101   format(t6,'M1',t11,'M2',t16,'LD',t21,'NT')
102   format(t5,i3,t10,i3,t15,i3,t20,i3)
                iu = nunit()
                open(iu,file='RES')
                do i = 1,nop
                  call twohlnn(xt(i),p,an,gan)
                  write(iu,*) xt(i),yt(i),an,an-yt(i)
                enddo
                close(iu)
            endif
* END OF PRINTING
            end
      subroutine granal(n,p,grad)
      implicit double precision (a-h,o-z)
      save
      include 'NET_PARAMS'
      include 'POINT_PARAMS'
      dimension p(n),gan(nt),grad(n)
      common/tograd/xt(nop),yt(nop)
* Regularization (Weight decay, quadratic)
            call getfla(2,str)    ! Regularization Penalty
            do ip =1,n
              grad(ip) = 0.d0
              do i = 1,nop
                call twohlnn(xt(i),p,an,gan)
                grad(ip) =  grad(ip) + 2*(an-yt(i))*gan(ip)/nop
              enddo
              grad(ip) = grad(ip) + 2.d0*str*p(ip)/nt
            enddo
      end

      subroutine twohlnn(x,p,an,gan)
! Implementation of a Neural Network with two hidden layers.
! INPUT:  x (the pattern), p (the network's weights).             
! OUTPUT: AN (the ANN), ANX (its x-gradient), ANXX (its x-Hessian)
      implicit double precision (a-h,o-z)
      include 'NET_PARAMS'
      dimension p(NT),gan(NT)
      dimension a(M2),w(M2,M1),c(M2),v(M1,LD),b(M1),x(LD)
      dimension A1(M1),A2(M2)
      dimension A1F(M1),A2F(M2),A1D(M1),A2D(M2),A1SD(M1),A2SD(M2)
      do j=1,M2                  ! Associate P(.) to a,w,c
        a(j) = p(j)
        c(j) = p(M2+M1*M2+j)
          do k = 1,M1
            w(j,k) = p(M2+M1*(j-1)+k) 
          enddo
      enddo
*      
      do k= 1,M1                 ! Associate P(.) to v,b
        b(k) = p(2*M2+M1*M2+M1*LD+k)
          do l = 1,LD
             v(k,l) = p(2*M2+M1*M2+LD*(k-1)+l)
          enddo
      enddo
* Now do all the work using a,w,c,v,b       
      do k = 1,M1                ! Store recurrent expressions A1,A1D
        A1(k) = 0.d0
        do l = 1,LD
          A1(k) = A1(k) + v(k,l)*x(l)
        enddo
        A1(k) = A1(k) + b(k)
        A1F(k) = f0(A1(k))
        A1D(k)  = f1(A1(k))
        A1SD(k) = f2(A1(k))
      enddo 
      do j = 1,M2                ! Store recurrent expressions A2,A2D
        A2(j) = 0.d0
        do k = 1,M1
          A2(j) = A2(j) + w(j,k)*A1F(k)
        enddo
        A2(j) = A2(j) + c(j)
        A2F(j) = f0(A2(j))
        A2D(j)  = f1(A2(j))
        A2SD(j) = f2(A2(j))
      enddo
      an   =  0.d0              ! Calculate the ANN 
      do j=1,M2
          an = an + a(j)*A2F(j)
      enddo   
* Calculate GAN(.), i.e. the ANN gradient wrt the weights p(.)
      do j=1,M2
        gan(j) = A2F(j)                    ! wrt a_j
        gan(M2+M1*M2+j) = a(j)*A2D(j)      ! wrt c_j
         do k = 1,M1
           ineq = M2+M1*(j-1)+k
           gan(ineq) = a(j)*A2D(j)*A1F(k)  ! wrt w_{jk}
         enddo
      enddo
      do k=1,M1
        AF1W = 0.d0
        do j=1,M2
          AF1W = AF1W + a(j)*w(j,k)*A2D(j)
        enddo
        do l = 1,LD
          ineq = 2*M2+M1*M2+LD*(k-1)+l
          gan(ineq) = A1D(k)*x(l)*AF1W
        enddo
        gan(2*M2+M1*M2+M1*LD+k) = A1D(k)*AF1W
      enddo
*
      end

      function f0(x)
      implicit double precision (a-h,o-z)
      f0 = 1.d0/(1.d0+exp(-x))
      end
      function f1(x)
      implicit double precision (a-h,o-z)
      z =  1.d0/(1.d0+exp(-x))
      f1 = z*(1.-z)
      end
      function f2(x)
      implicit double precision (a-h,o-z)
      z = 1.d0/(1.d0+exp(-x))
      f2 = z*(1.-z)*(1.-2*z)
      end
      function f3(x)
      implicit double precision (a-h,o-z)
      z = 1.d0/(1.d0+exp(-x))
      zp = z*(1-z)
      f3 = zp*(1.-6*zp)
      end
