********************************************
********************************************
       PROGRAM NEUOROGEN
***    In case of Buggs -->  AZIZ@tifr.res.in
*********************************************   
       PRINT*,'                    '
       PRINT*,' WHICH EXAMPLE WANT TO TRY OUT ?'
       PRINT*,'                    '
       PRINT*,' SIMPLE PATTERNS ? GIVE 1  '
       PRINT*,' FOR  PHYSICS BASED PATTERNS GIVE 2'
       PRINT*,'                    '
       READ*,ITRY
       IF(ITRY.EQ.0) THEN
       PRINT*,'YOU HAVE NOT GIVEN ANY CHOICE, LET US TAKE 1'
       ITRY =1
       ENDIF
       IF(ITRY.EQ.1) CALL NEUROFA
       IF(ITRY.EQ.2) CALL NEUROFB  
       END
*************************************************************
        SUBROUTINE NEUSETUP(ICHOICE)
        IF(ICHOICE.EQ.1) CALL NEUSETA
        IF(ICHOICE.EQ.2) CALL NEUSETB
        END
*************************************************************
        SUBROUTINE NEUROKID(XVAR,TRGT,OUTNET,ITEST)
*************************************************************
***    Version   1.0
***    A Program started for Heavy Flavoure analysis ********
***    Written  By T. Aziz  During  1991- onwards.   ********
***    Utility routines added later.                 ********
*************************************************************
***    A Simple and efficient Program Developed for
***    Event Classification in High Energy Physics Environment
***    Can be used anywhere for Pattern Classification
***    up to 10 different Types. More can be done by
***    simple modification. 
***    Two Options possible. One based on Functional Expansion.
***    Basically two layer net. Other Multilayer Perceptron which
***    is the more common option. The Functional Expansion avoides
***    Hidden Layers and thus it is similar to Flat Net but works
***    very close to MLP for some problems.  In Functional Expansion
***     Order of Fourier like  Transform are relevant. 
*******************************************************************
***    XVAR( )Given inputs equal to NVAR numbers                  *
***    TRGT( ) Target Values intended to be achieved              *
***    OUTNET( )  Outputs given by the Network for the            *
***    corresponding output nodes as many as active nodes         * 
***    ITEST  =0 if in Training MODE else = 1 in Testing MODE     *
***    T0 is the so called temperature parameter, more relevant   *
***    beter plotting of the network output between 0 and 1.      *
***    less relevant for discrimination power of the net.         *
*******************************************************************
        COMMON/NUECOM/PI,NFR,MLPR  
        COMMON/LAYERS/ILAYER(5),LAYERS 
        COMMON/NEUINIT/MCYCLE,NVAR,NODE,LUPD,ST0,ISELF
        DIMENSION XVAR(50),TRGT(20),OUTNET(20)
C*************************************************************
        IF(MLPR.EQ.0) THEN
        CALL NEUFUN(XVAR, TRGT,OUTNET,ITEST)
        ELSEIF(MLPR.EQ.1) THEN
        CALL NEUMLP(XVAR, TRGT,OUTNET,ITEST)
        ENDIF
        RETURN
        END
************************************************************
*************************************************************
        SUBROUTINE NEUSETA
        COMMON/USRCOM/NHIDDEN(3),NODEIN,NODEOUT,MLP,NFTR
        COMMON/NEUINIT/MCYCLE,NVAR,NODE,LUPD,ST0,ISELF
C
C       SETUP THE NEWTWORK
***     SET UP THE NETWOKR
        MLP = 1                !  MULTILAYER PERCEPTRON OPTION =1
                               ! = 0 IS FUNCTIONAL EXPANSION NET
********
        IF(NVAR.EQ.0) THEN
        PRINT*,' WHERE ARE INPUTS ************ ?'
        STOP
        ENDIF
                               !  THIS ALSO FIXES THE INPUT NODES
        NODEOUT = 1            !  ONE OUTPUT FROM THE NET
        PRINT*,'                    '
        PRINT*,'GIVE THE NUMBER OF NODES IN THE FIRST HIDDEN LAYER '
        PRINT*,'                    '
        READ*,NHID
        IF(NHID.EQ.0) THEN
        PRINT*,'HOW MANY HIDDEN UNITS ?'
        ENDIF
        
        NHIDDEN(1) = NHID      !  NODES IN THE FIRST HIDDEN LAYER
        PRINT*,'                    '
        PRINT*,'GIVE THE NUMBER OF NODES IN THE NEXT HIDDEN LAYER '
        PRINT*,'NORMALLY NOT NEEDED FOR MOST OF THE PROBLEMS, 0 ? '
        PRINT*,'                    '
        READ*,NHID2
        NHIDDEN(2) = NHID2     ! AND SO ON. USUALLY ONE HIDDEN LAYER
                               ! IS ENOUGH
        CALL NEUROIN           ! SET SOME MORE DEFAULTS THAT ALSO CAN
                               ! BE CHANGED BY EXPERIENCED USER
        END
*************************************************************
        SUBROUTINE NEUSETB 
        COMMON/USRCOM/NHIDDEN(3),NODEIN,NODEOUT,MLP,NFTR
        COMMON/NEUINIT/MCYCLE,NVAR,NODE,LUPD,ST0,ISELF
        
*****    SET UP THE NETWORK
        MLP = 0                !  FUNCTIONAL EXPANSION NET
*****   or = 1 IT IS MULTILAYER PERCEPTRON
********
        IF(NVAR.EQ.0) THEN
        PRINT*,' WHERE ARE INPUTS ************ ?'
        STOP
        ENDIF
*** 
        NODEOUT = 1            !  ONE OUTPUT FROM THE NET
        PRINT*,'                    '
        PRINT*,'GIVE THE NUMBER OF EXPANSION TERMS 0 TO 6 '
        PRINT*,'GENERALLY 4 IS FOUND TO BE GOOD, BUT TRY  '
        PRINT*,'                    '
        READ*,NFR
        NFTR    = NFR          !  NUMBER OF EXPANSION TERMS
*****
        CALL NEUROIN           ! SET SOME GENERAL DEFAULTS, THAT ALSO
                               ! CAN BE CHANGED BY EXPERIENCED USER
