! -------------------------------
! a simple 3 layers (1 hidden layer) neural network
! (c) 2011 Bui Hoang Hai
! -----------------------------------
  subroutine controller_init 
    Integer :: i,j, Mvars
    Real :: tmp
        print*,trim(input_file)
        open(1,file=trim(input_file),status='old')
        Read (1,*)
        Read (1,*) Nsample, Mvars ! Nin
        !Read (1,*); Read (1,*)
        Nin = Mvars-1; Nout=1
        print*,"==>", Nin, Nout, Nsample

        call ann_init   ! Init ann first

        allocate(Inputs(Nin,NSample));
        allocate(Outputs(Nout,NSample))
        allocate(min_in(Nin));
        allocate(max_in(Nin))
        allocate(min_out(Nout));
        allocate(max_out(Nout))

        if (training_mode) then ! training_mode = .T.
           do j=1,NSample
              read(1,*) (Outputs(i,j),i=1,Nout), (Inputs(i,j),i=1,Nin)
           enddo
           do i=1,Nin
              min_in(i)=MINVAL(Inputs(i,:))
              max_in(i)=MAXVAL(Inputs(i,:))
           enddo
           do i=1,Nout
              min_out(i)=MINVAL(Outputs(i,:))
              max_out(i)=MAXVAL(Outputs(i,:))
           enddo

        else                    ! training_mode = .F.

           do j=1,NSample
              read(1,*) tmp, (Inputs(i,j),i=1,Nin)
           enddo

           call controller_load_info
         !   call controller_print_info
        endif
        close(1)
  end subroutine controller_init

!---------------------------------------------    
  subroutine controller_train_data
     integer :: epoc, i, j

     do epoc=1,Nepoc
        mean_error=0.
        do j=1, Nsample
           do i=1,Nin
              Fin(i) = rescale(Inputs(i,j),min_in(i),max_in(i),-1.,1.)
           enddo
           Call ann_forward

           do i=1,Nout
              Ftarget(i) = rescale(Outputs(i,j),min_out(i),max_out(i),0.2,0.8)
           enddo
           Call ann_backward
           mean_error=mean_error+error/Nsample
        enddo
     enddo
     ! call controller_print_info
     call controller_save_info
  end subroutine controller_train_data

!---------------------------------------------    
  subroutine controller_output_data 
    Integer :: i, j
     open(3,file=output_file,status='unknown')

     write(3,"(3I8)") Nin, Nout, Nsample

     do j=1, Nsample
        do i=1,Nin
           Fin(i) = rescale(Inputs(i,j),min_in(i),max_in(i),-1.,1.)
        enddo
        Call ann_forward
        do i=1,Nout
           Ftarget(i) = rescale(Fout(i),0.2,0.8,min_out(i),max_out(i))
        enddo
        Yhq(j) = Ftarget(1)   ! Out TT
        write(3,"(30F12.4)") (Inputs(i,j),i=1,Nin),(Ftarget(i),i=1,Nout)
     enddo
     close (3)
  end subroutine controller_output_data

!---------------------------------------------    
  subroutine controller_print_info
    Integer :: i, j
        write(*,*)"Input,Hidden,Output"
        write(*,*) Nin,Nhidden,Nout
        write(*,*)'Weight_input-hidden:'
        do i=0,Nin
           write(*,*) (weight_ih(i,j),j=1,Nhidden)
        enddo
        write(*,*)'Weight_hidden-output:'
        do i=0,Nhidden
           write(*,*) (weight_ho(i,j),j=1,Nout)
        enddo
    write(*,*)"min_in:"
    write(*,*) min_in
    write(*,*)"Max_in:"
    write(*,*) max_in
    write(*,*)"min_out:"
    write(*,*) min_out
    write(*,*)"Max_out:"
    write(*,*) max_out
  end subroutine controller_print_info
