!*************************************************************************
!*                                                                       *
!*     Initialization of parameters for LMBNNR                           *
!*     (version 1.10, last modified 11.11.2021)                          *
!*                                                                       *
!*************************************************************************
!*
!*     Modules included:
!*
!*     initnn           ! Initialization of parameters for NNR.
!*     initlmbm         ! Initialization of LMBM -solver.
!*
 
MODULE initnn  ! Initialization of parameters for NN.

    USE r_precision, ONLY : prec  ! Precision for reals.
    IMPLICIT NONE

    ! Input and output files (change the name of the output file here): 
    CHARACTER(LEN=80), SAVE :: &
        infile, &                    ! Name of dataset file.
        outfile0 = 'result1.txt', &  ! Result file with indices, function values, and cpu-time.
        outfile1 = 'result2.txt'     ! Result file with predictions.

    ! Input real parameters.
    REAL, SAVE :: &
        tlimit    = 86400.0          ! Maximum CPU-time in seconds, from user.
    

    ! Input integer parameters
    INTEGER, SAVE :: &
        nft, &                       ! Number of input features in data, from user.
        nrecord, &                   ! Number of recorts in data, from user.
        ntrain, &                    ! Size of the training set, from user.
        nH_max, &                    ! Maximum number of nodes, from user, by default nH_max = 10.
        noutcom, &                   ! Output column, from user, by default noutcom = nft+1.
        istop = 1, &                 ! Automated stopping procedure is
                                     !   0 - not used,
                                     !   1 - used (recommended).
        isetx = 4, &                 ! Updating method for weights (see subroutine set_x below). 
                                     !   isetx = 4 by default.
        iscale = 1, &                ! Scaling is
                                     !   0 - not used,
                                     !   1 - used (recommended).  
        ipredi = 1                   ! Prediction:
                                     !   0 - only at the termination with complite model.
                                     !   1 - every step.
    
    
    ! Allocatable real arrays.
    REAL(KIND=prec), SAVE, DIMENSION(:,:), allocatable :: &
        a, &                         ! Data set, from input file, a(nrecord,nft+1).
        mX, &                        ! Input matrix, mX(nft,nrecord).
        ascaled                      ! Scaled data matrix.
    REAL(KIND=prec), SAVE, DIMENSION(:), allocatable :: & !
        mY, &                        ! Output vector, mY(nrecord).                             
        pred, &                      ! the values of the predicted outputs
        clin, &                      ! For scaling.
        dlin                         ! For scaling.
    
    ! Other Real parameters.
    REAL(KIND=prec), SAVE :: & !
        rho = 1.0, &                 ! Penalty parameter.
        asp = 0.00001, &             ! Stopping tolerance.
        fplain                       ! f without penalization.
        
    ! Other integer parameters.
    INTEGER, SAVE :: & !
        nH, &                        ! Current number of clusters, loops from 1 to nH_max
        ns, &                        ! Switch for auxiliary and real NN problem. ! NOT YET USED.
        maxdim                       ! Maximum number of variables in optimization, maxdim = nH_max*(nft+2)+1.
        

