c     Unbiased determination of the proton structure function
c     F2p with faithful uncertainty estimation
c     The NNPDF Collaboration: Luigi Del Debbio, Stefano Forte,
c     Jose Ignacio Latorre, Andrea Piccione, Joan Rojo 
c     
c     This code is organized as follows:
c     
c     - f2av(x,q2,it,f2,er) needs x, q2 and the kind of target
c       (1 proton, -1 deuteron, 0 proton minus deuteron) as inputs
c       and gives f2 and its error as an output. The result is given
c       as an average on 1000 neural networks.
c     - f2cor(xa,q2a,ita,xb,q2b,itb,f2a,era,f2b,erb,cor) needs, 
c       q2 and the kind of target for two points as inputs and 
c       gives f2, their errors and their correlation as an output.
c       The result is given as an average on 1000 neural networks.
c     - net(x,q2,it,inet,f2,itimep,info) needs x, q2 and the kind 
c       of target as inputs and gives f2 and its error as an output. 
c       The result is the output of the inet-th neural network.
c       The flag itimep prevents reading twice the neural network
c       parameters if it is not necessary.
c       Info gives information on the used kinematical range where
c       the fit may fail:
c            * info=0 x and q2 are inside the range
c            * info=1 the value of x is outside the range
c            * info=2 the value of q2 is outside the range
c            * info=3 the values of x and q2 are outside the range
c

      subroutine f2av(xa,q2a,ita,f2a,era)
      implicit double precision(a-h,o-z)
      parameter(mxp=650,mxsets=5001)
      dimension f2pdnet(mxsets)
      itimep=0
      nets=1000

      do kk=1,nets
         call net(xa,q2a,ita,kk,f2,itimep,info)
         f2pdnet(kk)=f2
      enddo

c Average on 1000 neural networks of central value and error 

      x=0d0
      x2=0d0
      do kk=1,nets
         x=x+f2pdnet(kk)
         x2=x2+f2pdnet(kk)**2
      enddo
      f2a=x/nets
      f22a=x2/nets
      era=sqrt(f22a-f2a**2)
      return
      end

      subroutine f2cor(xa,q2a,ita,xb,q2b,itb,f2a,era,f2b,erb,cor)
      implicit double precision(a-h,o-z)
      parameter(mxp=650,mxsets=5001)
      dimension f2pdneta(mxsets),f2pdnetb(mxsets)

      itimep=0
      nets=1000

      do kk=1,nets
         call net(xa,q2a,ita,kk,f21,itimep,info)
         f2pdneta(kk)=f21
      enddo

      itimep=0
      if(ita.eq.itb)itimep=1

      do kk=1,nets
         call net(xb,q2b,itb,kk,f22,itimep,info)
         f2pdnetb(kk)=f22
      enddo

c Average on 1000 neural networks of central values, 
c errors and correlations

      x=0d0
      x2=0d0
      y=0d0
      y2=0d0
      z=0d0
      do kk=1,nets
         x=x+f2pdneta(kk)
         x2=x2+f2pdneta(kk)**2
         y=y+f2pdnetb(kk)
         y2=y2+f2pdnetb(kk)**2
         z=z+f2pdneta(kk)*f2pdnetb(kk)
      enddo
      f2a=x/nets
      f22a=x2/nets
      era=sqrt(f22a-f2a**2)
      f2b=y/nets
      f22b=y2/nets
      erb=sqrt(f22b-f2b**2)
      cor=(z/nets-f2a*f2b)/era/erb
      return
      end

      subroutine net(x,q2,it,inet,f2,itimep,info)
      implicit real*8(a-h,o-z)
      integer tnl, nets,itimep,inet
      parameter(mxn=20,mxp=650,mxl=5)
      parameter(mxnets=5001,nets=1000)
      dimension wp(mxn,mxn,mxl,mxnets),thetap(mxn,mxl,mxnets)
      dimension zeta(mxn,mxl)
      dimension xmin(mxnets),q2min(mxnets)
      dimension xlogmin(mxnets),q2logmin(mxnets)
      dimension xspread(mxnets),q2spread(mxnets)
      dimension xlogspread(mxnets),q2logspread(mxnets)
      dimension outmin(mxnets),spreadout(mxnets)
      dimension n(mxl)
      save tnl,n,xspread,q2spread,xmin,q2min,wp,thetap
      save xlogmin,q2logmin,xlogspread,q2logspread,outmin,spreadout

c Reading of weights and thresholds

      info=0
      if(x.gt.0.75.or.x.lt.0.000006)info=1
      if(q2.gt.30000.or.q2.lt.0.045)info=2
      if(x.gt.0.75.or.x.lt.0.000006.and.
     .    q2.gt.30000.or.q2.lt.0.045)info=3

      if(itimep.eq.0) then
         if(it.gt.0)open(unit=77,type='old',
     .                            name="f2p.weights")
         if(it.lt.0)open(unit=77,type='old',
     .                            name="f2d.weights")
         if(it.eq.0)open(unit=77,type='old',
     .                            name="f2ns.weights")
         do  ixarxes=1,nets
    	     read(77,*) tnl
	     read(77,*) (n(kk),kk=1,tnl)	
             read(77,*) xmin(ixarxes),q2min(ixarxes),
     .                  xlogmin(ixarxes),q2logmin(ixarxes)              
	     read(77,*) xspread(ixarxes),q2spread(ixarxes),
     .                  xlogspread(ixarxes),q2logspread(ixarxes)         
	     read(77,*) outmin(ixarxes),spreadout(ixarxes)
             do  l=2,tnl
              do  j=1,n(l-1)
                do  i=1,n(l)
                  read(77,*) wp(i,j,l,ixarxes)
                enddo 
              enddo
             enddo
             do  l=2,tnl
              do  i=1,n(l) 
                  read(77,*)thetap(i,l,ixarxes)
              enddo
             enddo
         enddo
         itimep=1
         close(77)
      endif 
      f2=0

c Normalization of input parameters

         xnor=(x-xmin(inet))/xspread(inet)
         xnor=xnor*0.8+0.1
         q2nor=(q2-q2min(inet))/q2spread(inet)
         q2nor=q2nor*0.8+0.1

         xlog=log(x)
         xlognor=(xlog-xlogmin(inet))/xlogspread(inet)
         xlognor=xlognor*0.8+0.1

         q2log=log(q2)
         q2lognor=(q2log-q2logmin(inet))/q2logspread(inet)
         q2lognor=q2lognor*0.8+0.1

c Evaluation of the neural network output

        zeta(1,1)=xnor
        zeta(2,1)=q2nor
        zeta(3,1)=xlognor
        zeta(4,1)=q2lognor
        do   l=2,tnl-1
         do  i=1,n(l)
             h=0.
             do j=1,n(l-1)
               h=h+wp(i,j,l,inet)*zeta(j,l-1)
             enddo
             h=h-thetap(i,l,inet)
             zeta(i,l)=1./(1.+DExp(-h))
         enddo
        enddo 
        do i=1,n(tnl)
             h=0.
             do j=1,n(l-1)
               h=h+wp(i,j,l,inet)*zeta(j,l-1)
             enddo
             h=h-thetap(i,l,inet)
             zeta(i,tnl)=h
        enddo       
        f2=(zeta(1,tnl)-0.1)/0.8
        f2=spreadout(inet)*f2+outmin(inet)
      return
      end