! module for auxiliary subroutines to the genetic algorithm optimization
! code gaf90. Subroutines based on David Carrol's ga170

module ga_aux

!	use formats
	implicit none
	contains
	!--------------------------------------------------------------------
	subroutine initialize_parameters(nparam,idum,nchrome,microga,npossum, &
		& ig2sum,parmin,parmax,pardel,nposibl,g0,g1,ig2,outfp)
	!--------------------------------------------------------------------
	! This subroutine checks input variables and resets them if necessary
	!--------------------------------------------------------------------

		!! subroutine paramters
		integer, intent(in)	:: nparam
		integer, intent(in) 	:: microga
		integer, intent(inout) 	:: idum
		integer, intent(in) 	:: outfp
		integer, intent(out)	:: npossum
		integer, intent(out)	:: ig2sum
		integer, intent(out)	:: nchrome
		real,dimension(nparam),intent(in)		:: parmin
		real,dimension(nparam),intent(in)		:: parmax
		real,dimension(nparam),intent(out)		:: pardel
		real,dimension(nparam),intent(out)		:: g0
		real,dimension(nparam),intent(out)		:: g1
		integer,dimension(nparam),intent(out)	:: nposibl
		integer,dimension(nparam),intent(out) 	:: ig2

		!! local variables
		integer	:: i,j,n2j
		real	:: rand
		integer,parameter :: NMAX=30  !! max number of binary digits/parameter

		!! set minimum value and increment for each parameter
		g0=parmin
		pardel=parmax-parmin
		g1=pardel/dble(nposibl-1)

		!! check that enough posibilities exist given the range of the data
		do i=1,nparam
			do j=1,NMAX
				n2j=2**j
				if (n2j >= nposibl(i)) then
					ig2(i)=j
					exit
				end if
			end do
		end do

		!! count the total number of chromosomes (bits) required
		nchrome=0
		npossum=0
		ig2sum=0
		do i=1,nparam
			nchrome=nchrome+ig2(i)
			npossum=npossum+nposibl(i)
			ig2sum=ig2sum+(2**ig2(i))
		end do
!		if (npossum < ig2sum .and. microga /= 0) then
!			if(verb.gt.5) write(0,2100)
!			write(outfp,2100)
!		end if

		!! initialize random number generator
		call ran3(idum,rand)

!		include 'formats'
1600	format(1x,'ERROR: npopsiz > indmax.  Set indmax = ',i6)
1700	format(1x,'ERROR: nparam > nparmax.  Set nparmax = ',i6)
1800	format(1x,'ERROR: nchrome > nchrmax.  Set nchrmax = ',i6)
2100	format(1x,'WARNING: for some cases, a considerable performance'/ &
			&	1x,'   reduction has been observed when running a non-'/ &
			&	1x,'   optimal number of bits with the micro-GA.'/ &
			&	1x,'   If possible, use values for nposibl of 2**n,'/ &
			&	1x,'   e.g. 2, 4, 8, 16, 32, 64, etc.  See ReadMe file.')
1000	format('*** Parameter adjustment to individual     ',i4, &
			&	', parameter  ',i3,' ***')