****    THE NODE FOR WHICH SUMMARY IS REQUIRED !!!! IN CASE OF TWO CLASSES
****    THIS IS ALWAYS THE CASE.
        END
*************************************************************
       SUBROUTINE NEUROIN
*************************************************************
****  Set some defaults that can be altered by the user     *
****  in His  Calling routine bye using the common blocks.
*************************************************************
       COMMON/USRCOM/NHIDDEN(3),NODEIN,NODEOUT,MLP,NFTR
       COMMON/NEUINIT/MCYCLE,NVAR,NODE,LUPD,ST0,ISELF
       COMMON/LAYERS/ILAYER(5),LAYERS
       COMMON/NUECOM/PI,NFR,MLPR    
 
       PI = ACOS(-1.) 
***
       ISELF = 0             ! IT IS SUPERWISED NETWORK
                             ! = 1 SELF ORGANISED (TO BE DONE IN FUTURE)
       MLPR = 0              ! 1 FOR MULTILAYER PERCEPTRON 
                             ! 0 FOR FUNCTIONAL EXPANSION 
       IF(MLP.NE.0) MLPR =  MLP
       IF(MLPR.EQ.0) THEN
       NFR = NFTR 
       IF(NFTR.EQ.0) THEN
       PRINT*,'YOUR CHOICE IS SIMPLEST PERCEPTRON AND NO EXANSION'
       ENDIF
c**   Number  fourier transforms .  Usually four is reasonable
       ST0 =5.0        ! THE SIGMOID TEMPERATURE
       NODE = 1        ! THE NUMBER OF OUTPUT NODES ACTIVE
       IF(NODEOUT.NE.0) NODE = NODEOUT
C**     FOR A PROBLEM WITH ONLY TWO CLASSES 1 IS ENOUGH.
C**     FOR SEVERAL CLASSES MORE ARE NEED. A MAXIMUM OF 10 ARE ALLOWED.
c*  classification
****************************************************************** 
       ELSEIF(MLPR.EQ.1) THEN
       ST0 =2.0              ! MAY BE NEED DIFFERENT TEMPERATURE TERM
       LAYERS = 3
       ILAYER(1) = NVAR + 1     ! USER CHOICE COMES FROM INPUTS
**   KEEP SOME DEFAULTS THAT CHANGE ON CALL NEUROIN
       ILAYER(2) = 7            ! NODES IN THE FIRST HIDDE LAYER
       ILAYER(3) = 1            ! NODES IN THE OUTPUT LAYER
C       ILAYER(4) = 1           ! IF MORE LAYERS GO AHEAD
****************************************************************
       NLR = 1
       DO ILR =1, 3
       IF(NHIDDEN(ILR).NE.0) THEN
       NLR = NLR + 1
       ILAYER(NLR) = NHIDDEN(ILR)
       ENDIF
       ENDDO
       NLR = NLR +1
       IF(NODEOUT.NE.0) THEN
       ILAYER(NLR) = NODEOUT
       LAYERS = NLR
       ENDIF
       NODE      =  ILAYER(LAYERS)   ! FOR CONSISTENCY
**  THERE COULD BE OTHER LAYERS IF MORE THAN ONE HIDDEN LAYER IS NEEDED
       ENDIF
***************************************************************
       LUPD = 1         !  1 IF ETA PARAMETER TO BE UPDATED
C* IF LUPD =0  THEN ETA PARAMETER REMAINS FIXED , ESLE LUPD=1 ETA IS UPDATAED 
C* DURING TRAINING
****
       END
***********************************************************************
        SUBROUTINE NEUFUN(XVAR, TRGT,OUTNET,ITEST)
***********************************************************************
***  A Functional Expansion Network ( Flat Network) option
***  Simple and fast
***  XVAR( )Given inputs equal to NVAR numbers 
***  TRGT( ) Target Values to intended to be achieved
***  OUTNET( )  Outputs given by the Network for the corresponding
***  output nodes as many as active. Maximum 10.
***  ITEST  =0 if in Training MODE else = 1 in Testing MODE
***  T0 is the so called temperature parameter more relevant
***  beter plotting of the network output between 0 and 1.
***  less relevant for discrimination power of the net.
        COMMON/NEUINIT/MCYCLE,NVAR,NODE,LUPD,ST0,ISELF
        COMMON/NEUWT/WJI(10,1000), DWJI(10,1000), NWT, NODES
        DIMENSION FXX(1000)
        DIMENSION XVAR(50), TRGT(10), OUTNET(10)
        DIMENSION XIN(10), DELTAK(10)
********************************************
        LOGICAL FIRST
        DATA ETA/0.1/,ALFA/0.5/,T0/5.5/
        DATA KUPD/10000/
        DATA FIRST/.TRUE./    
c** calculate diff in  the calling routine where all patterns are considered.
c**    
             IF(FIRST) THEN
c** initialise
             NXI = NVAR
             IF(ST0.GT.0.5) T0 = ST0
             NODES = NODE
             IF(LUPD.EQ.1) ETA =ETA*5.

             IF(ITEST.LE.0) THEN
**   INITIALISE WTS
                   DO I=1,1000
                  DO IK=1,10
                 WJI(IK,I) = 2.*RNDM(DUM)-1.
                 DWJI(IK,I) = 0.    
                 ENDDO
              ENDDO
              ENDIF
                  IF(ITEST.NE.0) THEN
                CALL READWTS(70)              ! THERE IS NO NEED OF TRAINING
                ENDIF
             FIRST = .FALSE.
            KOUNT = 0 
          ENDIF  
C*  initialisation over
c*********************************
         KOUNT =KOUNT + 1    
**********************************
****      expand the inputs using sin and cos terms
          CALL NEUFXC(XVAR,FXX,NVAR,NKK)
             IF(ITEST.EQ.0) THEN    
            NWT = NKK 
           ELSEIF(ITEST.NE.0) THEN 

          IF( NWT.NE.NKK) THEN       

          PRINT*,' Network Not Trained for this Configuration !' 
          PRINT*,'OLD WTS , NEW WTS  = ', NWT, NKK 
          PRINT*,' THE JOB *****   STOPPED ********* '
          stop
             ENDIF
            ENDIF

               DO LN = 1, NODES
                XIN(LN) = 0.
              DO I = 1, NKK   
             XIN(LN) = XIN(LN) + WJI(LN,I)*FXX(I)
            ENDDO 
             XNET =XIN(LN)
            OUTNET(LN) = FTRANS(XNET,T0)
           ENDDO 

