!-----------------------------------------------------------------------------!
! program ga_main                                                             ! 
!-----------------------------------------------------------------------------!
!	Usage: ga_main.x [optional parameters]                                    ! 
!-----------------------------------------------------------------------------!
!	Description: General-purpose genetic algorithm inversion program.         !
!  		This code initializes a random sample of individuals with different   !
!		parameters to be optimized using the genetic algorithm approach, i.e. !
!  		evolution via survival of the fittest. The selection scheme used is   !
!  		tournament selection with a shuffling technique for choosing random   !
!  		pairs for mating. The routine includes binary coding for the          !
!  		individuals, jump mutation, creep mutation, and the option for        !
!  		single-point or uniform crossover. Niching (sharing) and an option    !
!  		for the number of children per pair of parents has been added.        !
!  		An option to use a micro-GA is also included.                         !
!-----------------------------------------------------------------------------!
!	Input parameters                                                          !
!                                                                             !
!  icreep   = 0 for no creep mutations                                        !
!           = 1 for creep mutations; creep mutations are recommended.         !
!  idum     The initial random number seed for the GA run. Must equal         !
!             a negative integer, e.g. idum=-1000.                            !
!  ielite   = 0 for no elitism (best individual not necessarily               !
!             replicated from one generation to the next).                    !
!           = 1 for elitism to be invoked (best individual replicated         !
!             into next generation); elitism is recommended.                  !
!  iend     = 0 for normal GA run (this is standard).                         !
!           = number of last population member to be looked at in a set       !
!             of individuals. Setting iend/=0 is only used for debugging      !
!             purposes and is commonly used in conjunction with iskip.        !
!  iniche   = 0 for no niching                                                !
!           = 1 for niching; niching is recommended.                          !
!  irestrt  = 0 for a new GA run, or for a single function evaluation         !
!           = 1 for a restart continuation of a GA run.                       !
!  iskip    = 0 for normal GA run (this is standard).                         !
!           = number in population to look at a specific individual or        !
!             set of individuals. Setting iskip-0 is only used for            !
!             debugging purposes.                                             !
!  itourny  No longer used. The GA is presently set up for only               !
!             tournament selection.                                           !
!  iunifrm  = 0 for single-point crossover                                    !
!           = 1 for uniform crossover; uniform crossover is recommended.      !
!  kountmx  = the maximum value of kount before a new restart file is         !
!             written; presently set to write every fifth generation.         !
!             Increasing this value will reduce I/O time requirements         !
!             and reduce wear and tear on your storage device                 !
!  maxgen   The maximum number of generations to run by the GA.               !
!             For a single function evaluation, set equal to 1.               !
!  fit_thr  = fitness threshold to assess convergence of the GA               !
!  microga  = 0 for normal conventional GA operation                          !
!           = 1 for micro-GA operation (this will automatically reset         !
!             some of the other input flags). I recommend using               !
!             npopsiz=5 when microga=1.                                       !
!  nchild   = 1 for one child per pair of parents (this is what I             !
!             typically use).                                                 !
!           = 2 for two children per pair of parents (2 is more common        !
!             in GA work).                                                    !
!  nichflg  = array of 1/0 flags for whether or not niching occurs on         !
!             a particular parameter.  Set to 0 for no niching on             !
!             a parameter, set to 1 for niching to operate on parameter.      !
!             The default value is 1, but the implementation of niching       !
!             is still controlled by the flag iniche.                         !
!  nowrite  = 0 to write detailed mutation and parameter adjustments          !
!           = 1 to not write detailed mutation and parameter adjustments      !
!  nparam   Number of parameters (groups of bits) of each individual.         !
!             Make sure that nparam matches the number of values in the       ! 
!             parmin, parmax and nposibl input arrays.                        !
!  npopsiz  The population size of a GA run (typically 100 works well).       !
!             For a single calculation, set equal to 1.                       !
!  nposibl  = array of integer number of possibilities per parameter.         !
!             For optimal code efficiency set nposibl=2**n, i.e. 2, 4,        !
!             8, 16, 32, 64, etc.                                             !
!  parmax   = array of the maximum allowed values of the parameters           !
!  parmin   = array of the minimum allowed values of the parameters           !
!  pcreep   The creep mutation probability. Typically set this                !
!             = (nchrome/nparam)/npopsiz.                                     !
!  pcross   The crossover probability. For single-point crossover, a          !
!             value of 0.6 or 0.7 is recommended. For uniform crossover,      !
!             a value of 0.5 is suggested.                                    !
!  pmutate  The jump mutation probability. Typically set = 1/npopsiz.         !
!                                                                             !
!-----------------------------------------------------------------------------!
!	Credits: For now, this program is heavily based on program ga170.f        !
!		written by David Carrol. Individual subrouitnes may be based on       !
!		other existing free code and are acknowledged in the corresponding    !
!		module file.                                                          !
!-----------------------------------------------------------------------------!
!  David Carroll's GA tips:	                                                  !
!                                                                             !
!  For single function evaluations, set npopsiz=1, maxgen=1, & irestrt=0.     !
!                                                                             !
!  My favorite initial choices of GA parameters are:                          !
!     microga=1, npopsiz=5, iunifrm=1, maxgen=200                             !
!     microga=1, npopsiz=5, iunifrm=0, maxgen=200                             !
!  I generally get good performance with both the uniform and single-         !
!  point crossover micro-GA.                                                  !
!                                                                             !
!  For those wishing to use the more conventional GA techniques,              !
!  my old favorite choice of GA parameters was:                               !
!     iunifrm=1, iniche=1, ielite=1, itourny=1, nchild=1                      !
!  For most problems I have dealt with, I get good performance using          !
!     npopsiz=100, pcross=0.5, pmutate=0.01, pcreep=0.02, maxgen=26           !
!  or                                                                         !
!     npopsiz= 50, pcross=0.5, pmutate=0.02, pcreep=0.04, maxgen=51           !
!                                                                             !
!  Any negative integer for idum should work.  I typically arbitrarily        !
!  choose idum=-10000 or -20000.                                              !
!-----------------------------------------------------------------------------!
program ga_main

	use sep
        !! genetic algorithm modules
	use ga_aux
	use ga_mutation
	use ga_crossover
	use ga_fitness
	use ga_niche
	use ga_microga
	use ga_mate
