!! module to implement the functions related to the fitness function evaluation.
!! The actual cost function is problem dependent and is implemented in a
!! separate module. Include everyting related to your function evaluation into
!! a module and replace the statement "use ga_Mauro_function, only: func" with 
!! the name of your module. This module will comunicate with your module only
!! through function func which should have the following interface:
!!					call func(npopsiz,nparam,j,funcval,parent)
!! where all parameters are input, except funcval which is the fitness value 
!! for the individual "j" in the "parent" population
module ga_fitnessPOINT
  use ga_aux
  use sourcemod
  implicit none
contains
  !------------------------------------------------------------------------
  subroutine evalout(maxgen,npopsiz,nchrome,nparam,iskip,iend,ibest,fbar,best,&
  & fit_thr,fitness,parent,iparent,g0,g1,ig2,verb,flat)
    !------------------------------------------------------------------------
    ! This subroutine evaluates the population, assigns fitness,
    ! and establishes the best individual.
    !------------------------------------------------------------------------

    !! subroutine arguments
    integer, intent(in) :: maxgen,verb
    integer, intent(in) :: npopsiz
    integer, intent(in) :: iskip
    integer, intent(in) :: iend
    integer, intent(in) :: nchrome
    integer, intent(in) :: nparam
    real, intent(in)	:: fit_thr
    real, intent(out)   :: fbar
    real, dimension(nparam,npopsiz), intent(inout)		:: parent
    integer, dimension(nchrome,npopsiz),intent(in)	:: iparent
    integer, dimension(nchrome), intent(out)		:: ibest
    real, dimension(npopsiz), intent(out)			:: fitness
    real, dimension(nparam),intent(in)				:: g0
    real, dimension(nparam),intent(in)				:: g1
    integer,dimension(nparam),intent(in)			:: ig2
    real,dimension(nparam+1)::tmp
    
    !! local variables
    integer		:: j,n,k,kk,stat
    integer		:: jbest ! member of the population with the best fitness
    integer,save :: count=1
    real		:: best,prev_best
    real		:: fitsum		!sum of fitness of parents
    real		:: fitsumsq		!sum of square fitness of parents
    real		:: std_dev		!standard deviation for a given generation
    logical:: flat
    save prev_best

    ! initialization
    fitsum=0.0d0
    fitsumsq=0.0d0
    best=100000.0d10

    !! loop over individuals in the population
    do j=1,npopsiz
       !! decode the current chromosome
       call decode(j,nparam,npopsiz,nchrome,parent,iparent,g0,g1,ig2)

       ! evaluate fitness of current individual
       call func(npopsiz,nparam,j,fitness(j),parent)

       ! update fitness and parameter sums
       fitsum=fitsum+fitness(j)
       fitsumsq=fitsumsq+fitness(j)**2

       ! check to see if fitness of this individual is the best fitness
       if (fitness(j) < best) then
          best=fitness(j)
          jbest=j
          ibest(1:nchrome)=iparent(1:nchrome,j)
       end if

    end do

    ! compute fitness average and standard deviation for this generation
    fbar=fitsum/dble(npopsiz)
    std_dev=sqrt((dble(npopsiz)*fitsumsq-fitsum**2)/(dble(npopsiz)* &
    & dble(npopsiz-1)))

    ! write output information
    if(verb.gt.5) write(0,*)'best=',best,' prev_best=',prev_best
    if(flat.gt.5) write(0,*)'fbar=',fbar,'std_dev=',std_dev
    if(flat) write(26,'(f19.6)') best
    if(flat) write(27,'(f19.6)') fbar
    if(flat) write(28,'(f19.6)') std_dev
    if(flat) write(29,'(f19.6)') fbar-best
    if(flat) write(30,'(4(f8.1,1x),2(i4,1x),f9.4)')parent(1:3,jbest),1./parent(4,jbest),parent(5,jbest)*180/3.141593,(parent(6,jbest)*180/3.141593),best
    if(exist_file("model")) then
       tmp(1)=parent(1,jbest)
       tmp(2)=parent(2,jbest)
       tmp(3)=parent(3,jbest)
       tmp(4)=1./parent(4,jbest)
       tmp(5)=nint(parent(5,jbest)*180/3.141593)
       tmp(6)=nint(parent(6,jbest)*180/3.141593)
       tmp(7)=best
       call srite("model",tmp,4*size(tmp))
    end if

    !! check termination criterion
    if ((best<fit_thr .and.(abs(best-fbar)<1.).and. &
    & (abs(prev_best-best)<0.00001)).or.count==maxgen) then
!    if( (abs(prev_best-best) .le. 0.00001) .or. count.eq.maxgen) then

       !! solution found.
       !! call the fitness function once more just to write the trace out
       call func(npopsiz,nparam,jbest,fitness(jbest),parent,flag=1)
       if(count/=maxgen) then
          call to_history("n2",count,"model")
          call sep_close()
          write(0,*) 'TERMINATION CRITERIUM SATISFIED. SOLUTION FOUND'
       else
          write(0,*) 'TERMINATION DUE TO GENERATION LIMIT maxgen:',maxgen
          write(0,*) '     POSSIBLE PISS-ANT SOLUTION'
          call to_history("n2",count,"model")
          call sep_close()
       end if
       write(0,*) 'total number of function evaluations =', &
       & (count+1)*npopsiz
       write(0,*) 'fitness value:',best,' average fitness value:',fbar
       open(unit=31,file='solution',status='unknown',action='write', &
       & iostat=stat)
       if (stat/=0) then
          write(0,*) 'ERROR: could not open file: solution' 
          stop
       else
          write(31,*)'x,y,z,s,v ' 
          do j=1,nparam
             write(31,*) parent(j,jbest)
          end do
          write(31,*) 1/parent(4,jbest)
       end if
       close(31)
       call reportdone(parent(:,jbest))
    end if
    !! update variables for next generation.
    prev_best=best
    count=count+1
    !		include 'formats'
    !! formats for ga_fitness
1075 format(i3,1x,1485i1,99(1x,f9.4),1x,f9.5)
    !1075	format(i3,1x,30i1,2(1x,f7.4),1x,f8.5)
1100 format(1x,'Average Function Value of Generation=',f8.5)
1200 format(1x,'Maximum Function Value              =',f8.5/)
1275 format(/' Average Values:',18x,99(1x,f9.4),1x,f8.5/)
    !1275	format(/' Average Values:',18x,2(1x,f7.4),1x,f8.5/)

    return

  end subroutine evalout
  !----------------------------------------------------------------------

end module ga_fitnessPOINT