1500 	format(i5,3x,30i2)

	end subroutine initialize_parameters

	!--------------------------------------------------------------------
	subroutine initialize_population(istart,npossum,ig2sum,nparam, &
		&	nowrite,microga,idum,irestrt,npopsiz,nchrome, outfp, &
		&	parent,iparent,nposibl,g0,g1,ig2,parmax)
	!--------------------------------------------------------------------
	! This subroutine sets up the program by generating the g0, g1 and
	! ig2 arrays, and counting the number of chromosomes required for the
	! specified input.  The subroutine also initializes the random number
	! generator, parent and iparent arrays (reads the ga.restart file).
	!--------------------------------------------------------------------

		!! subroutine parameters
		integer, intent(in)		:: nparam
		integer, intent(out)	:: istart
		integer, intent(in)		:: outfp  
		integer, intent(in)		:: microga
		integer, intent(inout)		:: idum
		integer, intent(in)		:: irestrt
		integer, intent(inout)	:: npopsiz
		integer, intent(out)	:: nchrome
		integer, intent(out)	:: npossum
		integer, intent(out)	:: ig2sum
		integer, intent(in)		:: nowrite
		integer,dimension(nparam),intent(out)			:: nposibl
		real,dimension(nparam),intent(out)				:: parmax
		real,dimension(nparam),intent(out)				:: g0
		real,dimension(nparam),intent(out)				:: g1
		real,dimension(nparam,npopsiz),intent(out)		:: parent
		integer,dimension(nchrome,npopsiz),intent(out)	:: iparent
		integer,dimension(nparam),intent(out)			:: ig2

		!! local variables
		integer :: i,j,stat,k,l,b
		real 	:: rand

		if (irestrt == 0) then

			!! initialize the random distribution of parameters in the 
			!! individual's parents.
			istart=1
			do i=1,npopsiz
				do j=1,nchrome
					call ran3(idum,rand)
					iparent(j,i)=1
					if(rand < 0.5d0) iparent(j,i)=0
				end do
			end do
			if (npossum < ig2sum) then
				call possibl(nparam,nchrome,npopsiz,nowrite,&
					&	parent,iparent,g0,g1,ig2,nposibl,parmax)
			end if

		else
			!! if irestrt.ne.0, read from restart file.
			open (unit=25, file='ga.restart', status='old', action='read', &
				& iostat=stat)
			if (stat/=0) then
				write(0,*)'ERROR: could not open restart file ga.restart'
				close(25)
				stop
			end if
			rewind 25
			read(25,*) istart,npopsiz
			do j=1,npopsiz
				read(25,*) k,(iparent(l,j),l=1,nchrome)
			end do
			close(25)
		end if
                b=idum-istart
		if(irestrt /= 0) call ran3(b,rand)

!		include 'formats'
1600	format(1x,'ERROR: npopsiz > indmax.  Set indmax = ',i6)
1700	format(1x,'ERROR: nparam > nparmax.  Set nparmax = ',i6)
1800	format(1x,'ERROR: nchrome > nchrmax.  Set nchrmax = ',i6)
2100	format(1x,'WARNING: for some cases, a considerable performance'/ &
			&	1x,'   reduction has been observed when running a non-'/ &
			&	1x,'   optimal number of bits with the micro-GA.'/ &
			&	1x,'   If possible, use values for nposibl of 2**n,'/ &
			&	1x,'   e.g. 2, 4, 8, 16, 32, 64, etc.  See ReadMe file.')
1000	format('*** Parameter adjustment to individual     ',i4, &
			&	', parameter  ',i3,' ***')
1500 	format(i5,3x,30i2)

	end subroutine initialize_population

	!--------------------------------------------------------------------
	subroutine possibl(nparam,nchrome,npopsiz,nowrite, &
		&	array,iarray,g0,g1,ig2,nposibl,parmax)
	!--------------------------------------------------------------------
	! This subroutine determines whether or not all parameters are within
	! the specified range of possibility. If not, the parameter is
	! randomly reassigned within the range. This subroutine is only
	! necessary when the number of possibilities per parameter is not
	! optimized to be 2**n, i.e. if npossum < ig2sum.
	!--------------------------------------------------------------------

		!! subroutine parameters
		integer, intent(in) :: npopsiz
		integer, intent(in) :: nparam
		integer, intent(in) :: nchrome
		integer, intent(in) :: nowrite
		real,dimension(nparam),intent(in)				:: g0
		real,dimension(nparam),intent(in)				:: g1
		integer,dimension(nparam),intent(in)			:: ig2
		integer,dimension(nchrome,npopsiz),intent(inout)	:: iarray
		real,dimension(nparam,npopsiz),intent(out)		:: array
		real,dimension(nparam),intent(in)				:: parmax
		integer,dimension(nparam),intent(in)			:: nposibl

		!! local variables
		integer :: i,j,irand,n2ig2j,b
		real :: rand
                b=1
		do i=1,npopsiz
			call decode(i,nparam,npopsiz,nchrome,array,iarray,g0,g1,ig2)
			do j=1,nparam
				n2ig2j=2**ig2(j)
				if(nposibl(j) /= n2ig2j .and. array(j,i) > parmax(j)) then
					call ran3(b,rand)
					irand=dint(dble(nposibl(j))*rand)
					array(j,i)=g0(j)+dble(irand)*g1(j)
					call code(i,j,nparam,npopsiz,nchrome,array,iarray,&
						&	g0,g1,ig2)