CONTAINS

    SUBROUTINE init_x(x0,n0)          ! Initial x. May be changed by user. Nevertheless, the initial x needs to be given here.
            
        USE param, ONLY : zero, one               ! Parameters.
        IMPLICIT NONE
        
        INTEGER, INTENT(IN) :: n0                 ! Initial number of variables.
        REAL(KIND=prec), INTENT(OUT) :: x0(n0)    ! Initial x.
    
        
        x0 = one/REAL(nft,prec) 
        x0(1) = zero                  ! bias from hidden layer to output.
        x0(n0) = one                  ! from hidden layer to output.
        x0(n0-1) = zero               ! bias from input layer to hidden layer.
                    
    END SUBROUTINE init_x
    
    SUBROUTINE set_x(x0,n0)           ! Updating x.
            
        USE param, ONLY : zero, one               ! Parameters.
        IMPLICIT NONE
        
        INTEGER, INTENT(IN) :: n0                 ! Number of variables.
        REAL(KIND=prec), INTENT(INOUT) :: x0(n0)  ! Next x: on input first n-(nft+2) values are the previous solution.
                                      ! About variables:    w_1BP   <--->  x(1)
                                      !                     w_0kj   <--->  x(1+(j-1)*(nft+2)+k)
                                      !                     w_0Bj   <--->  x(1+(j-1)*(nft+2)+nft+1)    
                                      !                     w_1jP   <--->  x(1+j*(nft+2))    (i.e.  x(1+(j-1)*(nft+2)+nft+2) )
                                      ! Note n0 = nH*(nft+2)+1
    
        INTEGER :: i
        
        SELECT CASE (isetx)
            CASE(1)
                ! previous result as obtained, new values are set 0 for bias, 1 from hidden to output and 1/nft to all others.
                DO i = n0-(nft+2)+1, n0
                    x0(i)=one/REAL(nft,prec) ! adds the desired value to "new" part of x.
                END DO
                 
                x0(n0-1) = zero ! new bias
                x0(n0) = one    ! from new hidden node to output
        
            CASE(2)
                ! no incremental, at every iteration the starting point is 0 for bias, 1 from hidden to output and 1/nft to all others. 
                DO i = 1, n0 
                    x0(i)=one/REAL(nft,prec) ! all values are 1/nft
                END DO
        
                x0(1) = zero
                DO i=1,nH  
                    x0(i*(nft+2)) = zero !bias terms
                    x0(i*(nft+2)+1) = one !hidden to output terms
                END DO
                x0(n0-1) = zero ! new bias
                x0(n0) = one    ! from new hidden node to output
                
        
            CASE(3)            
                ! previous result as obtained, new values are set to 0 but from new node to output that is 1.
                DO i = n0-(nft+2)+1, n0-1 ! all new values are 0.
                    x0(i) = zero
                END DO
                x0(n0) = one ! from new node to output.
        
            CASE(4)
                DO i = n0-(nft+2)+1, n0
                    x0(i-(nft+2)) = x0(i-(nft+2))/2.0_prec ! half values from input to hidden layer, changes "old" part of x (one weight).
                    x0(i) = x0(i-(nft+2)) ! adds the desired value to "new" part of x.
                END DO  
            
                x0(n0-(nft+2)) = 2.0_prec*x0(n0-(nft+2)) ! whole value from hidden to output
                x0(n0)=x0(n0-(nft+2))                    ! whole value from hidden to output
                
            CASE(5)
                DO i = n0-(nft+2)+1, n0
                    x0(i) = x0(i-(nft+2)) ! adds the desired value to the "new" part of x.
                END DO  
            
                x0(n0-(nft+2)) = x0(n0-(nft+2))/2.0_prec ! 1/2 from hidden to output
                x0(n0)=x0(n0-(nft+2))
            
            CASE(6)
                DO i = n0-(nft+2)+1, n0
                    x0(i) = -x0(i-(nft+2)) ! adds the desired value to "new" part of x.
                    x0(i-(nft+2)) = x0(i-(nft+2))*2.0_prec ! doubles values from input to hidden layer, changes "old" part of x (one weight).
                END DO  
                x0(n0-(nft+2)) = x0(n0-(nft+2))/2.0_prec ! whole value from hidden to output
                x0(n0)=x0(n0-(nft+2))
        END SELECT    
            
    END SUBROUTINE set_x
    

    SUBROUTINE init_nnpar()         ! Subroutine for initialization of parameters.
                               
        IMPLICIT NONE
        
        maxdim = nH_max*(nft+2)+1

    END SUBROUTINE init_nnpar

    !===================================================================
    !  Scaling
    !===================================================================

    SUBROUTINE scaling

        USE param, ONLY : zero, one               ! Parameters.
        IMPLICIT NONE

        ! Local variables
        REAL(KIND=prec) ::      &
            dm,var
        INTEGER ::  &
            i,j

        ! Intrinsic Functions
        INTRINSIC ABS,SQRT

        DO i=1,nft+1
            dm=zero
            DO j=1,nrecord   
                dm=dm+ascaled(j,i)
            END DO
            dm=dm/REAL(nrecord,prec)
            var=zero
            DO j=1,nrecord
                var=var+(ascaled(j,i)-dm)**2
            END DO
            var=var/REAL(nrecord-1,prec)
            var=SQRT(var)

            IF (var >= 1.0E-05_prec) THEN
                clin(i)=one/var
                dlin(i)=-dm/var
                DO j=1,nrecord
                    ascaled(j,i)=clin(i)*ascaled(j,i)+dlin(i)
                END DO
            ELSE
                IF (ABS(dm) <= one) THEN
                    clin(i)=one
                    dlin(i)=zero
                ELSE
                    clin(i)=one/dm
                    dlin(i)=zero
                END IF
                DO j=1,nrecord
                    ascaled(j,i)=clin(i)*ascaled(j,i)+dlin(i)
                END DO
            END IF
        END DO

        RETURN

    END SUBROUTINE scaling