!---------------------------------------------
  ! Save weights and scaling factors
  subroutine controller_save_info
    Integer :: i, j
     open(1,file=save_file,status='unknown')
     write(1,*) "Input,Hidden,Output:"
     write(1,"(3I8)") Nin,Nhidden,Nout
     write(1,*) "Weight_input-hidden:"
     do i=0,Nin
        write(1,"(30F12.7)") (weight_ih(i,j),j=1,Nhidden)
     enddo
     write(1,*) "Weight_hidden-output:"
     do i=0,Nhidden
        write(1,"(12F12.7)") (weight_ho(i,j),j=1,Nout)
     enddo
     write(1,*) "min_in:"
     write(1,"(30F12.7)") min_in
     write(1,*) "Max_in:"
     write(1,"(30F12.7)") max_in
     write(1,*) "min_out:"
     write(1,"(30F12.7)") min_out
     write(1,*) "Max_out:"
     write(1,"(30F12.7)") max_out
     close (1)
     ! print*, "Terminate controller_save_info"
  end subroutine controller_save_info

!---------------------------------------------
  subroutine controller_load_info
    integer :: Nin1, Nout1, i, j
    open(1,file=save_file,status='old')
        read(1,*)
        read(1,*) Nin1,Nhidden,Nout1
        if ((Nin.ne.Nin1).or.(Nout.ne.Nout1)) then
           write(*,*) "Number of inputs,outputs in save file are not equal to those in input file!"
           stop
        endif

        read(1,*)
        do i=0,Nin
           read(1,*) (weight_ih(i,j),j=1,Nhidden)
        enddo
        read(1,*)
        do i=0,Nhidden
           read(1,*) (weight_ho(i,j),j=1,Nout)
        enddo
    read(1,*)
    read(1,*) min_in
    read(1,*)
    read(1,*) max_in
    read(1,*)
    read(1,*) min_out
    read(1,*)
    read(1,*) max_out
  end subroutine controller_load_info

!---------------------------------------------
 subroutine controller_deinit
        if (Allocated(Inputs)) deallocate(Inputs);
        if (Allocated(Outputs)) deallocate(Outputs)
        if (Allocated(min_in)) deallocate(min_in);
        if (Allocated(max_in)) deallocate(max_in)
        if (Allocated(min_out)) deallocate(min_out);
        if (Allocated(max_out)) deallocate(max_out)
        if (Allocated(weight_ih)) deallocate(weight_ih)
        if (Allocated(weight_ho)) deallocate(weight_ho)
        if (Allocated(fin)) deallocate(fin)
        if (Allocated(fhid)) deallocate(fhid)
        if (Allocated(fout)) deallocate(fout)
        if (Allocated(ftarget)) deallocate(ftarget)
        if (Allocated(delta_o)) deallocate(delta_o)
        if (Allocated(delta_h)) deallocate(delta_h)
        if (Allocated(moment_ih)) deallocate(moment_ih)
        if (Allocated(moment_ho)) deallocate(moment_ho)
        if (debug) close(111)
 end subroutine controller_deinit

!---------------------------------------------
!  calc y from x such that [x1,x2]-->[y1,y2]
 real function rescale(x,x1,x2,y1,y2)
    real x,x1,x2,y1,y2
    rescale = y1 + (y2-y1)*(x-x1)/(x2-x1)
 end function rescale
!---------------------------------------------
!---------------------------------------------
    subroutine ann_init
      Integer :: iin, ihid
      real tmp

        allocate(weight_ih(0:Nin,Nhidden))
        allocate(weight_ho(0:Nhidden,Nout))
        allocate(moment_ih(0:Nin,Nhidden))
        allocate(moment_ho(0:Nhidden,Nout))

        ! init weight that takes values randomly from [-1,1]
        do iin=0,Nin
           do ihid=1,Nhidden
              call RANDOM_NUMBER(tmp)
              weight_ih(iin,ihid)= 2. - tmp
           enddo
        enddo
        do ihid=0,Nhidden
           do iout=1,Nout
              call RANDOM_NUMBER(tmp)
              weight_ho(ihid,iout)= 2. - tmp
           enddo
        enddo
        moment_ho = 0.
        moment_ih = 0.

        allocate(fin(0:Nin))
        allocate(fhid(0:Nhidden))
        allocate(fout(Nout))
        allocate(ftarget(Nout))
        allocate(delta_o(Nout))
        allocate(delta_h(Nhidden))
        !---bias
        fin(0)=1.
        fhid(0)=1.
        if (debug) open(111,file='debug.txt')
    end subroutine ann_init