!        use formats

        !! fwd modleing function
        use sourcemod

	implicit none

	! input variables  (described above)
	integer :: kountmx
	integer :: irestrt
	integer :: itourny
	integer :: ielite
	integer :: iunifrm
	integer :: npopsiz
	integer :: nparam
	integer :: maxgen
	integer :: idum
	integer :: icreep
	integer :: nowrite 
	integer :: iniche
	integer :: nchild
	integer :: microga
	integer :: iend
	integer :: iskip
	real    :: pcross
	real    :: pmutate
	real    :: pcreep
	real	:: fit_thr
        integer :: verbose

	! internal variables
	real 	:: best		! best fitness of current generation
	real 	:: diffrac	! fraction of total number of bits between the best and
						! the rest of the micro-ga population. 
	real 	:: fbar		! average fitness of population
	integer :: ig2sum	! sum of the number of possibilities in ig2 array
	integer :: istart	! the generation to be started from
	integer :: kount	! counter to control how often the restart file is written 
	integer :: mate1	! index of population member chosen as mate1
	integer :: mate2	! index of population member chosen as mate2
	integer :: nchrome	! number of chromosomes (binary bits) of each individual
	integer :: npossum	! sum of the number of possible values of all parameters

	real,allocatable,dimension(:)	:: fitness	! 1-D array of fitness values
	real,allocatable,dimension(:)	:: g0		! lower bound of parameters
												! to be optimized
	real,allocatable,dimension(:)	:: g1		! increment of each parameter
												! from lower bounds in g0
	real,allocatable,dimension(:)	:: parmax	! maximum allowed parameter values
	real,allocatable,dimension(:)	:: parmin	! minimum allowed parameter values
	real,allocatable,dimension(:)	:: pardel	! parmax-parmin
	real,allocatable,dimension(:)	:: geni		! fitness of generation array??
	real,allocatable,dimension(:)	:: genavg	! average fitness for each gen
	real,allocatable,dimension(:)	:: genmax	! max fitness for each generat
	real,allocatable,dimension(:,:)	:: parent	! 2-D array of parents
	real,allocatable,dimension(:,:)	:: child	! 2-D array of children
	integer,allocatable,dimension(:):: nposibl	! number of possibilities per 
												! parameter. 
												! for optimal code set to 2**n
	integer,allocatable,dimension(:):: nichflg	! nich flags for each parameter 
	integer,allocatable,dimension(:):: ibest	! best indivl's binary chromosomes
	integer,allocatable,dimension(:):: ig2		! number of bits per parameter
	integer,allocatable,dimension(:,:) :: iparent ! parents' binary chromosomes
	integer,allocatable,dimension(:,:) :: ichild  ! children's binary chromosomes

	!! auxiliary variables
	integer :: outfp
	integer :: i,j,k,stat
	integer :: ncross
	integer :: ipick
	real :: evals
	real,dimension(:),allocatable :: temp

	call initpar ()

	! get actual processing parameters from the command line and set their defaults
	call from_param("irestrt", irestrt, 0)
	call from_param("microga", microga, 1)
	call from_param("npopsiz", npopsiz, 4)
	call from_param("nparam", nparam, 4)
	call from_param("pmutate", pmutate, 0.02)
	call from_param("maxgen", maxgen, 200)
	call from_param("idum", idum, -1000)
	call from_param("pcross", pcross, 0.5)
	call from_param("itourny", itourny, 1)
	call from_param("ielite", ielite, 1)
	call from_param("icreep", icreep, 1)
	call from_param("pcreep", pcreep, 0.04)
	call from_param("iunifrm", iunifrm, 1)
	call from_param("iniche", iniche, 0)
	call from_param("nchild", nchild, 1)
	call from_param("iskip", iskip, 0)
	call from_param("iend", iend, 0)
	call from_param("nowrite", nowrite, 1)
	call from_param("kountmx", kountmx, 5)
	call from_param("diffrac", diffrac, 0.05)
	call from_param("fit_thr", fit_thr, 0.01)
	call from_param("verb", verbose, 1)
	
	!! open files to write processing info out 
	outfp=24
	open(unit=outfp, file='ga.out',status='unknown',action='write',iostat=stat)
	if (stat/=0) then
		write(0,*)'Error: Could not open output file ga.out'
		stop
	end if
	open(unit=26, file='best.out',status='unknown',action='write',iostat=stat)
	if (stat/=0) then
		write(0,*)'Error: Could not open output file best.out'
		stop
	end if
	open(unit=27, file='avg.out',status='unknown',action='write',iostat=stat)
	if (stat/=0) then
		write(0,*)'Error: Could not open output file avg.out'
		stop
	end if
	open(unit=28, file='std.out',status='unknown',action='write',iostat=stat)
	if (stat/=0) then
		write(0,*)'Error: Could not open output file std.out'
		stop
	end if
	open(unit=29, file='diff.out',status='unknown',action='write',iostat=stat)
	if (stat/=0) then
		write(0,*)'Error: Could not open output file diff.out'
		stop
	end if

	!! allocate arrays for input data
	allocate(parmin(nparam),parmax(nparam),nposibl(nparam), &
		&	nichflg(nparam),temp(nparam))

	!! read auxiliary header files for min and max par values, niching
	!! flags and number of (discrete) possibilities per gene