END MODULE initnn


MODULE initlmbm  ! Initialization of parameters for LMBM. There is, usually, no need to change these parameters.

    USE r_precision, ONLY : prec   ! Precision for reals.
    USE param, ONLY : zero, one    ! Parameters.
    
    IMPLICIT NONE

    ! Parameters
    INTEGER, PARAMETER :: &
        na      = 2, &             ! Maximum bundle dimension, na >= 2.
        mcu     = 15, &             ! Maximum number of stored corrections, mcu >=1.
        mcinit  = 7, &             ! Initial maximum number of stored corrections, mcu >= mcinit >= 3.
                                   ! If mcinit <= 0, the default value mcinit = 3 will be used.
                                   ! However, the value mcinit = 7 is recommented.
        inma    = 3, &             ! Selection of line search method:
                                   !   inma = 0, Armijo line search,
                                   !   inma = 1, nonmonotone Armijo line search.
                                   !   inma = 2, weak Wolfe line search.
                                   !   inma = 3, nonmonotone  weak Wolfe line search.
        mnma    = 10, &            ! Maximum number of function values used in nonmonotone line search.
        maxnin  = 20               ! Maximum number of interpolations, maxnin >= 0.
                                   ! The value maxnin = 2-20 is recommented with inma=0,
                                   ! maxnin >= 20 with inma=1 and 3, and maxnin =200 with inma=2.
                                   ! For example:
                                   !   inma = 0, maxin = 20.
                                   !   inma = 1, mnma = 20, maxin = 30.
                                   !   inma = 2, maxnin = 200.
                                   !   inma = 3, mnma=10, maxnin = 20.


    ! Real parameters (if parameter value <= 0.0 the default value of the parameter will be used).
    REAL(KIND=prec), SAVE :: &
        tolb    = zero, &         ! Tolerance for the function value (default = -large).
        tolf    = zero, &         ! Tolerance for change of function values (default = 1.0E-8).
        tolf2   = -10.0_prec, &   ! Second tolerance for change of function values.
                                  !   - If tolf2 < 0 the the parameter and the corresponding termination
                                  !   criterion will be ignored (recommended with inma=1,3).
                                  !   - If tolf2 = 0 the default value 1.0E+4 will be used.
        tolg    = 1.0E-5_prec, &  ! Tolerance for the termination criterion (default = 1.0E-5).
        tolg2   = 1.0E-3_prec, &  ! Tolerance for the second termination criterion (default = 1.0E-3).
        eta     = 1.0E-4_prec, &  ! Distance measure parameter, eta > 0.
                                  !   - If eta < 0  the default value 0.0001 will be used.
        epsl    = 0.24E+00, &     ! Line search parameter, 0 < epsl < 0.25 (default = 0.24).
        xmax    = 1000.0_prec     ! Maximum stepsize, 1 < XMAX (default = 1000).

    ! Integer parameters (if value <= 0 the default value of the parameter will be used).
    INTEGER, SAVE :: &
        mit     = 300, &          ! Maximun number of iterations.
        mfe     = 300, &          ! Maximun number of function evaluations.
        mtesf   =      0, &       ! Maximum number of iterations with changes of
                                  ! function values smaller than tolf (default = 10).
        iprint  =      0, &       ! Printout specification:
                                  !    -1  - No printout.
                                  !     0  - Only the error messages.
                                  !     1  - The final values of the objective function
                                  !          (default used if iprint < -1).
                                  !     2  - The final values of the objective function and the
                                  !          most serious warning messages.
                                  !     3  - The whole final solution.
                                  !     4  - At each iteration values of the objective function.
                                  !     5  - At each iteration the whole solution
        iscale  =     0           ! Selection of the scaling with LMBM:
                                  !     0  - Scaling at every iteration with STU/UTU (default).
                                  !     1  - Scaling at every iteration with STS/STU.
                                  !     2  - Interval scaling with STU/UTU.
                                  !     3  - Interval scaling with STS/STU.
                                  !     4  - Preliminary scaling with STU/UTU.
                                  !     5  - Preliminary scaling with STS/STU.
                                  !     6  - No scaling.

    REAL(KIND=prec), SAVE, DIMENSION(:), allocatable :: x  ! Vector of variables 
    REAL(KIND=prec), SAVE :: forig                 ! f at x0.
    INTEGER, SAVE :: n                             ! Number of variables.
    

