program unresag
      
      implicit none
      include 'GA.inc'
      include 'common.inc'
      include 'io.inc'
      real*8 :: WEIGHTS(18)
      integer*4 :: i,j,m,n,iloc,pos,iter,maxiter,numarg,tmplen
      integer*4 :: WybierzOsobnika,Najlepszy,ZnajdzPodobnego,FindWorst
      integer*4 :: status
      integer*4 :: mnoznik = 10
      character*16 :: mode,argn,argi,dir_nam
      character*15 :: filename
      real*8 :: loc,fc,sumfitness,sump,avrp,maxfc
      real :: rand, r
      real*8 :: maxfcever
      real*8 :: FunkcjaCeluB, FunkcjaCeluB2, FunkcjaCeluR
      character(80),dimension(5) :: wagi
      logical :: searching = .false.
      logical :: logdata = .true.
      logical :: fileex = .false.
      logical :: retry = .false.
      logical :: dir_ex = .false. 
      logical :: debug = .true.
      character*200 :: text, tmptext,tmptext2
      character(len=*), parameter :: FMT = "(19F10.5,E15.7,F10.5)"
      
c =======================================================================
c  Main program
c =======================================================================


      call ReadInput(status)
      if ((maxminstep.gt.0).and.(generation.eq.0)) then
c Do first score  
       do_fs=.true.
      endif
c  Is everything  OK?
      if (status.gt.0) then
       write(*,*) "There ware input file errors - check log."
       goto 2010
      endif

      call itime(timeArray)
      r = rand(timeArray(1)*10000+timeArray(2)*100+timeArray(3))
c     Fair enough?
c      r = 527865426


c-------------------------------------------------------
c bank & populacja arrays store this   
c 1-19 - weights
c   20 - energy obtained from zscore
c   21 - fittness 
c-------------------------------------------------------
      allocate (bank(banksize,21))
      allocate (populacja(BANK_MULTIPLIER*banksize,21))
      allocate (pairs(BANK_MULTIPLIER*banksize))
c
c Zero the arrays
c
      do i=1,BANK_MULTIPLIER*banksize
       do j=1,21
        populacja(i,j)=0
       enddo
       pairs(i)=0
      enddo
      do i=1,banksize
       do j=1,21
       bank(i,j)=0
       enddo
      enddo

c First time here?
      inquire(FILE=ostatefn,EXIST=fileex)