c** here return if not in training mode, since we do not need to do
c** any updates .

          IF( ITEST.GT.0) RETURN
          DO LN =1, NODES
       DELTAK(LN)=(TRGT(LN) - OUTNET(LN))*OUTNET(LN)*(1.-OUTNET(LN))/T0


                DO I = 1, NKK 
          DWJI(LN,I) = ETA*DELTAK(LN)*FXX(I) + ALFA*DWJI(LN,I)
          WJI(LN,I) = WJI(LN,I) + DWJI(LN,I)
           ENDDO

           ENDDO     ! ALL WTS UPDATED


              IF(LUPD.EQ.1) THEN     ! MODIFY ETA  ELSE LEAVE AS IT IS.
             IF(MOD(KOUNT,KUPD).EQ.0) THEN
             ETA = ETA*0.99   
           ENDIF            
          ENDIF

          RETURN
          END        
***************************************************************
         SUBROUTINE NEUFTX(XV,NV)       
         COMMON/NXAVG/XAV(50)
         DIMENSION XV(50)   
         DATA AX/1./     ! COULD BE ALTERED
            DO I=1, NV 
             IF(XAV(I).EQ.0.)XAV(I) =1.
            FX =(2./(1.+EXP(-AX*XV(I)/XAV(I))))-1.
            XV(I)=1.-2.*FX
           ENDDO
          RETURN
         END
         FUNCTION FTRANS(XX,T0)
           FTRANS = 1./(1. + EXP(-XX/T0 ))
C          FTRANS = 0.5 * (1. + TANH(XX/T0/2.))
          RETURN
          END
         DOUBLE PRECISION FUNCTION FTRXX(X1,X2)
         DOUBLE PRECISION X1,X2
         FTRXX = 1.D0/(1.D0+DEXP(-X1/X2))
         END  
********************************************************
      SUBROUTINE NEUDUMP(LUN)
C     DUMP WEIGHTS FOR FUTURE USE.
      COMMON/NEUWT/WJI(10,1000), DWJI(10,1000), NWT, NODES
      COMMON/NEUWTL/WJIL(4,30,30), DWJIL(4,30,30)
      COMMON/MPLNET/MLP,JLAYER(5),KLAYERS
      COMMON/NUECOM/PI,NFR,MLPR   
      open(lun,file ='dumpnet.dat',form='formatted',status ='unknown')
      IF(MLPR.LE.0) THEN
      WRITE (LUN,50) NODES, NWT
 50   format(I2,I4,10x,'Nodes and Wts in Func Exp Net') 
          do knode =1,nodes     
         do iwt=1,nwt
         write(lun,60) wji(knode, iwt), dwji(knode, iwt)
 60     format(2f15.9)
       enddo  
       enddo
         ELSE
         WRITE(LUN,70) KLAYERS,JLAYER 
 70      FORMAT(6(I2,1X), 5X,'LAYERS AND NODES IN MPL NET')
           DO I=2, KLAYERS
          DO K=1, JLAYER(I)
          DO J=1, JLAYER(I-1)
         WRITE(LUN,80) WJIL(I-1,J,K),DWJIL(I-1,J,K),(I-1),J,K
         ENDDO
        ENDDO
        ENDDO
80      FORMAT(2F15.9,3x,3I5)
      ENDIF
      close(lun)
      return
      end
***********************************************************
      SUBROUTINE READWTS(LUN)
      COMMON/NEUWT/WT(10,1000),DWT(10,1000),NWT,NODES
      COMMON/NEUWTL/WJIL(4,30,30), DWJIL(4,30,30)
      COMMON/MPLNET/MLP,JLAYER(5),KLAYERS
       COMMON/NUECOM/PI,NFR,MLPR
      print*,'mlpr = ',mlpr   
      open(lun,file ='dumpnet.dat',form='formatted',status ='unknown')
       print*,'mlpr = ', mlpr
       IF(MLPR.LE.0) THEN
      READ(LUN,50)NODES,NWT
 50   FORMAT(I2,I4)
            DO KNODE =1,NODES
          DO IWT=1, NWT
        READ(LUN , 60) WT(KNODE,IWT)
 60    FORMAT(F15.9)
       ENDDO 
      ENDDO
 
         ELSE
         READ(LUN,70) KLAYERS,JLAYER 
 70      FORMAT(6(I2,1X))
           print*,'klayer,jlayer =', klayers, jlayer
           DO I=2, KLAYERS
          DO K=1, JLAYER(I)
          DO J=1, JLAYER(I-1)
         READ(LUN,80) WJIL(I-1,J,K),DWJIL(I-1,J,K)
         ENDDO
        ENDDO
        ENDDO
80      FORMAT(2F15.9)
      ENDIF

       
      CLOSE(LUN)
      RETURN
      END 
!******************************************
      SUBROUTINE NEUFXC(XVAR,FXX,NVAR,NKK)
      COMMON/NUECOM/PI,NFR,MLPR 
      COMMON/FXKILL/KILLF(50)
      DIMENSION XVAR(50), FXX(1000)
            NKK =0      
         DO K=1,NVAR
        NKK = NKK + 1
        FXX(NKK) = XVAR(K)
          DO J = 1, NFR
         NKK = NKK + 1
           FXX(NKK) = SIN(PI*J*XVAR(K))
         NKK = NKK + 1
           FXX(NKK) = COS(PI*J*XVAR(K))
           ENDDO                 
        ENDDO               
         NKK = NKK + 1
        FXX(NKK) = 1.
       RETURN 
       END
***************************************************************
         subroutine varinit(xv,nv)
         COMMON/NXAVG/XAV(50)
         dimension xv(50)
         logical first
         data first/.true./
         if(first)then
         do i=1,nv
         xav(i)=0.
         enddo
         first=.false.
         fc=0.
         endif

         do i=1,nv
         xav(i)=(xav(i)*fc+xv(i))/(fc+1.)
         enddo
         fc=fc+1.
         return
         end