!					if (nowrite /= 0) write(0,1000) i,j
!					if (nowrite /= 0) write(24,1000) i,j
				end if
			end do
		end do

!		include 'formats'
1600	format(1x,'ERROR: npopsiz > indmax.  Set indmax = ',i6)
1700	format(1x,'ERROR: nparam > nparmax.  Set nparmax = ',i6)
1800	format(1x,'ERROR: nchrome > nchrmax.  Set nchrmax = ',i6)
2100	format(1x,'WARNING: for some cases, a considerable performance'/ &
			&	1x,'   reduction has been observed when running a non-'/ &
			&	1x,'   optimal number of bits with the micro-GA.'/ &
			&	1x,'   If possible, use values for nposibl of 2**n,'/ &
			&	1x,'   e.g. 2, 4, 8, 16, 32, 64, etc.  See ReadMe file.')
1000	format('*** Parameter adjustment to individual     ',i4, &
			&	', parameter  ',i3,' ***')
1500 	format(i5,3x,30i2)

	end subroutine possibl


	!--------------------------------------------------------------------
	subroutine decode(i,nparam,npopsiz,nchrome,array,iarray,g0,g1,ig2)
	!--------------------------------------------------------------------
	! This routine decodes a binary string to a real number.
	!--------------------------------------------------------------------

		!! subroutine arguments
		integer, intent(in) :: i
		integer, intent(in) :: nparam
		integer, intent(in) :: npopsiz
		integer, intent(in) :: nchrome
		real,dimension(nparam),intent(in)				:: g0
		real,dimension(nparam),intent(in)				:: g1
		integer,dimension(nparam),intent(in)			:: ig2
		integer,dimension(nchrome,npopsiz),intent(in)	:: iarray
		real,dimension(nparam,npopsiz),intent(out)		:: array

		!! local variables
		integer :: j,k,l,m,iparam

		l=1
		do k=1,nparam
			iparam=0
			m=l
			do j=m,m+ig2(k)-1
				l=l+1
				iparam=iparam+iarray(j,i)*(2**(m+ig2(k)-1-j))
			end do
			array(k,i)=g0(k)+g1(k)*dble(iparam)
		end do

	end subroutine decode

	!--------------------------------------------------------------------
	subroutine code(j,k,nparam,npopsiz,nchrome,array,iarray,g0,g1,ig2)
	!--------------------------------------------------------------------
	! This routine codes a parameter into a binary string.
	!--------------------------------------------------------------------

		!! subroutine parameters
		integer, intent(in) :: j
		integer, intent(in) :: k
		integer, intent(in) :: npopsiz
		integer, intent(in) :: nparam
		integer, intent(in) :: nchrome
		real,dimension(nparam),intent(in) 				:: g0
		real,dimension(nparam),intent(in) 				:: g1
		integer,dimension(nparam),intent(in)			:: ig2
		integer,dimension(nchrome,npopsiz),intent(out)	:: iarray
		real,dimension(nparam,npopsiz),intent(in)		:: array

		!! local variables
		integer :: i,m,iparam,istart,two2m

		!! check for bad input
		if (g1(k) == 0.0d0) return

		!! first, establish the beginning location of the parameter string of
		!! interest.
		istart=1
		do i=1,k-1
			istart=istart+ig2(i)
		end do

		!! find the equivalent coded parameter value, and back out the binary
		!! string by factors of two.
		m=ig2(k)-1
		two2m=2**m
		iparam=NINT((array(k,j)-g0(k))/g1(k))
		do i=istart,istart+ig2(k)-1
			iarray(i,j)=0
			if ((iparam+1) > two2m) then
				iarray(i,j)=1
				iparam=iparam-two2m
			end if
			m=m-1
		end do
		!! write(0,*)array(k,j),iparam,(iarray(i,j),i=istart,istart+ig2(k)-1)

	end subroutine code



	!--------------------------------------------------------------------
	subroutine restart(i,istart,kount,kountmx,maxgen,npopsiz,nchrome,iparent)
	!--------------------------------------------------------------------
	! This subroutine writes information already processed to the 
	! ga.restart file.
	!--------------------------------------------------------------------

		!! subroutine parameters
		integer, intent(in) :: i 
		integer, intent(in) :: istart 
		integer, intent(inout) :: kount 
		integer, intent(in) :: kountmx 
		integer, intent(in) :: npopsiz 
		integer, intent(in) :: nchrome 
		integer, intent(in) :: maxgen 
		integer, dimension(nchrome,npopsiz), intent(out):: iparent

		!! local variables
		integer :: j,l,stat

		kount=kount+1
		if(i == maxgen+istart-1 .or. kount == kountmx) then