!        call from_param("nichflg",nichflg,(/0,0,0,0,0/))
        do i=1,nparam
           nichflg(i)=0
        end do

        !!   t,x,y,z,s
        call from_param("parmin",parmin,  (/0.   ,-10000.  ,500.  ,0.0001818,-90.,-90.,-90./))
        call from_param("parmax",parmax,  (/4000.,0.       ,10000.,0.0006667,90. ,90. ,90. /))
        call from_param("nposibl",nposibl,(/50   ,102      ,10    ,100      ,90  ,90  ,90  /))
write(0,*)' parmax',parmax


!	call sreed("parmin",parmin,4*size(parmin))
!	call sreed("parmax",parmax,4*size(parmax))
!	call sreed("nichflg",temp,4*size(temp))
!	nichflg=INT(temp)
!	call sreed("nposibl",temp,4*size(temp))
!	nposibl=INT(temp)
	deallocate(temp)

!	nposibl=32768
        !! initialize travelpath (fwd modleing) operator
        call sourcemod_init(parmin,parmax,nposibl)
do i=1,nparam
	write(0,*)'i=',i,'pmin=',parmin(i),'pmax=',parmax(i),'nflg=',nichflg(i), &
		& 'npos=',nposibl(i)
end do

	!! allocate auxiliary arrays
	allocate(g0(nparam),g1(nparam),ig2(nparam),pardel(nparam))

	!! call subroutine to initialize and check GA parameters
	call initialize_parameters(nparam,idum,nchrome,microga,npossum,ig2sum, &
		&	parmin,parmax,pardel,nposibl,g0,g1,ig2,outfp)

	!! allocate population-related arrays
	allocate(parent(nparam,npopsiz),child(nparam,npopsiz),fitness(npopsiz),&
		&	geni(maxgen),genavg(maxgen),genmax(maxgen),iparent(nchrome,npopsiz), &
		&	ichild(nchrome,npopsiz),ibest(nchrome))

	!! perform necessary initialization and read the ga.restart file
	idum=1
        call initialize_population(istart,npossum,ig2sum,nparam, &
		&   nowrite,microga,idum,irestrt,npopsiz,nchrome,outfp, &
		&   parent,iparent,nposibl,g0,g1,ig2,parmax)


	!!!!!!!!!!!main loop: over generations !!!!!!!!!!!!!!!!!!!!!!
	kount=0
	do i=istart,maxgen+istart-1
		if(verbose.gt.5) write(0,1111) i
		write(24,1111) i
		write(24,1050)

		!! evaluate population, assign fitness, establish the best individual
		!! and write output information
		call evalout(maxgen,npopsiz,nchrome,nparam,iskip,iend,ibest,fbar,best, &
			& fit_thr,fitness,parent,iparent,g0,g1,ig2,verbose)

		!! check for dummy errors
		if(npopsiz == 1 .or. iskip /= 0) then
			write(0,*) 'ERROR: npopsize==1 or iskip /=0'
			close(24)
			stop
		endif

		!! save best and average fitness for this generation
		geni(i)=float(i)
		genavg(i)=fbar
		genmax(i)=best

		!! Implement "niching".
		if (iniche /= 0) then
			call niche (nparam,npopsiz,nchrome, &
        		& parent,iparent,fitness,nposibl,nichflg,parmax,parmin,pardel,verbose)
		end if

		!! Enter selection, crossover and mutation loop.
		ipick=npopsiz
		do j=1,npopsiz,nchild

			!! Randomly select a pair of parents to mate
			call selectn(ipick,j,itourny,nparam,npopsiz,nchrome,nchild, &
				& mate1,mate2,parent,iparent,child,ichild,fitness)

			!! Now perform crossover between the randomly selected pair.
			call crossover(ncross,j,mate1,mate2,nparam,npopsiz,nchrome, &
				& nchild,iunifrm,pcross,iparent,ichild)

		end do
		if(verbose.gt.5) write(0,1225) ncross
		write(24,1225) ncross

		!! Do random jump mutations. If running micro-GA, skip this
		!! (only creep mutations are allowed with micro-ga)
		!if (microga == 0) then
			call jump_mutate(npopsiz,nchrome,nowrite,pmutate,ichild,verbose)
		!end if

		!! Do creep mutations
		if (icreep /=0) then
			call creep_mutate(nparam,npopsiz,nchrome,nowrite,pcreep, &
   	    		& ichild,child,nposibl,g0,g1,ig2,parmax,parmin,verbose)
		end if
		
		!! Write child array back into parent array for new generation. Check
		!! to see if the best parent was replicated.
		call newgen(ielite,npossum,ig2sum,npopsiz,nparam,nchrome,nowrite, &
			& ibest,nposibl,parent,child,iparent,ichild,g0,g1,ig2,parmax)

		!! If using micro-ga, check for need to restart the population
		if (microga /= 0) then
			call gamicro(nparam,npopsiz,nchrome,nowrite,i,npossum,ig2sum, &
				& ibest,diffrac,parent,iparent,g0,g1,ig2,nposibl,parmax,verbose)
		end if

		!! Write to restart file. This prevents loosing all work if something
		!! goes wrong.
		call restart(i,istart,kount,kountmx,maxgen,npopsiz,nchrome,iparent)

	end do !!!!!! End of main generational processing loop
        call sep_write(picksout)

	!! output results
	write(24,3000)
	do i=1,maxgen
		evals=float(npopsiz)*geni(i)
		write(24,3100) geni(i),evals,genavg(i),genmax(i)
	end do

	!! close output files
	close (24); close (25); close (26); close (27)

	!! deallocate working space
	deallocate(parmin,parmax,nposibl,nichflg,g0,g1,ig2,pardel,ichild, &
		&	ibest,iparent,parent,child,fitness,geni,genavg,genmax)

	!! include file with formats for read and write statements
!	include 'formats'
        !! formats for ga_main
1050    format (1x, '  -      Binary Code', 16x, 'Param1  Param2  Fitness')
1111	format('1','-----------------  Generation',i 5,'  -----------------')
1225	format(/'  Number of Crossovers      =',i5)
1400	format(2x,'CPU time for all generations=',e12.6,' sec'/ &
			&	2x,'                             ',e12.6,' min')
3000	format(2x/'Summary of Output'/ &
			&	2x,'Generation   Evaluations   Avg.Fitness   Best Fitness')
3100	format(2x,3(e10.4,4x),e11.5)

		!! formats for ga_microgas
1375 	format(/'%%%%%%%  Restart micro-population at generation', &
			&	i5,'  %%%%%%%')

end program ga_main