**************************************************************
       SUBROUTINE NEUSTAT(ICALL, TRGT, OUTNET, NODE)
***  Utility routine    
***  Provide Purity and Efficiency for a Given NODE
***  and Print in Histogrammable form.    
       DIMENSION XNET(200),BNET(200),ALLNET(200)
       DIMENSION PUR(200),SUMB(200),SUMALL(200)
       DIMENSION TRGT(10),OUTNET(10)
C**********    THIS ROUTINE IS USED TO GET EFFICIENCY AND PURITY
        DATA NCHAN/25/
                    
        IF(ICALL.EQ.1) THEN
C* INITIALISE
        IBALL =0
        XCHN=FLOAT(NCHAN)
        DO IC=1,NCHAN
        XNET(IC)=0.
        BNET(IC)=0.
        ALLNET(IC)=0.
        ENDDO  
        RETURN
        ENDIF

C***********  PUT IN HISTOGRAM FORM
             IF(ICALL.EQ.2) THEN
             XOUT =OUTNET(NODE)
             ICHN=INT(XOUT*XCHN) +1
             IF(ICHN.GT.NCHAN)ICHN =NCHAN  
            ALLNET(ICHN)= ALLNET(ICHN) + 1.
           IF(TRGT(NODE).GT.0.5)THEN
          IBALL = IBALL + 1
          BNET(ICHN)=BNET(ICHN)+1.
          ENDIF

********
             IF(TRGT(NODE).LT.0.5) THEN
            XNET(ICHN)=XNET(ICHN)+1.
            ENDIF
                ENDIF             ! FILL ARRAYS


*****************************************************
      IF(ICALL.LT.3) RETURN
      SB=0.
      SL=0.
      do i=1,nchan
      pur(i)=0.
      sumb(i)=0.
      sumall(i)=0.
      enddo
      DO L=1,NCHAN
      M=NCHAN-L+1
      SB=SB+BNET(M)
      SUMB(M)=SB
      SL=SL+ALLNET(M) 
      SUMALL(M)=SL
      ENDDO
      DO K=1,NCHAN
      PUR(K)=0.
      IF(SUMALL(K).GT.0) THEN
      PUR(K)=SUMB(K)/SUMALL(K)
      ENDIF
      IF(IBALL.GT.0)THEN
      SUMB(K)=SUMB(K)/FLOAT(IBALL)
      ENDIF             
      ENDDO
c****** if it is data then purity and efficiency have absolutely no meening.

****    write out the last set information for plotting if necessary
      OPEN(60,FILE ='netout.dat',form ='formatted',status ='unknown')
c      write(60,1)nchan
c 1    format(5x,'total bins = ', i10)
      PRINT*,'******  Summary For Node  = ', NODE
      PRINT*,'RESPONSE 0-1  BINNED IN CHANNELS 1-',NCHAN 
      PRINT*,'    NETA,    BNETB,    ALLNET,    PURB,    EFFB,     CHAN'
      DO I=1, NCHAN      
      IX=I
      WRITE(6,2) XNET(I), BNET(I), ALLNET(I),PUR(I),SUMB(I),IX

      write(60,2) xnet(i),bnet(i),allnet(i),pur(i),sumb(i),IX
 2    format(3f10.2,2f10.4,i5)
      enddo
      close(60)                                               

***********************************
 
      CALL NEUDUMP(70)
      RETURN
      END
*************************************************************************
**************************************************************************
        SUBROUTINE NEUMLP(XVAR, TRGT,OUTN,ITEST)
***********************************************************************
***  A MultiLayer Perceptron option 
***  XVAR( )Given inputs equal to NVAR numbers 
***  TRGT( ) Target Values to intended to be achieved
***  OUTN( )  Outputs given by the Network for the corresponding
***  output nodes as many as active
***  ITEST  =0 if in Training MODE else = 1 in Testing MODE
***  T0 is the so called temperature parameter more relevant
***  beter plotting of the network output between 0 and 1.
***  less relevant for discrimination power of the net.
        COMMON/NEUINIT/MCYCLE,NVAR,NODE,LUPD,ST0,ISELF
        COMMON/LAYERS/ILAYER(5),LAYERS
        COMMON/NEUWTL/WJIL(4,30,30), DWJIL(4,30,30)
        COMMON/MPLNET/MLP,JLAYER(5),KLAYERS
        DIMENSION XVAR(30), TRGT(30), OUTNET(5,30),OUTN(30)
C** A maximum of 30 input variables are considered. 
        DIMENSION  DELTA(5,30)

        LOGICAL FIRST
        DATA ETA/0.1/,ALFA/0.5/,T0/5.5/
        DATA KUPD/10000/
        DATA FIRST/.TRUE./    
c** calculate diff in  the calling routine where all patterns are considered.
c**    
             IF(FIRST) THEN
c** initialise
             IF(ST0.GT.0.1 ) T0 = ST0
             IF(LUPD.EQ.1) ETA =ETA*5.

             IF(ITEST.LE.0) THEN
**   INITIALISE WTS
                    DO L=1,4
                   DO I=1,30
                  DO IK=1,30
                 WJIL(L,IK,I) =2.*RNDM(DUM) -1.
***
                 DWJIL(L,IK,I) = 0.    
                 ENDDO
                ENDDO
               ENDDO
****          
               KLAYERS = LAYERS              
              DO L=1,5
              JLAYER(L)=ILAYER(L)
              ENDDO

              ENDIF
** TEST IF TRAINED NETWORK FOR THIS CONFIGURATION
             IF(ITEST.NE.0) THEN
             CALL READWTS(70)     ! GET ALL THE WTS AND TRAINED CONFIG
***      OF THE NETWORK AND THEN MAKE SURE IT AGREES WITH THE 
***      CURRENT INITIALISATIONS (LAYERS AND NODES ARE THE SAME !)
             ICHK =0
             IF(LAYERS.NE.KLAYERS) THEN
             ICHK=1
             GOTO 33
             ENDIF
             IF(LAYERS.EQ.KLAYERS) THEN
               ICHK =0
               DO IK=1, LAYERS
               IF(ILAYER(IK).NE.JLAYER(IK)) ICHK=1
              ENDDO
             ENDIF