!			open (unit=25, file='ga.restart', status='OLD', action='write', &
!				& iostat=stat)
!			rewind 25
!			write(25,*) i+1,npopsiz
!			do j=1,npopsiz
!				write(25,1500) j,(iparent(l,j),l=1,nchrome)
!			end do
!			close(25)
			kount=0
		end if

!		include 'formats'
1600	format(1x,'ERROR: npopsiz > indmax.  Set indmax = ',i6)
1700	format(1x,'ERROR: nparam > nparmax.  Set nparmax = ',i6)
1800	format(1x,'ERROR: nchrome > nchrmax.  Set nchrmax = ',i6)
2100	format(1x,'WARNING: for some cases, a considerable performance'/ &
			&	1x,'   reduction has been observed when running a non-'/ &
			&	1x,'   optimal number of bits with the micro-GA.'/ &
			&	1x,'   If possible, use values for nposibl of 2**n,'/ &
			&	1x,'   e.g. 2, 4, 8, 16, 32, 64, etc.  See ReadMe file.')
1000	format('*** Parameter adjustment to individual     ',i4, &
			&	', parameter  ',i3,' ***')
1500 	format(i5,3x,30i2)

              end subroutine restart

	!-------------------------------------------------------------------
	subroutine ran3(jdum,rand)
	!-------------------------------------------------------------------
	! Returns a uniform random deviate between 0.0 and 1.0.  Set idum to
	! any negative value to initialize or reinitialize the sequence.
	! This function is taken from W.H. Press', "Numerical Recipes" p. 199.
	! Returns a uniform random deviate between 0.0 and 1.0.  Set idum to
	! any negative value to initialize or reinitialize the sequence.
	! This function is taken from W.H. Press', "Numerical Recipes" p. 199.
	!---------------------------------------------------------------------
	! According to Knuth, any large mbig, and any smaller (but still large)
	! mseed can be substituted for the values below.
	!---------------------------------------------------------------------

		!! subroutine parameters
		integer, intent(inout) :: jdum
		real, intent(out) :: rand

		!! local variables
!		double precision,save :: mbig,mseed,mz,fac,mj,mk
!		double precision, dimension(55),save :: ma
!		parameter (mbig=4000000.,mseed=1618033.,mz=0.,fac=1./4000000)
!		parameter (mbig=4000000.,mseed=1618033.,mz=0.,fac=1./mbig)
!               parameter (mbig=1000000000,mseed=161803398,mz=0,fac=1./mbig)
                integer,parameter::double=8
                real(kind=double),parameter::mbig=4000000.,mseed=1618033.,mz=0.,fac=0.00000025
		real(kind=double),save :: ma(55),mj,mk


		integer,save :: i,ii,k,inext,inextp,iff

		data iff / 0 /
		if (jdum < 0 .or. iff == 0) then
			iff=1
			mj=mseed-dble(iabs(jdum))
			mj=dmod(mj,mbig)
			ma(55)=mj
			mk=1
			do i=1,54
				ii=mod(21*i,55)
				ma(ii)=mk
				mk=mj-mk
				if(mk < mz) mk=mk+mbig
				mj=ma(ii)
			end do
			do k=1,4
				do i=1,55
					ma(i)=ma(i)-ma(1+mod(i+30,55))
					if(ma(i) < mz) ma(i)=ma(i)+mbig
				end do
			end do
			inext=0
			inextp=31
			jdum=1
		endif
		inext=inext+1
		if(inext == 56) inext=1
		inextp=inextp+1
		if(inextp == 56) inextp=1
		mj=ma(inext)-ma(inextp)
		if(mj < mz) mj=mj+mbig
		ma(inext)=mj
		rand=mj*fac

	end subroutine ran3
	!-----------------------------------------------------------------------

end module ga_aux