!---------------------------------------------
    subroutine ann_forward
      real :: asum
      Integer :: ihid, iout, iin
        ! Tinh dau ra cua lop an
        do ihid=1,Nhidden
           asum=0.
           do iin=0,Nin
              asum = asum + fin(iin)*weight_ih(iin,ihid)
           enddo
           fhid(ihid) = sigma(asum)
        enddo

        ! Tinh dau ra cua lop dau ra
        do iout=1,Nout
           asum =0.
           do ihid=0,Nhidden
              asum = asum + fhid(ihid)*weight_ho(ihid,iout)
           enddo
           fout(iout) = sigma(asum)
        enddo
        if (debug) then
           write(111,*)'FORWARD'
           write(111,*)'..fin:',fin
           write(111,*)'..weight_ih:'
           do j=1,Nhidden
              write(111,*)(weight_ih(i,j),i=0,Nin)
           enddo
           write(111,*)'..fhid',fhid
           write(111,*)'..weight_ho',weight_ho
           write(111,*)'..fout',fout
        endif

    end subroutine ann_forward

!---------------------------------------------
      !--- NOTE that you have already run forward_prop
      !--- Require ftarget
    subroutine ann_compute_error
       error=0.
       do iout=1,Nout
              delta_o(iout)=(ftarget(iout)-fout(iout))
              error=error+(delta_o(iout)**2)/Nout
       enddo
    end subroutine ann_compute_error

!---------------------------------------------
      !---NOTE that you have already run forward_prop
      !---Update weights for output layer
    subroutine ann_backward
      real :: delta_ih,delta_ho
      Integer :: iout, ihid, iin
     !---Calculate error, delta and update weights for the output layers
    error=0.
    do iout=1,Nout
       delta_o(iout)=(ftarget(iout)-fout(iout))
       error=error+(delta_o(iout)**2)/Nout         !--Mean square error
       delta_o(iout)=delta_o(iout)*fout(iout)*(1-fout(iout))  !--Derivative of Sigma:  y' = y(1-y)  
       do ihid=0,Nhidden
          delta_ho=learning_rate*delta_o(iout)*fhid(ihid)
          weight_ho(ihid,iout) = weight_ho(ihid,iout) + delta_ho + moment_ho(ihid,iout)
          moment_ho(ihid,iout)=momentum*delta_ho            !--momentum for next itteration
       enddo
    enddo
     !---Update weights for hidden layer
    do ihid=1,Nhidden
       delta_h(ihid)=0.
       do iout=1,Nout
          delta_h(ihid) = delta_h(ihid) + weight_ho(ihid,iout)*delta_o(iout)
       enddo
       delta_h(ihid) = delta_h(ihid)*fhid(ihid)*( 1- fhid(ihid))
       do iin=0,Nin
          delta_ih = learning_rate*delta_h(ihid)*fin(iin)
          weight_ih(iin,ihid)=weight_ih(iin,ihid) + delta_ih + moment_ih(iin,ihid)
      !--momentum for next itteration
          moment_ih(iin,ihid)=momentum*delta_ih
       enddo
    enddo

    if (debug) then
       write(111,*)'weight_ih',weight_ih
    endif
     !---
    if (debug) then
       write(111,*)'BACK'
       write(111,*)'--fout',fout
       write(111,*)'--ftarget',ftarget
       write(111,*)'delta_o',delta_o
       write(111,*)'weight_ih',weight_ih
       write(111,*)'weight_ho',weight_ho
    endif
  end subroutine ann_backward
!---------------------------------------------    
      !---Sigmoid function
    real function sigma(x)
    real x
        sigma=1./ (1. + exp(-x))
    end function sigma
!---------------------------------------------