**********    WARN
 33            IF(ICHK.NE.0) THEN
            PRINT*,'THE MLP NETWORK IS NOT TRAINED FOR '
            PRINT*,'THIS CONFIGURATION. CHECK BEFORE PROCEED '
            PRINT*,'THE PROGRAM IS STOPPED ****************** '
            STOP

           ENDIF

          ENDIF          ! FROM ITEST.NE.0 TEST
             FIRST = .FALSE.
            KOUNT = 0 
          ENDIF  
C*  initialisation over
c*********************************
         KOUNT =KOUNT + 1    
**********************************

*****   FIRST THE TARGET NODES OR THE LAST NODES
**  DEFINE ALL THE OUTNETS (FOR FIRST LAYER OUTNET =INPUTX

        DO  M =2, ILAYER(1)
        OUTNET(1,M) = XVAR(M)
        ENDDO
        OUTNET(1,1) = 1.
**  NOW GET THE OUTNET FOR OTHER LAYERS
           DO I =2, LAYERS
             DO K =1, ILAYER(I)     ! ALL NODES OF THIS LAYER
            XNET = 0.
            DO J =1, ILAYER(I-1)   ! NODES OF THE PREVIOUS LAYER
            XNET = XNET + WJIL(I-1,J,K)*OUTNET(I-1,J)
            ENDDO
             OUTNET(I,K) = FTRANS(XNET,T0)   ! USE TRANSFER FUN
             ENDDO
           ENDDO
**
***   NOW WORRY FOR UPDATES
          DO K =1, ILAYER(LAYERS)      ! THE LAST LAYER WHERE
***    TARGET IS SET FOR ERROR PROPAGATION
          OUTN(K) = OUTNET(LAYERS,K)
          DELTA(LAYERS, K)= (TRGT(K) - OUTNET(LAYERS,K))*
     & OUTNET(LAYERS,K)*(1. - OUTNET(LAYERS,K))/T0

** TAKE CARE OF THE DERIVATIVE SIGN.
       ENDDO
      IF(ITEST.EQ.1) RETURN    !  NO TRAINING REQUIRED

******     GET MORE DELTAS IF TRAINING IS TO PROCEED
        DO LM = LAYERS,2,-1
          DO I =1,ILAYER(LM-1)
         DELTAX =0.
         DO K=1, ILAYER(LM)
         DELTAX = DELTAX + DELTA(LM,K)*WJIL(LM-1,I,K)
         ENDDO
      DELTA(LM-1,I)=DELTAX*OUTNET(LM-1,I)*(1.-OUTNET(LM-1,I))/T0
 
          ENDDO
       ENDDO
**   NOW UPDATES
       DO LM=2,LAYERS
       DO I =1, ILAYER(LM-1)
       DO K =1, ILAYER(LM)
        DWJIL(LM-1,I,K)=ETA*DELTA(LM,K)*OUTNET(LM-1,I) +
     & ALFA*DWJIL(LM-1,I,K)
        WJIL(LM-1,I,K) =WJIL(LM-1,I,K) + DWJIL(LM-1,I,K)
        ENDDO
        ENDDO
       ENDDO
***   ALL THE WEGHTS UPDATED ....

*******************************************************************
*******************************************************************
              IF(LUPD.EQ.1) THEN     ! MODIFY ETA  ELSE LEAVE AS IT IS.
             IF(MOD(KOUNT,KUPD).EQ.0) THEN
             ETA = ETA*0.99        
           ENDIF            
          ENDIF

          RETURN
          END 
*****************************************************************
*****************************************************************
****  Here are some examples that one can try in order to understand
****  the performance of FEED-FORWARD-BACK-PROPAGATION NETWORKS
*****************************************************************
        SUBROUTINE NEUROFB
*****   AN EXAMPLE BASED ON L3 HCAL TYPE FAST SIMULATION FOR JETSET BASED
*****   EVENTS WHERE Z--> Q QBAR (FIVE FLAVOURES) ARE GENETAED AND ONE
*****   WOULD LIKE TAG b bbar events using only Calorimetric info.
*****   to be able to get reasonable purity neural network can not be
*****   avoided.    
        COMMON/USRCOM/NHIDDEN(3),NODEIN,NODEOUT,MLP,NFTR
        COMMON/NEUINIT/MCYCLE,NVAR,NODE,LUPD,ST0,ISELF
        common/fcom/xsbg(15,400000),ntot      ! The array where data stored
C**     THIS IS THE ARRAY CONTAINING DATA, YOU CAN HAVE YOUR OWN
C**     ARRAY OR READ FROM DATA FILES WHICHEVER CONVENIENT
        DIMENSION TRGT(10), OUTNET(10),XV(50)
        DATA BR/0.4/ , NCYCLE/150/ ,DC/0.3/ ,JCYCLE/5/
**********
        PRINT*,'                    '
        PRINT*,'GIVE NUMBER OF INPUTS '
        PRINT*,'FOR EXMAPLE 2 GOOD CHOICE IS 11, Try 10 and 12 also '
        PRINT*,'                    '
        READ*,NNVAR
        IF(NNVAR.EQ.0) THEN
        PRINT*,'HOW MANY INPUTS AND WHERE ARE THEY ?'
        STOP
        ENDIF
        NVAR   = NNVAR         !  NUMBER OF INPUT VARIABLES
********
        PRINT*,'                    '
        PRINT*,'WHICH NEURAL NET YOU WANT ? MLP OR FUNC EXPANSION '
        PRINT*,'TWO CHOICES, GIVE 1 FOR MLP, 2 FOR FUNC EXPANSION '
        PRINT*,'                    '
        READ*,ICHOICE
        IF(ICHOICE.EQ.1.OR.ICHOICE.EQ.2) THEN
        CALL NEUSETUP(ICHOICE)       ! SETUP THE NETWORK (1 FOR MLP, 2 FOR FUNCTIONAL EXPANSION)
        ELSE
        PRINT*,'YOU HAVE NOT MADE PROPER CHOICE, LET US TAKE MLP'
        ICHOICE =1
        CALL NEUSETUP(ICHOICE)       
        ENDIF
*****   
        IF(NVAR.GT.12) THEN
        PRINT*,'THE DATA GENERATED DO NOT HAVE ENOUGH VARAIBLES',NVAR
        PRINT*,' 11 IS THE RIGHT CHOICE, 12 CAN BE TRIED FOR FUN '
        PRINT*,' NO MORE EXCUTION '
        STOP
        ENDIF
        CALL QCDREAD
        KNODE = 1
****    THE NODE FOR WHICH SUMMARY IS REQUIRED !!!! IN CASE OF TWO CLASSES
****    THIS IS ALWAYS THE CASE.
c**    now train
       NTRAIN = NTOT/2
       NTEST = NTOT-NTRAIN
       PRINT*,'NTOT,NTRAIN,NTEST  =',NTOT,NTRAIN,NTEST
c***   set cycles
        DO KC =1, NCYCLE     !TRAINING CYCLES/SWEEPS
        DIFF = 0.        
        nbb=0
        nxb=0
        DO I = 1, NTRAIN                  ! START TRAINING
       IPAT = INT(RANF(DUM)*NTRAIN)+1     ! PICK UP PATTERNS RANDOMLY
         DO  J=1,NVAR
         XV(J)= XSBG(J,IPAT)
         ENDDO
         MCFL=XSBG(15,IPAT)     ! FLAVOURE CODE (1-4 for d,u,s,c, 5 b)
         IF(MCFL.LT.5) THEN
         TRGT(1) =0.01
c****   target value for the background
C***    SINCE BACKGROUND IS 4 TIMES, I TAKE ONLY PART OF IT.
C** EVEN THOUGH IT IS NOT VERY IMPORTANT BUT AVOIDES USELESS COMPUTING 
***  keep background and signal events at the same level for the purpose
***  of training.
        IF(RANF(DUM).LT.BR)THEN 
*** 
       CALL  NEUROKID(XV, TRGT,OUTNET,0)
C**                                   0  MEANS TRAINING

       DIFFX=ABS(TRGT(1) - OUTNET(1))
       DIFF= DIFF + DIFFX**2
       NXB=NXB+1
       ENDIF
        ELSE     ! LOOK FOR SIGNAL OR b quark events     
       TRGT(1) = 0.99 
c******  target value for the signal ..........

       CALL  NEUROKID(XV, TRGT,OUTNET, 0)

       DIFFX=ABS(TRGT(1) - OUTNET(1))
       DIFF= DIFF + DIFFX**2
       NBB=NBB+1
       ENDIF     
               ENDDO       ! TRAINING CYCLE OVER                    
****************************************************
       NEVT=NBB + NXB
       DIFF= DIFF/FLOAT(NEVT)
       DIFF=SQRT(DIFF)
       IF(BR.GE.0.99) THEN   
       PRINT*,'nbb, nxb,nevt,diff, cycle = ',NBB, NXB, NEVT, DIFF, KC
       ENDIF
        IF(MOD(KC,JCYCLE).EQ.0) THEN     ! GIVE THE STATUS AFTER NTH CYCLE
C**    DURING LEARNING IT IS WORTHWHILE TO WATCH THE NETWORK PERFORMANCE
c*     test where we are                         
*********
        CALL NEUSTAT(1, TRGT, OUTNET, KNODE)    !  BOOK THE STAT         
***  FOR TESTING PURITY & EFFICIENCY
************ Remeber to read wts in case one needs when
******* tarining was done before and now one is trying to test
            DIFF =0.
           DO I = NTRAIN, NTOT    
           DO J = 1, NVAR
           XV(J)=XSBG(J,I)
           ENDDO
         MCFL=XSBG(15,I)
         CALL  NEUROKID(XV, TRGT,OUTNET, 1)
          IF(MCFL.EQ.5) THEN
          TRGT(1) =0.99
          ELSEIF(MCFL.LT.5) THEN
          TRGT(1)=0.01
          ENDIF
        DIFFX=ABS(TRGT(1) - OUTNET(1))
        DIFF= DIFF + DIFFX**2     
         CALL NEUSTAT(2, TRGT, OUTNET, KNODE)   ! FILL THE STAT
         ENDDO  ! FROM NTRAIN TO NTOT 
**************
************** 
        CALL NEUSTAT(3, TRGT, OUTNET, KNODE)    !  REPORT THE OUTCOM
        DIFF =DIFF/FLOAT(NTOT-NTRAIN)
        DIFF =SQRT(DIFF)
        PRINT*,'ERROR ON TEST PATS & CYCLES ', DIFF, KC
         ENDIF          ! FROM MOD(KC,JCYCLE).
         ENDDO       !KC = CYCLES CHANGE LOOP
*** AT THE END OF THE JOB DUMP THE WEIGHTS 
       CALL NEUDUMP(70)       ! DUMP THE WTS AT THE END FOR FUTURE USE
       RETURN
       END
**************************************************************************
*****************************************************************
        SUBROUTINE NEUROFA
        COMMON/USRCOM/NHIDDEN(3),NODEIN,NODEOUT,MLP,NFTR
        COMMON/NEUINIT/MCYCLE,NVAR,NODE,LUPD,ST0,ISELF
        common/fcom/xsbg(15,400000),ntot      ! The array where data stored
C**     THIS IS THE ARRAY CONTAINING DATA, YOU CAN HAVE YOUR OWN
C**     ARRAY OR READ FROM DATA FILES WHICHEVER CONVENIENT
        DIMENSION TRGT(10), OUTNET(10),XV(50)
        DIMENSION RVEC(14)
        DATA BR/1.00/ , NCYCLE/150/ ,DC/0.3/ ,JCYCLE/5/
        DATA NTRAN/100000/,NVRN/10/,WID1/1.0/, WID2/2.0/
C**     NTRAN IS TOTAL # OF RANDOMLY GENERATED PATTERNS OF TWO CLASSES

********
        PRINT*,'                    '
        PRINT*,'GIVE NUMBER OF INPUTS '
        PRINT*,'                    ' 
        PRINT*,'FOR EXAMPLE 1 MAXIMUM IS 14,, Vary between 5-12 '
        PRINT*,'                    '
        READ*,NNVAR
        IF(NNVAR.EQ.0) THEN
        PRINT*,'HOW MANY INPUTS AND WHERE ARE THEY ?'
        STOP
        ENDIF
        NVAR   = NNVAR         !  NUMBER OF INPUT VARIABLES
*****************************************************************
        CALL NEUSETUP(1)     ! SETUP MLP NETWORK
        NVRN = NVAR
        IF(NVRN.GT.14) THEN
        PRINT*,'TO MANY INPUT VARIABLES, 14 IS MAXIMUM '
        STOP
        ENDIF
        PRINT*,'                    '
        PRINT*,' GIVE THE NUMBER OF PATTERNS TO BE GENERATED '
        PRINT*,'SHOULD BE LARGE, ABOUT 100K '
        PRINT*,'                    '
        READ*,NRAN
        IF(NRAN.NE.0)NTRAN= NRAN
*****   
        KNODE = 1
****    THE NODE FOR WHICH SUMMARY IS REQUIRED !!!! IN CASE OF TWO CLASSES
****    THIS IS ALWAYS THE CASE. THERE MAY BE OCCASION WHEN YOU NEED MORE
****    GENERATE PATTERNS
        DO IK =1, NTRAN
        CALL RNORML(RVEC,NVRN)
C*** RVEC A VECTOR OF NORMALLY DISTRIBUTED RANDOM NUMBERS ; DIMENSION NVRN
C*** ACTS AS A PATTERN VECTOR OF NVRN DIMENSIONS
C*** THE TWO CLASSES HAVE SIGMA =1. AND 2.0 ; MEAN IS ZERO FOR BOTH CLASSES
        IF(RANF(DUM).LT.0.5) THEN
        DO KK =1,NVRN
        XSBG(KK,IK)=WID1*RVEC(KK)
        ENDDO
        XSBG(15,IK)=0.01   ! TARGET VALUE
        ELSE
        DO KK =1,NVRN
        XSBG(KK,IK)=WID2*RVEC(KK)
        ENDDO
        XSBG(15,IK)=0.99   ! OTHER TARGET VALUE
        ENDIF
        ENDDO   ! ALL PATTERNS GENERATED
                
***
       NTOT = NTRAN
       PRINT*,'NTRAN, NTOT, RVEC',NTRAN,NTOT,RVEC
c***   now train
       NTRAIN = NTOT/2
       NTEST = NTOT-NTRAIN
       PRINT*,'NTOT,NTRAIN,NTEST  =',NTOT,NTRAIN,NTEST
C**    TRAIN THE NETWORK
       DO  KC=1,NCYCLE       !TRAINING CYCLES/SWEEPS/EPOCHS
       DIFF =0        
       NEVT =0
       DO  I=1,NTRAIN         ! RUN OVER TRAINING EVENTS
       IPAT = INT(RANF(DUM)*NTRAIN)+1     ! PICK UP RANDOMLY
         DO J=1,NVAR
         XV(J)=XSBG(J,IPAT)
         ENDDO
         TRGT(1) =XSBG(15,IPAT)
c****   target value for SIGNAL AND BACKGROUND
*** 
       CALL  NEUROKID(XV, TRGT,OUTNET,0)
C***                                  0 MEANS TRAIN

       diffx=abs(TRGT(1) - OUTNET(1))
       diff= diff + diffx**2
       NEVT = NEVT+1 
        ENDDO       ! TRAINING PART                   
****************************************************
       diff= diff/float(nevt)
       DIFF=SQRT(DIFF)
    
       PRINT*,'TRAIN NEVT, DIFF, CYCLE = ',NEVT,DIFF, KC
c*   go to next cycle and see if want to change any parameter etc like eta, alfa       
             IF(MOD(KC,JCYCLE).EQ.0) THEN   ! GIVE THE STATUS AFTER JCYCLES
c************************************************

********                           
c**  intialise the two routines for efficiency and muon counting
        CALL NEUSTAT(1, TRGT, OUTNET, KNODE)    !  BOOK THE STAT         
*********************
***  testing the stuff for purinty etc
************ Remeber to read wts in case one needs when
******* tarining was done before and now one is trying to test
           DIFF =0.
           DO I = NTRAIN, NTOT    
           DO J = 1, NVAR
           XV(J)=XSBG(J,I)
           ENDDO
        CALL  NEUROKID(XV, TRGT,OUTNET, 1)
C***                                    1 MEANS TEST
        TRGT(1)= XSBG(15,I)
        DIFFX=ABS(TRGT(1) - OUTNET(1))
        DIFF= DIFF + DIFFX**2     
************  get the efficiency 

          CALL NEUSTAT(2, TRGT, OUTNET, KNODE)   ! FILL THE STAT
           ENDDO  ! FROM NTRAIN TO NTOT STUFF
c******   get the relevant statistics etc......
 
        CALL NEUSTAT(3, TRGT, OUTNET, KNODE)    !  REPORT THE OUTCOM
        DIFF =DIFF/FLOAT(NTOT-NTRAIN)
        DIFF =SQRT(DIFF)
        PRINT*,'ERROR ON TEST PATS ', DIFF
        ENDIF          ! FROM MOD(KC,JCYCLE).
        ENDDO          !KC = CYCLES CHANGE LOOP
*** AT THE END OF THE JOB DUMP THE WEIGHTS 
        CALL NEUDUMP(70)       ! DUMP THE WTS AT THE END FOR FUTURE USE
        CALL NEUROFL(NVRN,WID1,WID2) ! TEST THE LIKLIHOOD FOR CLASSES
        RETURN
        END
**************************************************************************
       SUBROUTINE QCDREAD
       common/fcom/xsbg(15,400000),ntot 
       DIMENSION XV(15)
       CHARACTER*20, filename
       PRINT*,'                    '
       PRINT*,'GIVE FILENAME WHERE simulated DATA RESIDES '
       PRINT*,'Three possible ; qcdgen9.dat, qcdgen11.dat, qcdgen13.dat'
       PRINT*,'                    '
       READ*,filename
       open(30, file=filename, form='formatted',status='old')
       PRINT*,'FILE OPENED ', filename
       PRINT*,'                    '
       PRINT*,'BE PATIENT, DATA BEING READ '
       PRINT*,'                    '
       PRINT*,'=====****======'
       PRINT*,'THE 11 VARIABLES ARE : '
       PRINT*,'ENERGETIC CLUSTER IN JET 1,2 '
       PRINT*,'BOOSTS GAMMABETA  JET 1,2 '
       PRINT*,'ENERGY GAPES1,2 FOR  1&4 CLSUTER IN JET 1,2'
       PRINT*,'DIRECTED SPHERICITY JET 1,2 '
       PRINT*,'4 CLUSTER BOOST JET1,JET2   '
       PRINT*,'FRACTIONAL ENERGY OUTSIDE TWO JETS '
       PRINT*,'======******======'
       PRINT*,'                    '
       NT=0
20     READ(30,100, END=200) IP,XV
100    FORMAT(I10,7F10.6/8F10.6)

       NT = NT + 1
       DO K=1,15
       XSBG(K,NT) =  XV(K)
       ENDDO
       GOTO 20
 200   CONTINUE
       CLOSE(30)
        NTOT = NT
        PRINT*,'DATA REDING OVER, NT,NTOT =', IP,NTOT
        CALL NEUSCALE(12)
        END         
        SUBROUTINE NEUSCALE(NVAR)
        COMMON/FCOM/XSBG(15,400000),NTOT
        COMMON/NXAVG/XAV(50) 
        DIMENSION XV(15) 
        DO I=1,NTOT
        XSBG(12,I) = XSBG(15,I)
**      FOR SIMPLICITY           USE THIS FOR TRIAL TO TEST THE NET
        DO J=1,NVAR
        XV(J)=XSBG(J,I)
        ENDDO
        CALL VARINIT(XV,NVAR)
        ENDDO
********************************************* Scale factors
       DO I=1,11
       PRINT*,'VARIABLE AVERAGE  =', I, XAV(I)
       ENDDO
***    Now Scale ....
           DO I=1, NTOT
             DO J=1, NVAR
           XSBG(J,I) =XSBG(J,I)/XAV(J)
          ENDDO
         ENDDO
C***************************************************
        END  
*********************************************************
       SUBROUTINE NEUROFL(NVRN,WID1,WID2)
       common/fcom/xsbg(15,400000),ntot      ! The array where data stored
       DIMENSION RVEC(15), TRGT(10), OUTNET(10)
       CHARACTER*3 RESPONSE
C***   CHECKS THE LIKELIHOOD AND GIVES RESPONSE THAT CAN BE
C***   COMPARED WITH NEURAL NET PERFORMANCE IN THE SIMILAR FORMAT
       PRINT*,'                    '
       PRINT*,' DO YOU WANT TO TEST LIKELIHOOD RESPONSE ?'
       PRINT*,' IF YES, THEN SAY YES OR yes '
       PRINT*,'                    '
       READ*, RESPONSE
       IF(RESPONSE.EQ.'YES'.OR.RESPONSE.EQ.'yes') THEN
       NRIGHTA =0
       NRIGHTB =0
       NWRONGA =0
       NWRONGB =0
       NTEST = NTOT/2.
        CALL NEUSTAT(1, TRGT, OUTNET, 1)  ! INITIALISE
        DO IK =NTEST, NTOT
        DO JK =1, NVRN
        RVEC(JK)= XSBG(JK,IK)
        ENDDO
        TRGT(1) = XSBG(15,IK)
        PROBA = GSVEC(RVEC,WID1,NVRN)
        PROBB = GSVEC(RVEC,WID2,NVRN)
        OUT = (PROBB-PROBA)/(PROBB+PROBA)
        OUTNET(1)=((OUT + 1.)/2.)
C        PRINT*,'OUT, OUTNET ', OUT, OUTNET(1),TRGT(1)
        CALL NEUSTAT(2, TRGT, OUTNET, 1)   ! FILL
C        PRINT*,'TRGT, PROBA,PROBB =',TRGT,PROBA,PROBB
C        PRINT*,'RVEC =', RVEC
        IF(TRGT(1).GT.0.5) THEN
        IF(PROBB.GT.PROBA) THEN
        NRIGHTB= NRIGHTB + 1
        ELSE
        NWRONGB= NWRONGB + 1
        ENDIF
        ENDIF
*********
        IF(TRGT(1).LT.0.5) THEN
        IF(PROBA.GT.PROBB) THEN
        NRIGHTA = NRIGHTA + 1
        ELSE
        NWRONGA = NWRONGA + 1
        ENDIF
        ENDIF
        ENDDO   ! ALL PATTERNS ATTENDED
        PRINT*,'                    '
        PRINT*,'THIS IS LIKELYHOOD RESPONSE  '
        PRINT*,'***************************** '
        CALL NEUSTAT(3, TRGT, OUTNET, 1)   ! REPORT OUTCOME
        PRINT*,'NRIGHTA, NWRONGA = ', NRIGHTA, NWRONGA
        PRINT*,'NRIGHTB, NWRONGB = ', NRIGHTB, NWRONGB
        PRINT*,'NRIGHT, NWRONG ALL= ',NRIGHTA+NRIGHTB, NWRONGA+NWRONGB
        ENDIF     ! NO NEED TO WORK ON THIS
        END
        FUNCTION GSVEC(RVEC,SIG, NVR)
        DIMENSION RVEC(15)
        GG =1.
        DO IK =1, NVR
        XX = RVEC(IK)
        GG = GG*GSS(XX,SIG)
        ENDDO
        GSVEC=GG
        END     
        FUNCTION GSVEC2(RVEC,SIG, NVR)
        DIMENSION RVEC(15)
        GG =0.
        DO IK =1, NVR
        XX = RVEC(IK)
        GG = GG+ALOG(GSS(XX,SIG))
        ENDDO
        GSVEC2=GG
        END        
        FUNCTION GSS(XX,SIG)
***     MEAN IS ASSUMED ZERO
        DATA PI/3.141592654/
        TERM1 =SQRT(2.*PI)*SIG
        TERM2 =EXP(-XX*XX/(2.*SIG*SIG))
        GSS = TERM2/TERM1
        END
*********