CONTAINS

    SUBROUTINE defaults()  ! Default values for parameters.

        USE param, ONLY: small, large, zero, one, half
        USE initnn, ONLY : nft, nrecord, nH
        
        IMPLICIT NONE

        IF (iprint < -1) iprint  = 1               ! Printout specification.
        IF (mit   <= 0) mit      = 5000            ! Maximum number of iterations.
        IF (mfe   <= 0) mfe      = n*mit           ! Maximum number of function evaluations.
        IF (tolf  <= zero) tolf  = 1.0E-08_prec    ! Tolerance for change of function values.
        IF (tolf2 == zero) tolf2 = 1.0E+04_prec    ! Second tolerance for change of function values.
        IF (tolb  == zero) tolb  = -large + small  ! Tolerance for the function value.
        IF (tolg  <= zero) tolg  = 1.0E-05_prec    ! Tolerance for the termination criterion.
        IF (tolg2  <= zero) tolg  = 1.0E-03_prec   ! Tolerance for the second termination criterion.
        IF (xmax  <= zero) xmax  = 1.5_prec        ! Maximum stepsize.
        IF (eta   <  zero) eta   = half            ! Distance measure parameter
        IF (epsl  <= zero) epsl  = 1.0E-04_prec    ! Line search parameter,
        IF (mtesf <= 0) mtesf    = 10              ! Maximum number of iterations with changes
                                                   ! of function values smaller than tolf.
        IF (iscale > 6 .OR. iscale < 0) iscale = 0 ! Selection of the scaling.

        IF (nft*nrecord < 10000) THEN !
            IF (nH <= 70) THEN
                mit = 300
            ELSE IF (nH <= 100) THEN
                mit = 100
            ELSE
                mit = 50    
            END IF
        ELSE IF (nft*nrecord < 1000000) THEN !
            IF (nH <= 20) THEN
                mit = 300
            ELSE IF (nH <= 50) THEN
                mit = 50
            ELSE IF (nH <= 100) THEN
                mit = 20
            ELSE
                mit = 10    
            END IF
        ELSE IF (nft*nrecord < 100000000) THEN !
            IF (nH <= 10) THEN
                mit = 50
            ELSE IF (nH <= 50) THEN
                mit = 10
            ELSE  
                mit = 5    
            END IF
        ELSE !
            IF (nH <= 5) THEN
                mit = 50
            ELSE IF (nH <= 10) THEN
                mit = 20
            ELSE IF (nH <= 50) THEN
                mit = 10 
            ELSE 
                mit = 5   
            END IF
        END IF

        mfe=mit
        PRINT*,'hihu'

    END SUBROUTINE defaults

    SUBROUTINE init_lmbmpar()  ! User supplied subroutine for further initialization of parameters
                               ! (when needed) for LMBM. May be left empty.

        IMPLICIT NONE

    END SUBROUTINE init_lmbmpar

END MODULE initlmbm