c No. Prepere next generation
      if (fileex) then
       open(ostate,file=ostatefn) 
       read(ostate,'(I4)') generation
       read(ostate,'(L2)') do_optima
       read(ostate,'(L2)') do_ga
       read(ostate,'(F7.5)') avrd
       close(ostate)
       if (do_ga) then
        call write2log("Doing GA in this step")
       endif
       if (.not.do_optima) then
        maxminstep=0
        call write2log("ZSCORE weights minimalization disabled for now")
        generation=generation+1
       endif
       write(tmptext,'(I4)') generation 
       call write2log("This is genaration "//tmptext)
       call ReadOptimaW(BANK_MULTIPLIER*banksize,populacja)

c Yes. Generate random zero-population
      else
       call GenPopulation(BANK_MULTIPLIER*banksize,populacja)
       write(tmptext,'(I4)') generation 
       call write2log("This is genaration "//tmptext)
       if (do_optima) then
        do_ga=.false.
       endif 

c End of "first time here?" code
      endif

c
c Do we actual use a Genetic Algorithm?
c

      if (do_ga.eqv..true.) then

c Yes. This is only done when doing second pass with ZSCORE (without weight minimalizaton)
c or with weight minimalization disabled (maxmin=0) just to obtain final score)

c \\---//================================================================c
c  \\-//                                                                 c
c   "//  Genetic Algorithm code starts here                              c
c   //-\\                                                                c
c =//---\\===============================================================c

      call ReadZEnergy(BANK_MULTIPLIER*banksize,populacja)
      do i=1,BANK_MULTIPLIER*banksize
c       populacja(i,20)=FunkcjaCeluR(populacja(i,20))
       write(*,*) "Ind ",i,populacja(i,20)
      enddo

      call write2log("==[ Doing GA now ]=============================")

c
c Populate bank depending on algorithm
c
      select case(alg)
       case('simple')
        call GetNBest(populacja,bank,banksize)
       case('better')
        call GetNBest(populacja,bank,banksize)
       case('betterv2')
        call GetNBest(populacja,bank,banksize)
       case('csa')
c --- debug begin ---
	if (debug) then
	 write(*,*) "generation: ",generation
	 write(*,*) "maxminstep: ",maxminstep
	 write(*,*) "do_ga: ",do_ga
	 write(*,*) "do_optima: ",do_optima
 	 write(*,*) "do_fs: ",do_fs
c	 write(*,*) "cicutoff: ", cicutoff         
	endif
c --- debug end ---

c        csacutoff=cicutoff

c
c  Fill the bank just after the first time we get the score 
c
        if (((generation.eq.1).and.(maxminstep.eq.0)).or.((generation.eq&
     &.1).and.(maxminstep.gt.0))) then

c --- debug begin ---
	 if (debug) then 
	  write(*,*) "Calling GetNBest..."
	 endif
c --- debug end ---
         call GetNBest(populacja,bank,banksize)
         call CalcAvgDist(bank,avrd)
         write(tmptext,'(F7.5)') avrd
         call write2log("Average distance between individuals in initial&
     & bank is "//trim(tmptext))
c
c Cutoff 
c
         csacutoff=(maxco*avrd)-generation*avrd*(maxco-minco)/maxgen
         write(tmptext,'(F7.5)') csacutoff
         call write2log("CSA cutoff is now set to "//trim(tmptext))
c         csacutoff=maxco*avrd         

         

c --- debug begin ---
	 if (debug) then
	  do i=1,banksize
	   write(*,FMT) bank(i,:)
	  enddo
	 endif
c --- debug end ---

         call CalcFitness(bank,sumfitness)
         call WriteBank(bank)  

c
c Do normal CSA
c        
        else

c --- debug begin ---
	 if (debug) then 
	  write(*,*) "Szukam podobnego"
	 endif
c --- debug end ---

         call ReadBank(bank)
         write(tmptext,'(F7.5)') avrd
         call write2log("Average distance in bank is "//trim(tmptext))

         csacutoff=maxco*avrd-generation*avrd*(maxco-minco)/maxgen
         write(tmptext,'(F7.5)') csacutoff
         call write2log("CSA cutoff is now set to "//trim(tmptext))
  
         do i=1,BANK_MULTIPLIER*banksize
          write(tmptext,'(I4)') i
          call write2log("Checking ind "//trim(tmptext))  
          j=ZnajdzPodobnego(banksize,bank,populacja(i,:),csacutoff)
          if (j.gt.0) then  ! W banku jest podobny
           if (populacja(i,20).lt.bank(j,20)) then
            write(tmptext,'(I4)') j
            write(tmptext2,'(I4)') i
            call write2log("  Swaping ind"//trim(tmptext)//" from bank t&
     &o ind "//trim(tmptext2)//" from population")
            write(tmptext2,'(19F8.5,E15.7,F8.5)') bank(j,:)
            call write2log("  BANK"//trim(tmptext)//":"//trim(tmptext2))

            write(tmptext,'(I4)') i
            write(tmptext2,'(19F8.5,E15.7,F8.5)') populacja(i,:)
            call write2log("  POP "//trim(tmptext)//":"//trim(tmptext2))
      
            bank(j,:)=populacja(i,:)
           else
            call write2log("  Found simialar but not better")
           endif
          else              ! W banku nie ma podobnego 
           j=FindWorst(banksize,bank)
           write(tmptext2,'(I4)') j
           if (populacja(i,20).lt.bank(j,20)) then
            call write2log("  Worst in bank is "//trim(tmptext2))
            write(tmptext,'(I4)') i
            call write2log("  Swaping worst ind in bank to "//trim(tmpte&
     &xt))
            write(tmptext,'(I4)') j
            write(tmptext2,'(19F8.5,E15.7,F8.5)') bank(j,:)
            call write2log("  BANK"//trim(tmptext)//":"//trim(tmptext2))
            write(tmptext,'(I4)') i
            write(tmptext2,'(19F8.5,E15.7,F8.5)') populacja(i,:)
            call write2log("  POP "//trim(tmptext)//":"//trim(tmptext2))
            bank(j,:)=populacja(i,:)
           else
            call write2log("  The worst in bank is better then this Ind"&
     &)
           endif
          endif
         enddo

c -----------------------------------------------------------------------
c      Calculate fitness for reproduction
c -----------------------------------------------------------------------
          call CalcFitness(bank,sumfitness)

c --- debug begin ---
	if (debug) then
	 write (*,*) "Dumping bank: "        
	 do i=1,banksize
	  write(*,FMT) bank(i,:)
	 enddo
	 write (*,*) "Sumfitness: ", sumfitness   
	endif
c --- debug end ---

         call WriteBank(bank)

        endif
       case('cluster')
        write(*,*) "Well this is not implemented yet"
        goto 2010
      end select

c      call WritePopSum()


c Have we done last generation ?
      if (generation.eq.maxgen) then 
       call write2log("Maximum number of genarations reached. QUITING.")
       goto 2010
      endif
c Prapering for  next generation so increase counter 
        iter=iter+1
 
c -----------------------------------------------------------------------
c	Reproduction
c -----------------------------------------------------------------------
      call write2log("==[ Reproducting ]============================")
      select case(alg)
       case('simple','csa')
        do i=1,BANK_MULTIPLIER*banksize
        m=WybierzOsobnika(banksize,bank)
        write(tmptext,'(I4)') m
        write(tmptext2,'(I4)') i
        call write2log("Reproducting individual "//trim(tmptext)//" from&
     & bank as new individual "//trim(tmptext2))
        populacja(i,:)=bank(m,:)
        end do
       case('better','betterv2')
        populacja(1,:)=bank(1,:)
        call write2log("Reproducting first (best) individual")
        do i=2,BANK_MULTIPLIER*banksize
        m=WybierzOsobnika(banksize,bank)
        write(tmptext,'(I4)') m
        write(tmptext2,'(I4)') i
        call write2log("Reproducting individual "//trim(tmptext)//" from&
     & bank as new individual "//trim(tmptext2))
        populacja(i,:)=bank(m,:)
        end do
      end select

c -----------------------------------------------------------------------
c	Pairing
c -----------------------------------------------------------------------
      do i=1,BANK_MULTIPLIER*banksize
       pairs(i)=0
      end do
      
      do i=1,BANK_MULTIPLIER*banksize
       if (pairs(i).eq.0) then
210     pos=1+int(rand(0)/(1/real(BANK_MULTIPLIER*banksize)))
        if (pos.le.i) then
         goto 210
        endif
        do j=1,(i-1)
         if (pairs(j).eq.pos) then
          goto 210
         endif
        end do
        pairs(i) = pos
        pairs(pos) = i
       endif
      end do  
      call write2log("==[ Pairing ]==============================")
      do i=1,BANK_MULTIPLIER*banksize
       write(tmptext,'(A11,I3.3,A5,I3.3,A21)') "Individual ",i," and ",p&
     &airs(i)," will be crossed over"
       call write2log(trim(tmptext))
      end do

c -----------------------------------------------------------------------
c 	Crossover
c -----------------------------------------------------------------------
      call write2log("==[ Crossing ]==============================")
      do i=1,BANK_MULTIPLIER*banksize
       if (pairs(i) .gt. i) then
        call CrossoverS(populacja(i,:),populacja(pairs(i),:))
       endif
      end do
       
      
       call write2log("         WLONG   WSCP    WELEC   WBOND   WANG    &
     &WSCLOC  WTOR    WTORD   WCORRH  WCORR4  WCORR5  WCORR6  WEL_LOC WT&
     &URN3  WTURN4  WTURN6  WVDWPP  WHPB    WSCCOR")
      do i=1,BANK_MULTIPLIER*banksize
       tmplen = 9
       text = ""
       write(text,'(A4,I3.3,A2)') "Ind",i,": "
       do j=1,19
        write(tmptext,'(F7.5)') populacja(i,j)
        text =  text(1:tmplen) // tmptext
        tmplen = tmplen + 8
       end do
       call write2log(trim(text))
      end do

c -----------------------------------------------------------------------
c 	Mutation
c -----------------------------------------------------------------------

      call write2log("==[ Mutating ]==============================") 
      do i = 1,19
       if (IGNORE_WEIGHT(i).eq.1) then
        write(tmptext,'(I2)') i
        call write2log("Parameter "//trim(WNAME(i))//" will be ignored")
       endif
      end do  
      do i=1,BANK_MULTIPLIER*banksize
       call MutationV(populacja(i,:),i,MUTRATIO)
      end do
c for betterv2 and csa restore best individual
      select case(alg)
       case('betterv2','csa')
       populacja(1,:)=bank(1,:)
      end select

c
c Genetic Algorithm blokcode stops here
c
      endif

c ---------------------------------
c  Do some final stuff 
c ---------------------------------
c
c Set the variables before writing to state file 
c
      if ((.not.do_ga).and.(.not.do_optima)) then
       if (do_fs) then
        do_ga=.false.
        do_optima=.true.
       else
        do_ga=.true.
        do_optima=.false.
       endif
      else
       if (do_fs) then
        do_optima=.not.do_optima
        if (generation.eq.0) then
         do_ga=.false.
        else
         do_ga=.not.do_ga
        endif
       else
        do_ga=.true.
        do_optima=.false.
       endif 
      endif


c
c Write te current state and population
c
      call WriteState()
c      call WritePopSum()
c
c Create the inputs
c 
      write(tmptext,'(I3)') generation+1
      call write2log("Preparing inputs for next generation ("//trim(adju&
     &stl(tmptext))//")")
      call CreateInputs(BANK_MULTIPLIER*banksize,populacja)
c
c All done? Then let's quit.
c
      goto 2010



c ======================================================================
c
c    MAIN CODE ENDS HERE
c
c ======================================================================



2000  deallocate(populacja)
      deallocate(pairs)
2010  end program


c ======================================================================
c  Clustering subroutines
c ======================================================================

      include 'cluster.inc'


c ======================================================================
c  Write2log subroutine
c ======================================================================
c  Creates logs  
c ----------------------------------------------------------------------

      subroutine Write2log(text)
       include 'common.inc'
       include 'io.inc'
       character*200 :: tmptext
       character(len=*) :: text
       logical :: ex
       call itime(timeArray)
       inquire(file=ologfn,exist=ex)
       if (.not.ex) then
       open(olog, file=ologfn, access="append")
       write(olog,"(A)") "==============================================&
     &=========="
       tmptext = "= UNRES weights optimizer version "//version 
       write(olog,"(A)") tmptext
       write(olog,"(A)") info
       write(olog,"(A)") "==============================================&
     &=========="
       close(olog)
       endif
       open(olog, file = ologfn , access = "append")
       write(olog,"(I2.2,A1,I2.2,A1,I2.2,A2,A)") timeArray(1),":",timeAr&
     &ray(2),":",timeArray(3),": ",text
       close(olog)
      end subroutine Write2log

c ======================================================================
c  GetNBest subroutine
c ======================================================================
c Fills up the bank population up to banksize with individuals from
c inputpop having the best(lowest) score
c ----------------------------------------------------------------------

      subroutine GetNBest(inputpop,bank,banksize)
       integer*4 :: banksize, i,j,idx
       real*8 :: inputpop, bank, best, last
       include 'GA.inc'
       dimension inputpop(banksize*BANK_MULTIPLIER,21)
       dimension bank(banksize,21)
       
       last=0.0
       idx=1
       do j=1,banksize
       best=huge(0.0d0)
        do i=1,banksize*BANK_MULTIPLIER
         if ((inputpop(i,20).lt.best).and.(inputpop(i,20).gt.last)) then
          best=inputpop(i,20)
          idx=i
         endif
        end do
        bank(j,:)=inputpop(idx,:)
        last=best
       end do
      end subroutine GetNBest

c ======================================================================
c  FindWorst subroutine
c ======================================================================
c  Finds the worst (highest score) individual in population
c ----------------------------------------------------------------------

      integer*4 function FindWorst(rozmiar,pop)
       integer*4 :: rozmiar, i, idx
       real*8 :: pop,last
       dimension pop(rozmiar,21)
       
       last=0.0
       idx = 0
       do i=1,rozmiar
        if (pop(i,20).gt.last) then
         idx = i
         last=pop(i,20)
        endif
       end do
       FindWorst=idx
       return
      end function

c ======================================================================
c  GenPopulation subroutine
c ======================================================================
c  Generates a population (pop) of nind individuals with 19 random 
c  wights from a set of ranges (see GA.inc).
c ----------------------------------------------------------------------

      subroutine GenPopulation(nind,pop)
       integer*4 :: i,j,nind
       real*8 :: pop
       dimension pop(nind,21) 
       include 'GA.inc'
       do j=1,nind
        do i=1,19
         select case(i)
          case (1)
           pop(j,1)=WLONG_LOW+(rand(0)*(WLONG_HI-WLONG_LOW))
          case (2)
           pop(j,2)=WSCP_LOW+(rand(0)*(WSCP_HI-WSCP_LOW))
          case (3)
           pop(j,3)=WELEC_LOW+(rand(0)*(WELEC_HI-WELEC_LOW))
          case (4)
           pop(j,4)=WBOND_LOW+(rand(0)*(WBOND_HI-WBOND_LOW))
          case (5)
           pop(j,5)=WANG_LOW+(rand(0)*(WANG_HI-WANG_LOW))
          case (6)
           pop(j,6)=WSCLOC_LOW+(rand(0)*(WSCLOC_HI-WSCLOC_LOW))
          case (7)
           pop(j,7)=WTOR_LOW+(rand(0)*(WTOR_HI-WTOR_LOW))
          case (8)
           pop(j,8)=WTORD_LOW+(rand(0)*(WTORD_HI-WTORD_LOW))
          case (9)
           pop(j,9)=WCORRH_LOW+(rand(0)*(WCORRH_HI-WCORRH_LOW))
          case (10)
           pop(j,10)=WCORR4_LOW+(rand(0)*(WCORR4_HI-WCORR4_LOW))
          case (11)
           pop(j,11)=WCORR5_LOW+(rand(0)*(WCORR5_HI-WCORR5_LOW))
          case (12)
           pop(j,12)=WCORR6_LOW+(rand(0)*(WCORR6_HI-WCORR6_LOW))
          case (13)
           pop(j,13)=WEL_LOC_LOW+(rand(0)*(WEL_LOC_HI-WEL_LOC_LOW))
          case (14)
           pop(j,14)=WTURN3_LOW+(rand(0)*(WTURN3_HI-WTURN3_LOW))
          case (15)
           pop(j,15)=WTURN4_LOW+(rand(0)*(WTURN4_HI-WTURN4_LOW))
          case (16)
           pop(j,16)=WTURN6_LOW+(rand(0)*(WTURN6_HI-WTURN6_LOW))
          case (17)
           pop(j,17)=WVDWPP_LOW+(rand(0)*(WVDWPP_HI-WVDWPP_LOW))
          case (18)
           pop(j,18)=WHPB_LOW+(rand(0)*(WHPB_HI-WHPB_LOW))
          case (19)
           pop(j,19)=WSCCOR_LOW+(rand(0)*(WSCCOR_HI-WSCCOR_LOW))
         end select
        end do
       end do
      end subroutine GenPopulation
     
c ======================================================================
c  CrossoverS subroutine
c ======================================================================
c  Does a simple value crossing-over on a pair
c  of two given individuals, giving two new individuals.
c ----------------------------------------------------------------------

      subroutine CrossoverS(ind1,ind2)
       real*8 :: ind1(19),ind2(19),temp(19)
       integer*4 :: loc
       loc = 1 + int(rand(0)/(1/18.0))
       temp = ind2
       do i=(loc),19
        ind2(i)=ind1(i)
        ind1(i)=temp(i)
       end do 
      end subroutine CrossoverS
     
c ======================================================================
c  MutationV subroutine
c ======================================================================
c  Does a value mutation for the given individual (ind) with the 
c  given mutation ratio (mratio).
c ----------------------------------------------------------------------

      subroutine MutationV(ind,id,mratio)
       include 'GA.inc'
c       real*8 :: rand
       real :: mratio
       real*8 :: ind(19),ti
       integer*4 :: i,id
       character*100 :: tmptext
       character*150 :: logtext
       
       do i = 1,19
        if (IGNORE_WEIGHT(i).eq.1) then
c            write(tmptext,'(I2)') i
c	    call write2log("Skipping weight "//trim(tmptext))
            goto 2300
        endif 
        if (rand(0)<mratio) then
          select case(i)
            case (1)
             ti=ind(1)
             ind(1)=WLONG_LOW+(rand(0)*(WLONG_HI-WLONG_LOW))
            case (2)
             ti=ind(2)
             ind(2)=WSCP_LOW+(rand(0)*(WSCP_HI-WSCP_LOW))
            case (3)
             ti=ind(3)
             ind(3)=WELEC_LOW+(rand(0)*(WELEC_HI-WELEC_LOW))
            case (4)
             ti=ind(4)
             ind(4)=WBOND_LOW+(rand(0)*(WBOND_HI-WBOND_LOW))
            case (5)
             ti=ind(5)
             ind(5)=WANG_LOW+(rand(0)*(WANG_HI-WANG_LOW))
            case (6)
             ti=ind(6)
             ind(6)=WSCLOC_LOW+(rand(0)*(WSCLOC_HI-WSCLOC_LOW))
            case (7)
             ti=ind(7)
             ind(7)=WTOR_LOW+(rand(0)*(WTOR_HI-WTOR_LOW))
            case (8)
             ti=ind(8)
             ind(8)=WTORD_LOW+(rand(0)*(WTORD_HI-WTORD_LOW))
            case (9)
             ti=ind(9)
             ind(9)=WCORRH_LOW+(rand(0)*(WCORRH_HI-WCORRH_LOW))
            case (10)
             ti=ind(10)
             ind(10)=WCORR4_LOW+(rand(0)*(WCORR4_HI-WCORR4_LOW))
            case (11)
             ti=ind(11)
             ind(11)=WCORR5_LOW+(rand(0)*(WCORR5_HI-WCORR5_LOW))
            case (12)
             ti=ind(12)
             ind(12)=WCORR6_LOW+(rand(0)*(WCORR6_HI-WCORR6_LOW))
            case (13)
             ti=ind(13)
             ind(13)=WEL_LOC_LOW+(rand(0)*(WEL_LOC_HI-WEL_LOC_LOW))
            case (14)
             ti=ind(14)
             ind(14)=WTURN3_LOW+(rand(0)*(WTURN3_HI-WTURN3_LOW))
            case (15)
             ti=ind(15)
             ind(15)=WTURN4_LOW+(rand(0)*(WTURN4_HI-WTURN4_LOW))
            case (16)
             ti=ind(16)
             ind(16)=WTURN6_LOW+(rand(0)*(WTURN6_HI-WTURN6_LOW))
            case (17)
             ti=ind(17)
             ind(17)=WVDWPP_LOW+(rand(0)*(WVDWPP_HI-WVDWPP_LOW))
            case (18)
             ti=ind(18)
             ind(18)=WHPB_LOW+(rand(0)*(WHPB_HI-WHPB_LOW))
            case (19)
             ti=ind(19)
             ind(19)=WSCCOR_LOW+(rand(0)*(WSCCOR_HI-WSCCOR_LOW))
        end select
         write(tmptext,'(I4,A13,I2,A7,F7.5,A5,F7.5)') id," at position "&
     &,i," from: ",ti," to: ",ind(i)
         logtext = "Mutation occured in individual "//tmptext
         call write2log (logtext)
 
        endif
2300   end do
      end subroutine 

c ======================================================================
c  FunkcjaCeluB function
c ======================================================================
c  Simple (dummy) scoring function
c ----------------------------------------------------------------------

      real*8 function FunkcjaCeluB(osobnik)
       include 'TEST.inc'
       real*8 :: osobnik(20)
       integer*4 :: i
       FunkcjaCeluB = 20 -exp(-abs(WLONG-osobnik(1)))                   &
     & - exp(-abs(WSCP-osobnik(2)))                                     &
     & - exp(-abs(WELEC-osobnik(3)))                                    &
     & - exp(-abs(WBOND-osobnik(4)))                                    &
     & - exp(-abs(WANG-osobnik(5)))                                     &
     & - exp(-abs(WSCLOC-osobnik(6)))                                   &
     & - exp(-abs(WTOR-osobnik(7)))                                     &
     & - exp(-abs(WTORD-osobnik(8)))                                    &
     & - exp(-abs(WCORRH-osobnik(9)))                                   &
     & - exp(-abs(WCORR4-osobnik(10)))                                  &
     & - exp(-abs(WCORR5-osobnik(11)))                                  &
     & - exp(-abs(WCORR6-osobnik(12)))                                  &
     & - exp(-abs(WEL_LOC-osobnik(13)))                                 &
     & - exp(-abs(WTURN3-osobnik(14)))                                  &
     & - exp(-abs(WTURN4-osobnik(15)))                                  &
     & - exp(-abs(WTURN6-osobnik(16)))                                  &
     & - exp(-abs(WVDWPP-osobnik(17)))                                  &
     & - exp(-abs(WHPB-osobnik(18)))                                    &
     & - exp(-abs(WSCCOR-osobnik(19)))
       return
      end function

c ======================================================================
c  FunkcjaCeluB2 function
c ======================================================================
c  Simple (dummy) scoring function
c ----------------------------------------------------------------------

      real*8 function FunkcjaCeluB2(osobnik)
       include 'TEST.inc'
       real*8 :: osobnik(20)
       integer*4 :: i
       FunkcjaCeluB2 = 1 + 4*(WLONG-osobnik(1))**2                      &
     & + 4*(WSCP-osobnik(2))**2                                         &
     & + 4*(WELEC-osobnik(3))**2                                        &
     & + 4*(WBOND-osobnik(4))**2                                        &
     & + 4*(WANG-osobnik(5))**2                                         &
     & + 4*(WSCLOC-osobnik(6))**2                                       &
     & + 4*(WTOR-osobnik(7))**2                                         &
     & + 4*(WTORD-osobnik(8))**2                                        &
     & + 4*(WCORRH-osobnik(9))**2                                       &
     & + 4*(WCORR4-osobnik(10))**2                                      &
     & + 4*(WCORR5-osobnik(11))**2                                      &
     & + 4*(WCORR6-osobnik(12))**2                                      &
     & + 4*(WEL_LOC-osobnik(13))**2                                     &
     & + 4*(WTURN3-osobnik(14))**2                                      &
     & + 4*(WTURN4-osobnik(15))**2                                      &
     & + 4*(WTURN6-osobnik(16))**2                                      &
     & + 4*(WVDWPP-osobnik(17))**2                                      &
     & + 4*(WHPB-osobnik(18))**2                                        &
     & + 4*(WSCCOR-osobnik(19))**2
       return
      end function

c ======================================================================
c  Weights2str subroutine
c ======================================================================
c  Writes weights of individual(osobnik) to string (lista)
c ----------------------------------------------------------------------

      subroutine Weights2Str(osobnik,coff,lista)   
      real*8, dimension(19) :: osobnik
      real*8 :: coff
      character(80) :: a
      integer :: i
      character(80) :: lista
      dimension lista(5)
c	Zerowanie tablicy
      do i=1,5
       lista(i)=''
      enddo
c
      !write(*,*) coff
      write(a,"(A6,F7.5,A6,F7.5,A7,F7.5,A7,F7.5,A6,F7.5,A)") "WLONG=",os&
     &obnik(1)," WSCP=",osobnik(2)," WELEC=",osobnik(3)," WBOND=",osobni&
     &k(4)," WANG=",osobnik(5),"            &"
      lista(1)=a
      write(a,"(A7,F7.5,A6,F7.5,A7,F7.5,A8,F7.5,A8,F7.5,A)") "WSCLOC=",o&
     &sobnik(6)," WTOR=",osobnik(7)," WTORD=",osobnik(8)," WCORRH=",osob&
     &nik(9),"WCORR5=",osobnik(11),"        &"
      lista(2)=a       
      write(a,"(A7,F7.5,A9,F7.5,A8,F7.5,A8,F7.5,A8,F7.5,A)") "WCORR6=",o&
     &sobnik(12)," WEL_LOC=",osobnik(13)," WTURN3=",osobnik(14)," WTURN4&
     &=",osobnik(15),"WTURN6=",osobnik(16),"    &"
      lista(3)=a
      write(a,"(A7,F7.5,A6,F7.5,A8,F7.5,A)") "WVDWPP=",osobnik(17)," WHP&
     &B=",osobnik(18)," WSCCOR=",osobnik(19),"                          &
     &           &"
      lista(4)=a
      write(a,"(A7,F7.5,A8,F7.5)") "CUTOFF=",coff," WCORR4=",osobnik(10)
      lista(5)=a
      
      end subroutine

c ======================================================================
c  FunkcjaCeluR function
c ======================================================================
c  Real scoring function
c ----------------------------------------------------------------------

      real*8 function FunkcjaCeluR(zscore)
       real*8 :: zscore
       if (zscore.eq.0) then
        FunkcjaCeluR = 0
       else
        FunkcjaCeluR = 1/zscore
       endif
       return
      end function

c ======================================================================
c  ReadInput subroutine
c ======================================================================

      subroutine ReadInput(status)
      include "io.inc"
      include 'common.inc'
      include 'GA.inc'
      logical :: ex
      integer :: stat
      integer*4 :: status
      character*100 :: wiersz,tmp
      inquire(FILE=inpfn,EXIST=ex) 
      if (ex) then
       status = 0
       call write2log('')
       call write2log('==[ Reading main input ]======================')
       call write2log('')
       open(inp,file=inpfn)
       do
         read(inp, '(A)', iostat=stat) wiersz
         if (stat /= 0) exit
c PDB=
         if ((wiersz(1:4).eq.'PDB=').or.(wiersz(1:4).eq.'pdb=')) then
          npdb=1
          tmp = wiersz(5:len_trim(wiersz))
          do i=1,len_trim(tmp)
           if (tmp(i:i).eq.' ') then
            npdb=npdb+1
           endif
          end do
          if (npdb.gt.maxnpdb) then
           call write2log("Number of input PDB exceeds maxnpdb!")
           status = 1
           exit
          endif
          do i=1,npdb
           if (index(trim(tmp),' ').gt.0) then
            pdbfiles(i)=tmp(1:index(trim(tmp),' '))
           else
            pdbfiles(i)=tmp(1:len_trim(tmp))
           endif
           tmp=tmp(index(trim(tmp),' ')+1:len_trim(tmp))
          end do
         endif ! Koniec czytania "PDB="

c ALG=  
         if ((wiersz(1:4).eq.'ALG=').or.(wiersz(1:4).eq.'alg=')) then
          alg=wiersz(5:len_trim(wiersz))
          select case(alg)
           case('simple')
            call write2log ("Using 'simple' algorithm")
           case('better')
            call write2log ("Using 'better' algorithm")
           case('betterv2')
            call write2log ("Using 'betterv2' algorithm")
           case('csa')
            call write2log ("Using 'CSA' algorithm")
           case('cluster')
            call write2log ("Using 'Clustering CSA' algorithm")
           case default
            alg="csa"
            call write2log ("Unknown algorithm. Using 'csa' as default")
          end select
         endif 

c GENERATIONS=
         select case (wiersz(1:12))
          case('GENERATIONS=','generations=')
           tmp = wiersz(13:len_trim(wiersz))
           read(tmp,'(I4)') maxgen 
         end select 

c CICUTOFF=
c         select case(wiersz(1:9))
c          case('CICUTOFF=','cicutoff=')
c          tmp = wiersz(10:len_trim(wiersz))
c          read(tmp(1:len_trim(tmp)),'(F7.5)') cicutoff
c          call write2log("Initial CSA cutoff is set to "//tmp)
c         end select

c MINCO=
         select case(wiersz(1:6))
          case('MINCO=','minco=')
          tmp = wiersz(7:len_trim(wiersz))
          read(tmp(1:len_trim(tmp)),'(F7.5)') minco
          call write2log("Minimal CSA cutoff factor is set to "//tmp)
         end select

c MAXCO=
         select case(wiersz(1:6))
          case('MAXCO=','maxco=')
          tmp = wiersz(7:len_trim(wiersz))
          read(tmp(1:len_trim(tmp)),'(F7.5)') maxco
          call write2log("Maximal CSA cutoff factor is set to "//tmp)
         end select


c POPULATION= 
         select case(wiersz(1:11))
          case('POPULATION=','population=')
           tmp = wiersz(12:len_trim(wiersz))
           read(tmp,'(I4)') banksize
           call write2log("Bank size is set to "//tmp)
         end select

c WHAMTEMPLATE
         select case(wiersz(1:13))
          case('WHAMTEMPLATE=','whamtemplate=')
           tmp = wiersz(14:len_trim(wiersz))
           ntwham=1
           do i=1,len_trim(tmp)
            if (tmp(i:i).eq.' ') then
             ntwham=ntwham+1
            endif
           end do
           if (ntwham.gt.maxnpdb) then
            call write2log("Number of WHAM templates exceeds maxnpdb!")
             status = 1
             exit
           endif
           do i=1,ntwham
            if (index(trim(tmp),' ').gt.0) then
              whamtemplate(i)=tmp(1:index(trim(tmp),' '))
            else
              whamtemplate(i)=tmp(1:len_trim(tmp))
            endif
              tmp=tmp(index(trim(tmp),' ')+1:len_trim(tmp))
            end do
         end select

c MREMDTEMPLATE=
         select case(wiersz(1:14))
          case('MREMDTEMPLATE=','mremdtemplate=')
          tmp = wiersz(15:len_trim(wiersz))
          ntmremd=1
          do i=1,len_trim(tmp)
           if (tmp(i:i).eq.' ') then
           ntmremd=ntmremd+1
           endif
          end do
          if (ntmremd.gt.maxnpdb) then
            call write2log("Number of MREMD templates exceeds maxnpdb!")
            status = 1
            exit
          endif
          do i=1,ntmremd
            if (index(trim(tmp),' ').gt.0) then
             mremdtemplate(i)=tmp(1:index(trim(tmp),' '))
            else
             mremdtemplate(i)=tmp(1:len_trim(tmp))
            endif
            tmp=tmp(index(trim(tmp),' ')+1:len_trim(tmp))
          end do
         end select

c MAXMIN= 
         select case(wiersz(1:7))
          case('MAXMIN=','maxmin=')
          tmp = wiersz(8:len_trim(wiersz))
          read(tmp,'(I4)') maxminstep
          if (maxminstep.gt.0) then
           call write2log("ZSCORE weights minimalization set to "//tmp)
           do_optima=.true.
          else
           call write2log("ZSCORE weights minimalization disabled")
           do_optima=.false.
          endif
         end select

c SCRIPTS=
         select case(wiersz(1:8))
          case('SCRIPTS=','scripts=')
          nscripts=1
          tmp = wiersz(9:len_trim(wiersz))
          do i=1,len_trim(tmp)
           if (tmp(i:i).eq.' ') then
            nscripts=nscripts+1
           endif
          end do
          if (nscripts.gt.maxscripts) then
           call write2log("Number of scripts exceeds maxscripts!")
           status = 1
           exit
          endif
          do i=1,nscripts
           if (index(trim(tmp),' ').gt.0) then
            scripts(i)=tmp(1:index(trim(tmp),' '))
           else
            scripts(i)=tmp(1:len_trim(tmp))
           endif
           call write2log("Script will be used: "//scripts(i)) 
           tmp=tmp(index(trim(tmp),' ')+1:len_trim(tmp))
          end do
         end select ! Koniec czytania "scripts="
        end do

        close(inp)

c write additional info to log
        write(tmp,'(F8.5)') MUTRATIO
        call write2log("Mutation ratio is set to "//tmp)
        write(tmp,'(I2)') BANK_MULTIPLIER
        call write2log("Bank multiplier is set to "//tmp)
        write(tmp,'(I5)') BANK_MULTIPLIER*banksize
        tmp=adjustl(tmp)
        call write2log("Trail population size will be "//trim(tmp)//" in&
     &dividuals")
      else
        call write2log("Input file unresga.inp does not exist!!")
        write(*,*) "Input file unresga.inp does not exist!!"
        status = 1
      endif

c Checking set parameters after reading input

c PDB= ?
      if (npdb.eq.0) then
        call write2log("Can't find 'PDB=' entry in input file")
        status = 1
      else
       do i=1,npdb
        inquire(FILE=trim(pdbfiles(i)),EXIST=ex)
        if (.not.ex) then
         call write2log("Can't find file "//trim(pdbfiles(i)))
         status = 1
        endif
       enddo
      endif

c GENERATIONS= ?
      if (maxgen.eq.0) then
       call write2log("Couldn't find GENERATIONS= entry in input file")
       status = 1 
      endif

c POPULATION= ?
      if (banksize.eq.0) then
       call write2log("Can't find POPULATION= entry in input file")
       status=1
      elseif(mod(banksize,2).eq.1) then
       call write2log("POPULATION has to be a even number.")
       status=1
      endif

c WHAMTEMPLATE= ?
      if (ntwham.ne.npdb) then
       call write2log("Number of WHAM templates and PDB files is not equ&
     &al!")
       status=1
      endif

c MREMDTEMPLATE= ?
      if (ntmremd.ne.npdb) then
       call write2log("Number of MREMD templates and PDB files is not eq&
     &ual")
       status=1
      endif

c SCRIPTS= ?
      if (nscripts.gt.0) then
       do i=1,nscripts
        inquire(FILE=trim(scripts(i)),EXIST=ex)
        if (.not.ex) then
         call write2log("Can't find file "//trim(scripts(i)))
         status = 1
        endif
       enddo
      endif


      end subroutine ReadInput


c ======================================================================
c  CreateInputs subroutine
c ======================================================================
c  Creates input files for MREMD, WHAM and ZSCORE. 
c ----------------------------------------------------------------------

      subroutine CreateInputs(rozmiar,pop)
       include 'io.inc'
       include 'GA.inc'
       include 'common.inc'
       character(80),dimension(5) :: wagi
       integer*4 :: rozmiar,i,j,k
       real*8 :: pop
       dimension pop(rozmiar,21)
       !character*50 :: pdbfiles(maxnpdb)
       character*50 :: command,tmptext, prefix
       character*256 :: line,plik
       integer :: gdzie = 0
       character(len=*), parameter :: FMT = "(F9.5,F8.5,F8.5,F8.5,F8.5,F&
     &8.5,F8.5,F8.5,F8.5,F8.5,F8.5,F8.5,F8.5,F8.5,F8.5,F8.5,F8.5,F8.5,F8&
     &.5,F20.5)"
       !write(tmptext,'(I10)') banksize
       !call write2log(trim(tmptext))

        do j=1,npdb
         prefix=pdbfiles(j)(1:len_trim(pdbfiles(j))-4)
c Create MREMD input
         open(imremd, file = mremdtemplate(j))
         plik=trim(prefix)//"_"//omremdfn
         call write2log(trim("Creating MREMD input: "//plik))         
         open(omremd, file = plik)
3000     read(imremd, '(A)', end = 3010) line
         gdzie = index (line, '{PREFIX}')
         if (gdzie .ne. 0) then
          line=trim(line(1:gdzie-1))//trim(prefix)//trim(line(gdzie+8:))
          write(omremd,"(A)") trim(line)
         elseif (index(line,'{POPSIZE}') .ne. 0) then
          gdzie = index (line, '{POPSIZE}')
          write(tmptext,'(I4)') rozmiar
          line=trim(line(1:gdzie-1))//trim(adjustl(tmptext))//trim(line(
     &gdzie+9:))
          write(omremd,"(A)") trim(line)
         elseif (index(line,'{WEIGHTS}') .ne. 0) then
          gdzie = index (line, '{WEIGHTS}')
          do i=1,rozmiar
           call Weights2Str(pop(i,:),CUTOFF,wagi)
           write(omremd,"(A)") wagi(1)
           write(omremd,"(A)") wagi(2)
           write(omremd,"(A)") wagi(3)
           write(omremd,"(A)") wagi(4)
           write(omremd,"(A)") wagi(5)
          enddo
         else
          write(omremd,"(A)") trim(line)
         endif

         goto 3000
3010     continue
         close(imremd)
         close(omremd)

c Create WHAM input
         open(iwham, file = whamtemplate(j))
         plik=trim(prefix)//"_"//owhamfn
         call write2log(trim("Creating WHAM input: "//plik))         
         open(owham, file = plik)
3015     read(iwham, '(A)', end = 3020) line
         gdzie = index (line, '{PREFIX}')
         if (gdzie .ne. 0) then
          line=trim(line(1:gdzie-1))//trim(prefix)//trim(line(gdzie+8:))
          write(owham,"(A)") trim(line)
         elseif (index(line,'{POPSIZE}') .ne. 0) then
          gdzie = index (line, '{POPSIZE}')
          write(tmptext,'(I4)') rozmiar
          line=trim(line(1:gdzie-1))//trim(adjustl(tmptext))//trim(line(&
     &gdzie+9:))
          write(owham,"(A)") trim(line)
         elseif (index(line,'{WEIGHTS}') .ne. 0) then
          gdzie = index (line, '{WEIGHTS}')
          do i=1,rozmiar
           call Weights2Str(pop(i,:),CUTOFF,wagi)
           write(owham,"(A)") wagi(1)
           write(owham,"(A)") wagi(2)
           write(owham,"(A)") wagi(3)
           write(owham,"(A)") wagi(4)
           write(owham,"(A)") wagi(5)
           write(owham,"(A)") ''
          enddo
         else
           write(owham,"(A)") trim(line)
         endif

         goto 3015
3020     continue
         close(iwham)
         close(owham)
         enddo
c Create ZSCORE input
         open(izscore, file = izscorefn)
         call write2log(trim("Creating ZSCORE input: "//ozscorefn))
         open(ozscore, file = ozscorefn)
3025     read(izscore, '(A)', end = 3030) line
         gdzie = index (line, '{PREFIX')
         if (gdzie .ne. 0) then
          read(line(gdzie+7:gdzie+9),'(I2.2)') k 
          prefix=pdbfiles(k)(1:len_trim(pdbfiles(k))-4)  
          line=trim(line(1:gdzie-1))//trim(prefix)//trim(line(gdzie+10:)&
     &)

          write(ozscore,"(A)") trim(line)
         elseif (index(line,'{PARAM') .ne. 0) then
          gdzie=index(line,'{PARAM')
          read(line(gdzie+6:gdzie+8),'(I2.2)') k 
          prefix=pdbfiles(k)(1:len_trim(pdbfiles(k))-4)
          do i=1,rozmiar
           write(tmptext,'(I3.3)') i
           line=trim(line(1:gdzie-1))//trim(prefix)//"_wham_par"//tmptex&
     &t
           write(ozscore,"(A)") trim(line)
          enddo
         elseif (index(line,'{POPSIZE}') .ne. 0) then
          gdzie = index (line, '{POPSIZE}')
          write(tmptext,'(I4)') rozmiar
          line=trim(line(1:gdzie-1))//trim(adjustl(tmptext))//trim(line(&
     &gdzie+9:))
          write(ozscore,"(A)") trim(line)
         elseif (index(line,'{WEIGHTS}') .ne. 0) then
          gdzie = index (line, '{WEIGHTS}')
          do i=1,rozmiar
           call Weights2Str(pop(i,:),CUTOFF,wagi)
           write(ozscore,"(A)") wagi(1)
           write(ozscore,"(A)") wagi(2)
           write(ozscore,"(A)") wagi(3)
           write(ozscore,"(A)") wagi(4)
           write(ozscore,"(A)") trim(wagi(5))
           write(ozscore,"(A)") ''
          enddo
         elseif (index(line,'{MAXMIN}') .ne. 0) then
          gdzie = index (line, '{MAXMIN}')
          write(tmptext,'(I4)') maxminstep
          line=trim(line(1:gdzie-1))//trim(adjustl(tmptext))//trim(line(&
     &gdzie+9:))
          write(ozscore,"(A)") trim(line)
         else
          write(ozscore,"(A)") trim(line)
         endif

         goto 3025
3030     continue
         close(izscore)
         close(ozscore)

c Copy scripts
         call write2log('Copying scripts is disabled.')
c	 do i=1,nscripts
c	  ! Open output file
c	  plik=trim(prefix)//"/"//trim(scripts(i))
c	  open(oscript, file = plik)
c         ! Open input file
c	  plik=trim(scripts(i))
c	  open(iscript, file = plik)
c          ! rewrite from input to output
c3030	  read(iscript, '(A)', end = 3035) line
c	  write(oscript,'(A)') trim(line)
c	  goto 3030
c3035	  continue
c	  close(oscript)
c	  close(iscript)
c	 enddo
c        enddo

c ==========================================
c  Dump weights
c ==========================================

       
c       write(ow,'(I4)') generation
c       write(ow,'(200A)') "# WLONG   WSCP    WELEC   WBOND    WANG   WSC&
c     &LOC  WTOR    WTORD   WCORRH  WCORR4  WCORR5  WCORR6  WEL_LOC WTURN&
c     &3  WTURN4  WTURN6  WVDWPP  WHPB    WSCCOR"
c	do I=1,rozmiar
c         write(ow,FMT) pop(i,1),pop(i,2),pop(i,3),pop(i,4), pop(i,5),pop&
c     &(i,6),pop(i,7),pop(i,8),pop(i,9),pop(i,10),pop(i,11),pop(i,12),pop&
c     &(i,13),pop(i,14),pop(i,15),pop(i,16),pop(i,17),pop(i,18),pop(i,19)
c        enddo
c       close(ow)

      end subroutine CreateInputs

  
c ======================================================================
c  ZnajdzPodobnego function
c ======================================================================
c  Finds the most similar individual to osobnik in population(populacja)
c  and returns his position in the population
c ----------------------------------------------------------------------

      integer*4 function ZnajdzPodobnego(rozmiar,populacja,osobnik,csacu&
     &toff)
      integer*4 :: rozmiar,pozycja
      real*8 :: delta,mindelta,csacutoff
      character*30 :: tmp
      real*8, dimension(rozmiar,21) :: populacja
      real*8, dimension(21) :: osobnik 

      mindelta = csacutoff
c mindelta = 10000.0
      pozycja = 0

      do i=1,rozmiar
        delta = 0.0
        do j=1,19
          delta=delta+(populacja(i,j)-osobnik(j))**2
        enddo
c        write(tmp,'(2F7.5 )') delta,csacutoff 
c        call write2log("Delta is "//tmp)
        
        if (delta.lt.mindelta) then
          mindelta=delta
          pozycja=i
          write(tmp,'(F7.5,X,I4)') delta,pozycja
          call write2log("  Delta is "//tmp)
        endif
      end do

      ZnajdzPodobnego = pozycja
      end function
  
c ======================================================================
c  WybierzOsobnika function
c ======================================================================
c  Selects individual from population using roulette method
c ----------------------------------------------------------------------

      integer*4 function WybierzOsobnika(rozmiar,populacja)
      integer*4 :: osobnik, rozmiar 
      real*8 :: los, partsum

      real*8, dimension(rozmiar,21) :: populacja
       
      partsum = 0.0
      osobnik = 0
      los = rand(0)
       do
         osobnik = osobnik + 1
         partsum = partsum + populacja(osobnik,21)
         if ((osobnik.eq.rozmiar).or.(partsum.ge.los)) exit
       end do
      
       WybierzOsobnika = osobnik
       return
      end function

c ======================================================================
c  Najlepszy function
c ======================================================================
c  Selects the individual with the best score from population
c ----------------------------------------------------------------------

      integer*4 function Najlepszy(rozmiar, maks, populacja)
      integer*4 :: rozmiar
      real*8 :: maks
      real*8, dimension(rozmiar,21) :: populacja
       Najlepszy = 0
      do i=1,rozmiar
       if (populacja(i,20).eq.maks) then
        Najlepszy = i
       endif
      end do
       return
      end function

c ======================================================================
c  WriteState subroutine
c ======================================================================
c  Writes current state of calculations in case jobs have to be
c  restartet
c ----------------------------------------------------------------------

      subroutine WriteState()
       include 'GA.inc'
       include 'io.inc'
       include 'common.inc'

       open(ostate,file=ostatefn)
       write(ostate,'(I4)') generation
       write(ostate,'(L2)') do_optima
       write(ostate,'(L2)') do_ga
       write(ostate,'(F7.5)') avrd
       close(ostate)
      end subroutine WriteState
 
c ======================================================================
c  ReadPopulation subroutine
c ======================================================================
c  Reads population from population.summary file
c ----------------------------------------------------------------------

      subroutine ReadPopulation(nind,pop)
       integer*4 :: i,nind
       real*8 :: pop
       character*8 :: dummy
       dimension pop(nind,21) 
       include 'GA.inc'
       include 'io.inc'
       character(len=*), parameter :: FMT = "(F9.5,F8.5,F8.5,F8.5,F8.5,F&
     &8.5,F8.5,F8.5,F8.5,F8.5,F8.5,F8.5,F8.5,F8.5,F8.5,F8.5,F8.5,F8.5,F8&
     &.5,F20.5)"
       read(opopsum,'(A)') dummy
       do i=1,nind
       read(opopsum,FMT) pop(i,1),pop(i,2),pop(i,3),pop(i,4), pop(i,5),p&
     &op(i,6),pop(i,7),pop(i,8),pop(i,9),pop(i,10),pop(i,11),pop(i,12),p&
     &op(i,13),pop(i,14),pop(i,15),pop(i,16),pop(i,17),pop(i,18),pop(i,1&
     &9),pop(i,20)
       end do
      end subroutine ReadPopulation

c ======================================================================
c  ReadOptimaW subroutine
c ======================================================================
c  Reads weights from zscore generated weights_opt.zscorepar_XXX files
c  
c  nind - number of individuals/files to read
c  pop - array where weights are stored  
c ----------------------------------------------------------------------

      subroutine ReadOptimaW(nind,pop)
       integer*4 :: i,nind,stat

       real*8 :: pop
       character*40 :: filename
       character*200 :: tmptxt
       dimension pop(nind,21)
       logical :: err 
       include 'GA.inc'
       include 'io.inc'
   
c Cosmetic stuff for pretty log
       call write2log("==[ Reading weights from zscore files ]=======")
       call write2log("# WLONG    WSCP     WELEC    WBOND    WANG     WS&
     &CLOC   WTOR     WTORD    WCORRH   WCORR4   WCORR5   WCORR6   WEL_L&
     &OC  WTURN3   WTURN4   WTURN6   WVDWPP   WHPB     WSCCOR")

c Actual weight read and error handling
       do i=1,nind
        err=.false.
        write(tmptxt,'(I3.3)') i
        filename = "weights_opt.zscorepar_"//trim(tmptxt)
        open(izsoptw, file = filename)
        read(izsoptw,"(A4,F7.5,A6,F7.5,A7,F7.5,A7,F7.5,A6,F7.5,A)",iosta&
     &t=stat) tmptxt,pop(i,1),tmptxt,pop(i,2),tmptxt,pop(i,3),tmptxt,pop&
     &(i,4),tmptxt,pop(i,5)
        if (stat/=0) then
         err=.true.
         call write2log("-- ERROR while reading "//trim(filename)//"--")
         do j=1,21 
          pop(i,j)=0.0
         enddo
         goto 4000
        endif
        read(izsoptw,"(A7,F7.5,A6,F7.5,A7,F7.5,A8,F7.5,A8,F7.5,A)",iosta&
     &t=stat) tmptxt,pop(i,6),tmptxt,pop(i,7),tmptxt,pop(i,8),tmptxt,pop&
     &(i,9),tmptxt,pop(i,11)
        if (stat/=0) then
         err=.true.
         call write2log("-- ERROR while reading "//trim(filename)//"--")
         do j=1,21 
          pop(i,j)=0.0
         enddo
         goto 4000
        endif
        read(izsoptw,"(A7,F7.5,A9,F7.5,A8,F7.5,A8,F7.5,A8,F7.5,A)",iosta&
     &t=stat) tmptxt,pop(i,12),tmptxt,pop(i,13),tmptxt,pop(i,14),tmptxt,&
     &pop(i,15),tmptxt,pop(i,16)
         if (stat/=0) then
         err=.true.
         call write2log("-- ERROR while reading "//trim(filename)//"--")
         do j=1,21 
          pop(i,j)=0.0
         enddo
         goto 4000
        endif
        read(izsoptw,"(A7,F7.5,A8,F7.5,A6,F7.5,A)",iostat=stat) tmptxt, &
     &pop(i,19),tmptxt,pop(i,17),tmptxt,pop(i,18),tmptxt
         if (stat/=0) then
         err=.true.
         call write2log("-- ERROR while reading "//trim(filename)//"--")
         do j=1,21 
          pop(i,j)=0.0
         enddo
         goto 4000
        endif
        read(izsoptw,"(A22,F7.5)",iostat=stat) tmptxt, pop(i,10)
        if (stat/=0) then
        err=.true.
        call write2log("-- ERROR while reading "//trim(filename)//"--")
        do j=1,21 
         pop(i,j)=0.0
        enddo
        endif

        if (.not.err) then
         write(tmptxt,'(19F9.5)') (pop(i,j),j=1,19)
         call write2log(trim(tmptxt))
        endif 

4000    close(izsoptw)
       end do
      end subroutine ReadOptimaW


c ======================================================================
c  ReadZEnergy subroutine
c ======================================================================
c  Reads Energy from ZSCORE output files zscore_out_parXXX_GB000
c ----------------------------------------------------------------------

      subroutine ReadZEnergy(nind,pop)
       integer*4 :: i,nind
       integer*4 :: stat
       real*8 :: pop
       character*40 :: filename,tmptxt,tmptxt2
       character*100 :: wiersz
       dimension pop(nind,21) 
       include 'io.inc'
    
       call write2log("==[ Reading Final Function Value from zscore ]===&
     &====")
       do i=1,nind
        write(tmptxt,'(I3.3)') i
        pop(i,20)=0
        filename=trim(izenergyfn)//trim(tmptxt)//"_GB000"
        open(izenergy, file = filename)
        do
         read(izenergy, '(A)', iostat=stat) wiersz
         if (stat /= 0) exit
         if (wiersz(1:22).eq.' Final function value:') then
          tmptxt=trim(adjustl(wiersz(23:48)))
          read(tmptxt,'(F16.5)') pop(i,20)
          write(tmptxt2,'(E15.7)') pop(i,20) 
          call write2log("Reading score from "//filename//" : "//tmptxt/&
     &/" : "//tmptxt2)
         endif
        end do
        if (pop(i,20).eq.0) then
c Setting FFV to a realy big valiue or otherwise CalcFittness will have problems
         pop(i,20)=huge(0.0d0)
         call write2log("ERROR while reading FFV from "//filename)
        endif 
        close(izenergy)
       end do

      end subroutine ReadZEnergy

c =======================================================================
c  WritePopSum subroutine
c =======================================================================
c  Writes whole individuals (weights, score, ID, generation) to file 
c  population.summary  This file can be used to plot evolution of the 
c  weights and scores over time. 
c -----------------------------------------------------------------------
      subroutine WritePopSum()
       logical :: jestplik = .false.
       character*5 :: tmptext
c       integer*4 :: nind
c       real*8 :: bank 
c       dimension bank(nind,21)
       include 'io.inc'
       include 'common.inc'
       

c      write(*,*)  banksize

      write(tmptext,'(I3.3)') generation
      inquire(file=opopsumfn, exist=jestplik) 
      if (.not. jestplik) then
       open(opopsum, file = opopsumfn)
       write(opopsum, "(A189)") "# WLONG   WSCP    WELEC   WBOND   WANG &
     &   WSCLOC  WTOR    WTORD   WCORRH  WCORR4  WCORR5  WCORR6  WEL_LOC&
     & WTURN3  WTURN4  WTURN6  WVDWPP  WHPB    WSCCOR                SCO&
     &RE    IND      GEN"
       close(opopsum)
      endif
   
      open(opopsum, file = opopsumfn, access = "append")
c       call ReadPopulation(banksize,populacja)
      write(opopsum, "(A189)") "# ==== Generation "//trim(tmptext)//" ==&
     &==================================================================&
     &==================================================================&
     &================================="
c       call write2log(banksize)
       do i=1,banksize
c      write(*,*) i 
c        write(opopsum,"(F9.5,2F8.5)") populacja(i,1),populacja(i,2),popu&
c     &lacja(i,3) 
      write(opopsum,"(F9.5,F8.5,F8.5,F8.5,F8.5,F8.5,F8.5,F8.5,F8.5,F8.5,&
     &F8.5,F8.5,F8.5,F8.5,F8.5,F8.5,F8.5,F8.5,F8.5,F20.5,I8,I8)") bank(i&
     &,1),bank(i,2),bank(i,3),bank(i,4),bank(i,5),bank(i,6),bank(i,7),ba&
     &nk(i,8),bank(i,9),bank(i,10),bank(i,11),bank(i,12),bank(i,13),bank&
     &(i,14),bank(i,15),bank(i,16),bank(i,17),bank(i,18),bank(i,19),bank&
     &(i,20),i,iter
       enddo
      write(opopsum,*) ""
      close(opopsum)

      end subroutine WritePopSum

c ======================================================================
c  WriteBank subroutine
c ======================================================================
c  Writes CSA BANK to file 
c ----------------------------------------------------------------------
      subroutine WriteBank(b)
       include 'io.inc'
       include 'common.inc'
       real*8,dimension(banksize,21) :: b  
       character*250 :: tmptext
       character*250 :: header   


       header="#    WLONG      WSCP     WELEC     WBOND      WANG    WSC&
     &LOC      WTOR     WTORD    WCORRH    WCORR4    WCORR5    WCORR6   &
     &WEL_LOC    WTURN3    WTURN4    WTURN6    WVDWPP      WHPB    WSCCO&
     &R            FFV   FITNESS" 
 
       call write2log("Writing Bank to file "//iobankfn)
       open(iobank, file = iobankfn)
       write(iobank, "(A)") trim(adjustl(header))
       
       do i=1,banksize
        write(iobank,'(19F10.5,E15.7,F10.5)') (b(i,j),j=1,21)
        write(tmptext,'(19F10.5,E15.7,F10.5)') (b(i,j),j=1,21)
        call write2log(trim(tmptext))  
       enddo 
       close(iobank)
      end subroutine

c ======================================================================
c  WriteBankHistory subroutine
c ======================================================================
c  Writes CSA BANK History to file 
c ----------------------------------------------------------------------
      subroutine WriteBankHistory(b)
       include 'io.inc'
       include 'common.inc'
       real*8,dimension(banksize,21) :: b  
       character*250 :: tmptext
       character*250 :: header   


       header="#    WLONG      WSCP     WELEC     WBOND      WANG    WSC&
     &LOC      WTOR     WTORD    WCORRH    WCORR4    WCORR5    WCORR6   &
     &WEL_LOC    WTURN3    WTURN4    WTURN6    WVDWPP      WHPB    WSCCO&
     &R            FFV   FITNESS" 
 
       call write2log("Writing Bank history to file "//obankhfn)
       open(obankh, file = obankhfn)
       write(iobank, "(A)") trim(adjustl(header))
       
       do i=1,banksize
        write(obankh,'(19F10.5,E15.7,F10.5)') (b(i,j),j=1,21)
        write(tmptext,'(19F10.5,E15.7,F10.5)') (b(i,j),j=1,21)
        call write2log(trim(tmptext))  
       enddo 
       close(iobank)
      end subroutine

c
c ======================================================================
c  ReadBank subroutine
c ======================================================================
c  Read CSA BANK from file
c ----------------------------------------------------------------------
      subroutine ReadBank(b)
       include 'io.inc'
       include 'common.inc'
       real*8,dimension(banksize,21) :: b 
       character*250 :: tmptext 

       call write2log("Reading Bank from file "//iobankfn)  
       open(iobank, file = iobankfn)
c skip header
       read(iobank,"(A)") tmptext 
       do i=1,banksize
        read(iobank,"(19F10.5,E15.7,F10.5)") (b(i,j),j=1,21)
        write(tmptext,'(19F10.5,E15.7,F10.5)') (b(i,j),j=1,21)
        call write2log(trim(tmptext))
       enddo 
       close(iobank)
      end subroutine

c ======================================================================
c  CalcFitness subroutine
c ======================================================================
c  Calculates fitness of every individual in the bank
c ----------------------------------------------------------------------
      subroutine CalcFitness(b,fitn)
       include 'common.inc'
       real*8,dimension(banksize,21) :: b
       real*8 :: fitn, sumfitn

       fitn = 0.0
       sumfitn = 0.0

       do i=1,banksize
        fitn=fitn+(1/b(i,20))
        sumfitn=sumfitn+b(i,20)
       end do
       do i=1,banksize
        b(i,21)=1/(b(i,20)*fitn)
c        b(i,21)=b(i,20)/fitn
       end do
       fitn=sumfitn 
      end subroutine

c ======================================================================
c  CalcAvgDist subroutine
c ======================================================================
c  Calculates average distance between individuals in the bank
c ----------------------------------------------------------------------
      subroutine CalcAvgDist(b,avgd)
       include 'common.inc'
       real*8,dimension(banksize,21) :: b
       real*8 :: d,avgd
       integer*4 :: nd,w

       d=0.0                          ! distance
       nd = (banksize-1)*banksize/2   ! number of distances to calculate

       do i=1,banksize-1
        do j=i+1,banksize
         do w=1,19
          d=d+(b(i,w)-b(j,w))**2
         end do
        end do 
       end do
       avgd=sqrt(d)/nd
      end subroutine

c -----------------------------------------------------------------------